Kibble-Zurek mechanism in a quenched ferromagnetic Bose-Einstein condensate
Hiroki Saito1, Yuki Kawaguchi2, and Masahito Ueda2,3
Department of Applied Physics and Chemistry, The University of Electro-Communications, Tokyo 182-8585, Japan
Department of Physics, Tokyo Institute of Technology, Tokyo 152-8551, Japan
Macroscopic Quantum Control Project, ERATO, JST, Bunkyo-ku, Tokyo 113-8656, Japan
(Dated: October 31, 2018)
The spin vortices are shown to be created through the Kibble-Zurek (KZ) mechanism in a quantum
phase transition of a spin-1 ferromagnetic Bose-Einstein condensate, when the applied magnetic field
is quenched below a critical value. It is shown that the magnetic correlation functions have finite
correlation lengths, and magnetizations at widely separated positions grow in random directions,
resulting in spin vortices. We numerically confirm the scaling law that the winding number of spin
vortices is proportional to the square root of the length of the closed path, and for slow quench,
proportional to τ
with τQ being the quench time. The relation between the spin conservation
and the KZ mechanism is discussed.
PACS numbers: 03.75.Mn, 03.75.Lm, 73.43.Nq, 64.60.Ht
I. INTRODUCTION
Spontaneous symmetry breaking in a phase transition
produces local domains of an order parameter. If do-
mains are separated by such a long distance that they
cannot exchange information, local domains grow ini-
tially with random phases and eventually give rise to
topological defects when they overlap. This scenario
of topological-defect formation in continuous-symmetry
breaking is known as the Kibble-Zurek (KZ) mecha-
nism [1, 2], which originally predicted the cosmic-string
and monopole formation in the early Universe [1], and
has since been applied to a wide variety of systems. Ex-
perimentally, the KZ mechanism has been examined in
liquid crystals [3, 4], superfluid 4He [5] and 3He [6, 7],
an optical Kerr medium [8], Josephson junctions [9, 10],
and superconducting films [11].
Recently, spontaneous magnetization in a spinor Bose-
Einstein condensate (BEC) has attracted much interest
as a new system for studying the KZ mechanism [12, 13,
14, 15]. In the experiment performed by the Berkeley
group [12], a BEC of F = 1 87Rb atoms are prepared in
the m = 0 state, where F is the hyperfine spin and m
is its projection on the direction of the magnetic field.
By quench of the magnetic field, say in the z direction,
magnetization appears in the x-y plane. Since the spinor
Hamiltonian is axisymmetric with respect to the z axis,
the magnetization in the x-y direction breaks the U(1)
symmetry in the spin space. Thus, local domain forma-
tion is expected to lead to topological defects — spin
vortices — through the KZ mechanism.
However, the origin of the spin vortices observed after
the quench in the Berkeley experiment [12] cannot be at-
tributed to the KZ mechanism. In fact, in Ref. [12], the
spin correlation extends over the entire system (at least
in the x direction) and the domains are not independent
with each other. We have shown in Ref. [13] that the ori-
gin of the observed spin vortices is initial spin correlation
due to the residual m = ±1 atoms, which forms domain
structure followed by spin vortex creation [16]. In order
to realize the KZ mechanism in this system, i.e., in order
to ensure that the magnetic domains grow independently,
the size of the system must be much larger than the spin
correlation length and the long-range correlation in the
initial spin state must be absent. The aim of the present
paper is to show that under these conditions spin vortices
are generated through the KZ mechanism.
In the present paper we will consider 1D-ring and 2D-
disk geometries. We will show that in the 1D ring the
average spin winding number after the quench is propor-
tional to the square root of the system size, which is in
agreement with the KZ prediction [2]. In 2D the winding
number along a path with radius R is also proportional
to R1/2 as long as R is much larger than the vortex spac-
ing, while it is proportional to R for small R. When
the magnetic field is quenched slowly, the winding num-
ber is shown to be proportional to τ
Q with τQ being
the quench time. This power law can be understood by
Zurek’s simple discussion [2].
The spinor BEC is different from the other systems in
which the KZ mechanism has been observed, in that the
total spin is conserved when the quadratic Zeeman en-
ergy q is negligible. This fact is seemingly incompatible
with the KZ postulate, since the magnetic domains must
be correlated with each other so that the total magnetiza-
tion vanishes. We will show that for q = 0 small magnetic
domains are aligned to cancel out the local spin averaged
over the correlation length, and that they are indepen-
dent with each other over a greater length scale; the spin
conservation is thus compatible with the KZ mechanism.
The present paper is organized as follows. Section II
analyzes spontaneous magnetization of a spin-1 BEC and
the resultant magnetic correlation functions using the
Bogoliubov approximation. Section III performs numeri-
cal simulations of the dynamics of quenched BECs in 1D
and 2D, and shows that the KZ mechanism does emerge
in the present system. Section IV provides conclusions.
http://arxiv.org/abs/0704.1377v1
II. BOGOLIUBOV ANALYSIS OF A
QUENCHED FERROMAGNETIC
BOSE-EINSTEIN CONDENSATE
A. Hamiltonian for the spin-1 atoms
We consider spin-1 bosonic atoms with mass M con-
fined in a potential Vtrap(r). The noninteracting part of
the Hamiltonian is given by
Ĥ0 =
ψ̂†m(r)
2 + Vtrap(r)
ψ̂m(r),
where ψ̂m(r) annihilates an atom in magnetic sublevel m
of spin at a position r.
The interaction between atoms with s-wave scattering
is described by
Ĥint =
2(r) + c1F̂
:, (2)
where the symbol :: denotes the normal order and
ρ̂(r) =
ψ̂†m(r)ψ̂m(r), (3)
F̂ (r) =
ψ̂†m(r)fmm′ ψ̂m′(r), (4)
with f = (fx, fy, fz) being the spin-1 matrices. The in-
teraction coefficients in Eq. (2) are given by
4πh̄2
a0 + 2a2
, (5a)
4πh̄2
a2 − a0
, (5b)
where aS is the s-wave scattering lengths for two colliding
atoms with total spin S.
When magnetic field B is applied, the linear Zeeman
effect rotates the spin around the direction of B at the
Larmor frequency. Since Ĥ0 and Ĥint are spin-rotation
invariant and we assume the uniform magnetic field, the
linear Zeeman term has only a trivial effect on spin dy-
namics — uniform rotation of spins about B — which is
therefore ignored. The quadratic Zeeman effects for an
F = 1 87Rb atom is described by
Ĥq =
ψ̂†m(r)
(B · f)2
ψ̂m′(r), (6)
where µB is the Bohr magneton and Ehf > 0 is the hy-
perfine splitting energy between F = 1 and F = 2. The
total Hamiltonian is given by the sum of Eqs. (1), (2),
and (6),
Ĥ = Ĥ0 + Ĥq + Ĥint. (7)
B. Time evolution in the Bogoliubov
approximation
In the initial state, all atoms are prepared in them = 0
state. We study the spin dynamics of the system using
the Bogoliubov approximation with respect to this initial
state. For simplicity, we assume Vtrap = 0 in this section.
In the Bogoliubov approximation, the BEC part in the
field operator is replaced by a c-number. In the present
case, we write the m = 0 component of the field operator
ψ̂0(r) = e
−ic0nt/h̄
n+ δψ̂0(r)
, (8)
where n is the atomic density. We expand ψ̂±1(r) as
ψ̂±1(r) = e
−ic0nt/h̄
eik·r â±1,k, (9)
where V is the volume of the system and â±1,k is the
annihilation operator of an atom in the m = ±1 state
with wave vector k. Keeping only up to the second order
of δψ̂0(r) and ψ̂±1(r) in the Hamiltonian, we obtain the
Heisenberg equation of motion for â±1,k as
dâ±1,k(t)
= (εk + q + c1n)â±1,k(t) + c1nâ
∓1,−k(t),
where εk = h̄
2k2/(2M) and q = µ2BB
2/(4Ehf). The mag-
netic field is assumed to be applied in the z direction. The
solution of Eq. (10) is obtained as
â±1,k(t) =
− i εk + q + c1n
â±1,k(0)
∓1,−k(0), (11)
where
(εk + q)(εk + q + 2c1n). (12)
When Ek is imaginary, the corresponding modes are
dynamically unstable and grow exponentially. Since c1 <
0 and q > 0 for F = 1 87Rb atoms, the exponential
growth occurs for
q < 2|c1|n ≡ qc. (13)
This critical value of q agrees with the phase boundary
between the polar phase and the broken-axisymmetry
phase [17, 18]. When q ≤ qc/2, the wave number of
the most unstable mode is
kmu = ±
, (14)
and when qc/2 < q < qc, kmu = 0.
C. Fast quench
We consider the situation in which q is much larger
than the other characteristic energies for t < 0, and q is
suddenly quenched below qc at t = 0. During t < 0, the
time evolution in Eq. (11) is â±1,k(t) ≃ e−iqt/h̄â±1,k(0),
and the m = ±1 state remains in the vacuum state. For
t > 0, we obtain the time evolution of the density of the
m = ±1 component as
±1(r, t)ψ̂±1(r, t)
16|Ek|2
e2|Ek|t/h̄,(15)
where the expectation value is taken with respect to the
vacuum state of the m = ±1 component. In the sec-
ond line of Eq. (15), we have kept the unstable modes
alone with k < kc ≡
2M(qc − q)/h̄ by assuming that
|Ek|t/h̄ ≫ 1. This result indicates that the m = ±1
components grow exponentially after the quench.
Since the operator ψ̂0 in Eq. (4) is replaced by
the Bogoliubov approximation, the magnetization oper-
ator F̂+ = F̂
− = F̂x + iF̂y has the form,
F̂+(r) =
1(r) + ψ̂−1(r)
. (16)
Using Eq. (11), the time evolution of the correlation func-
tion is calculated to be
F̂+(r, t)F̂−(r
′, t)
εk + q
eik·(r−r
(17a)
qc − q − εk
e2|Ek|t/h̄+ik·(r−r
′), (17b)
where in the second line we have kept the unstable modes
alone.
From the exponential factor in Eq. (17b), we see that
the sum is contributed mostly from k around the mode
with maximum |Ek|. The denominator in the summand
of Eq. (17b) is much smoother than the exponential fac-
tor if q is not close to qc, and then we approximate εk
with εmu = h̄
2k2mu/(2M) in the denominator. We expand
2|Ek|t/h̄ around kmu in the exponent as
2|Ek|t
ξ2∆k2 − 1
Ξ4∆k4
+O(∆k6),
where ∆k = k − kmu. It is clear that τ sets the time
scale for the exponential growth. The magnetization is
observed when it sufficiently grows, i.e., t ∼ τ . Replacing
the summation with the Gaussian integral in Eq. (17b),
we find that ξ represents the correlation length. For q <
qc/2, kmu is given by Eq. (14), and
, (19)
qc − 2q
. (20)
For qc/2 < q < qc, kmu = 0 and
q(qc − q)
, (21)
2q − qc
q(qc − q)
. (22)
At q = qc/2, Eqs. (20) and (22) vanish, and the ∆k
4 term
in Eq. (18) becomes important, with
Ξ = 4
2M2q2c
. (23)
We first consider a 1D system with the periodic bound-
ary condition, i.e., the 1D ring geometry. We assume that
the radius of the ring R is much larger than the domain
size, and the curvature of the ring does not affect the
dynamics.
For q < qc/2, the magnetic correlation function is cal-
culated to be
F̂+(θ, t)F̂−(θ
′, t)
cos[kmuR(θ − θ′)]
×et/τ−τR
2(θ−θ′)2/(tξ2), (24)
where τ and ξ are given by Eqs. (19) and (20), and θ and
θ′ are azimuthal angles. For qc/2 < q < qc, we obtain
F̂+(θ, t)F̂−(θ
′, t)
qc − q
et/τ−τR
2(θ−θ′)2/(tξ2)
with Eqs. (21) and (22). At q = qc/2, the correlation
function reads
F̂+(θ, t)F̂−(θ
′, t)
qc − q
τR4(θ − θ′)4
R2(θ − θ′)2
τR4(θ − θ′)4
,(26)
where Γ is the Gamma function and
0F2(a, b, z) =
Γ(a)Γ(b)
Γ(a+ j)Γ(b + j)
is the generalized hypergeometric function. Equa-
tion (26) is shown in Fig. 2(a), where Ξ gives a char-
acteristic width of the correlation function.
Next, we consider the 2D geometry. For qc/2 < q <
qc, and then kmu = 0, the integral can be performed
analytically, giving
F̂+(r, t)F̂−(r
′, t)
2πξ2t
qc − q
et/τ−τ |r−r
′|2/(tξ2),
where τ and ξ are given in Eqs. (21) and (22). For other
q, we can perform only the angular integral as
F̂+(r, t)F̂−(r
′, t)
qc − q − εmu
kJ0(k|r − r′|)e2|Ek|t/h̄dk,
where J0 is the Bessel function. If the exponential factor
is much sharper than the Bessel function around kmu,
the correlation function (29) is approximated to be ∝
J0(kmu|r − r′|) [14, 15].
As shown above, the correlation function (17b) has a
finite correlation length, and the magnetization at posi-
tions widely separated from each other grow with inde-
pendent directions in the x-y plane. Thus, the growth
of the magnetic domains is expected to leave topological
defects through the KZ mechanism.
D. Slow quench
In the previous sections, we have assumed that the
magnetic field is suddenly quenched to the desired value
at t = 0 and q is held constant for t > 0. We assume here
that for t > 0 the magnetic field is gradually quenched
q(t) = qc
. (30)
The magnetic correlation can be estimated to be
F̂+(r, t)F̂−(r
′, t)
dk exp
2|Ek(t)|t
dt+ ik · (r − r′)
. (31)
Since we are interested in the vicinity of the critical
point where correlation starts to grow, we expand |Ek(t)|
around kmu = 0 and keep the terms up to the order of
k2. For the 1D ring, we obtain
F̂+(θ, t)F̂−(θ
′, t)
∝ ef(t)−R
2(θ−θ′)2/ξ2Q , (32)
and for the 2D geometry,
F̂+(r, t)F̂−(r
′, t)
∝ ef(t)−|r−r
′|2/ξ2
Q , (33)
where
f(t) =
tan−1
τQ − t
1− 2t
, (34)
t(τQ − t)
. (35)
For t≪ τQ, f(t) can be expanded as
f(t) =
, (36)
and from f(t) ∼ 1, the time scale for magnetization is
given by
Q . (37)
Substitution of tQ into Eq. (35) yields
Q . (38)
The same power law is obtained in Ref. [14].
It is interesting to note that the results (37) and
(38) are easily obtained also by the simple discussion by
Zurek [2]. Since q(t) depends on time, τ and ξ given in
Eqs. (21) and (22) are time dependent, and hence they
are regarded as the growth time and correlation length
at each instant of time. The local magnetization is de-
veloped after a time tQ has elapsed such that
τ(tQ) ∼ tQ. (39)
Using
τ(t) =
t(τQ − t)
, (40)
we obtain tQ in Eq. (37). Substituting this tQ into
ξ2(t) =
τQ − 2t
t(τQ − t)
yields Eq. (38).
III. NUMERICAL SIMULATIONS AND THE
KIBBLE-ZUREK MECHANISM
A. Gross-Pitaevskii equation with quantum
fluctuations
The multicomponent Gross-Pitaevskii (GP) equation
is obtained by replacing the field operators ψ̂m with the
macroscopic wave function ψm in the Heisenberg equa-
tion of motion:
∇2 + Vtrap + q + c0ρ
F∓ψ0 ± Fzψ±1
, (42a)
∇2 + Vtrap + c0ρ
(F+ψ1 + F−ψ−1) , (42b)
where ρ and F are defined using ψm instead of ψ̂m in
Eqs. (3) and (4). The wave function is normalized as
|ψm|2 = N, (43)
with N being the number of atoms in the condensate.
Suppose that all atoms are initially in the m = 0 state.
It follows then from Eq. (42a) that ψ±1 will remain zero
in the subsequent time evolution. This is because quan-
tum fluctuations in the transverse magnetization that
trigger the growth of magnetization are neglected in the
mean-field approximation. We therefore introduce an ap-
propriate initial noise in ψ±1 so that the mean-field ap-
proximation reproduces the quantum evolution.
Let us write the initial state as
ψ±1(r) =
eik·ra±1,k(0), (44)
where a±1,k are c-numbers. We assume that the c-
number amplitudes a±1,k(0) are stochastic variables
whose average values vanish,
〈a±1,k(0)〉avg = 0, (45)
where by 〈· · · 〉avg we denote the statistical average over
an appropriate probability distribution. The linear ap-
proximation of the GP equation with respect to a±1,k
gives the same time evolution as Eq. (11), in which the
operators are replaced by the c-numbers. We thus obtain
F+(r, t)F−(r
′, t) =
εk + q
e−ik·(r−r
′)|a1,k(0)|2
+eik·(r−r
′)|a−1,−k(0)|2
. (46)
Comparing Eq. (46) with Eq. (17a), we find that they
have the same form if the variance of a±1,k(0) satisfies
〈|a±1,k(0)|2〉avg =
for all k.
0 100 200 300
t [ms]
q=qc/2
0 π 2π
 |F+| /ρ arg F+
q=qc/ 20
0 π 2π
 |F+| /ρ arg F+
(b) (c)
FIG. 1: (Color online) (a) Time evolution of the auto cor-
relation function given in Eq. (49) for the 1D ring geometry.
(b) Magnitude of the normalized magnetization |F+|/ρ (solid
curve, left axis) and direction of the magnetization argF+
(dashed curve, right axis) at t = 70 ms for q = 0 and (c) for
q = qc/2. The radius of the ring is R = 50 µm, the atomic
density is n = 2.8 × 1014 cm−3, and the number of atoms is
N = 106.
In the following, we will perform numerical simulation
of spontaneous magnetization using the GP equation and
show that the ensuing dynamics exhibits defect formation
similar to the KZ mechanism. As the initial state of the
m = ±1 wave functions, we use Eq. (44) with
a±1,k(0) = αrnd + iβrnd, (48)
where αrnd and βrnd are random variables following the
normal distribution p(x) =
2/π exp(−2x2). Equation
(48) then satisfies Eqs. (45) and (47).
B. 1D ring geometry
Let us first investigate the 1D ring system. Experi-
mentally this geometry can be realized, e.g., by an opti-
cal trap using a Laguerre-Gaussian beam [19]. We reduce
the GP equation (42) to 1D by assuming that the wave
function ψm depends only on the azimuthal angle θ. The
average density of atoms is assumed to be n = 2.8× 1014
cm−3. When the radius of the ring is R = 50 µm and
the radius of the small circle is 2 µm, the total number
of atoms is N ≃ 106.
Figure 1 illustrates a single run of time evolution for
an initial state given by Eqs. (44) and (48). Figure 1
(a) shows time evolution of the auto correlation function
defined by
F̄ (t) =
|F+(θ, t)|2
ρ2(θ, t)
. (49)
For both q = 0 and q = qc/2, the transverse magneti-
zation grows exponentially with a time constant ∼ τ =
h̄/qc ≃ 8 ms. Snapshots of the transverse magnetization
at t = 70 ms are shown in Figs. 1 (b) and 1 (c) for q = 0
and q = qc/2, respectively. We define the spin winding
number as
2i|F+|2
, (50)
which represents the number of rotation of the spin vector
in the x-y plane along the circumference, and of course
w is an integer. The spin winding numbers are w = 7 in
Fig. 1 (b) and w = −1 in Fig. 1 (c).
Figure 2 (a) shows the ensemble average of the nor-
malized correlation function,
〈Fcorr(δθ)〉avg =
dθF+(θ)F−(θ + δθ)
dθρ(θ)ρ(θ + δθ)
, (51)
at t = 70 ms. For q = qc/2, the correlation function has
the characteristic width of ∼ Ξ in Eq. (23), indicating
that the ring is filled with magnetic domains with an av-
erage size of ∼ Ξ. According to the KZ theory, the mag-
netic domains with random directions give rise to the spin
winding, which is estimated to be w ∼ (R/Ξ)1/2. This R
dependence of w is clearly seen in Fig. 2 (b). The ensem-
ble average of the winding number, 〈w〉avg , vanishes due
to the random nature of the initial noise, and the square
root of its variance, 〈w2〉1/2avg, should be regarded as a typ-
ical winding number. The variance is expected to obey
the χ2 distribution with 1000 degrees of freedom, and
hence we show the 95% confidence interval to estimate
the statistical errors in Fig. 2. As shown in the inset of
Fig. 2 (b), the typical winding number changes in time,
since the ferromagnetic energy is converted to the kinetic
energy and the system exhibits complicated dynamics.
The situation is different for q = 0, in which the cor-
relation function oscillates with a Gaussian envelope as
shown in Fig. 2 (a). This form of the correlation function
gives us the answer to the question as to how the KZ
mechanism manifests itself in spin conserving systems.
The finite correlation length for q = 0 indicates that the
spin is conserved not only globally but also locally, that
is, the locally integrated spin over the correlation length
|δr|<∼ξ
F (r + δr)dδr, (52)
is held to be zero for any r. This local spin conservation
is due to formation of staggered domain or helical spin
structures whose periodic length is much smaller than
−π/2 0 π/2
q=qc/2
10 100 1000
R [µm]
q=qc/2
0 100
t [ms]
q=qc/2
FIG. 2: (Color online) (a) Numerically obtained correlation
function given in Eq. (51) at t = 70 ms (solid curves), and
theoretical fits (dashed curves) from Eqs. (24) and (26). Other
parameters are the same as those in Fig. 1. (b) R dependence
of the variance of the spin winding number, where the number
of atoms is related to R as N = 106×R [µm] /50. The dashed
lines are semi-log fits to the numerical data. The inset shows
the time dependence of 〈w2〉avg for R = 50 µm. The data
in (a) and (b) are averages over 1000 runs of simulations for
different initial states produced by random numbers. The
error bars in (b) represent the 95% confidence interval of the
χ2 distribution.
ξ. Thus, the neighboring domains tend to have opposite
magnetizations to cancel out the spin locally, and the
domains far from each other grow independently; the spin
conservation and the KZ mechanism are thus compatible.
The oscillation in the correlation function originates
from the fact that the most unstable modes have nonzero
wave numbers ±kmu. Each correlated region of size ∼
ξ = [8h̄2/(Mqc)]
1/2 contains spin waves of eikmuRθ and
e−ikmuRθ. If there is an imbalance between these modes,
the winding number monotonically increases or decreases
in each region of ∼ ξ. This is the reason why 〈w2〉avg is
larger for q = 0 than for q = qc/2 in Fig. 2 (b). It follows
from this consideration that for kmuξ ≫ 1 the winding
0.1 1
1−2q/qc
slope = 3/2
FIG. 3: (Color online) Dependence of the variance of the
spin winding number on q. Except for q, the parameters are
the same as those in Fig. 1. The dashed line is proportional
to (1−2q/qc)
3/2. The plots show the averages over 1000 runs
of simulations for different initial states produced by random
numbers. The error bars represent the 95% confidence inter-
val of the χ2 distribution.
number is proportional to
w ∼ kmuξ
= kmu
1− 2q
, (53)
where Eqs. (14) and (20) are used. Figure 3 shows the
averaged variance of the winding number versus 1−2q/qc.
For small q, 〈w2〉avg is proportional to (1 − 2q/qc)3/2, in
agreement with Eq. (53). When q is close to qc/2, the
spin winding within the correlated region, kmuξ, becomes
small, and then the winding number reduces to the value
shown in Fig. 2 (b), i.e., 〈w2〉avg ≃ 4.
We next discuss the results of simulations of slow
quench as in Eq. (30). Since the winding number for
the slow quench is small compared with the fast quench,
we take a large ring of R = 400 µm. Figure 4 shows
the variance of the winding number as a function of the
quench time. We can clearly see that 〈w2〉avg has a power
law of τ
Q within the statistical error, which is in agree-
ment with ξ−1Q ∼ τ
Q , with ξQ being given in Eq. (38).
Thus, the present system follow the quench-time scaling
of Zurek [2]. We note that, in the slow quench, the wind-
ing number converges to an almost constant value for
varying quench time τQ, as shown in the inset of Fig. 4.
This is because little excess energy other than for exciting
spin vortices is available for the slow quench.
C. 2D disk geometry
When the confinement in the z direction is tight, the
system is effectively 2D. For simplicity, we ignore the den-
sity dependence in the z direction, and assume that the
τQ[s]
slope=−1/3
0 100 200 300 400
t [ms]
0.8 s
1.6 s
3.2 s
6.4 s
FIG. 4: (Color online) Dependence of the variance of the spin
winding number at t = 400 ms on the quench time τQ, where q
is varied as in Eq. (30). The radius of the ring is R = 400 µm,
the atomic density is n = 2.8 × 1014 cm−3, and the number
of atoms is N = 8 × 106. The dashed line is proportional
. The inset shows time evolution of 〈w2〉avg. The
data are averages over 1000 runs of simulations for different
initial states produced by random numbers. The error bars
represent the 95% confidence interval of the χ2 distribution.
2D GP equation has the same form as Eq. (42). We as-
sume that the wave function vanishes at the wall located
at (x2 + y2)1/2 = Rw = 100 µm, and that the potential
is flat inside of the wall. Then the density n = 2.8× 1014
cm−3 is almost constant except within the healing length
{3/[8πn(a0 + 2a2)]}1/2 ≃ 0.16 µm near the wall. When
the thickness in the z direction is ≃ 1 µm, the number of
atoms is N ≃ 107. Such a system will be realized using
an optical sheet and a hollow laser beam.
The initial state of ψ0 is a stationary solution of the GP
equation, and the initial state of ψ±1 is given by Eq. (44)
with random variables (48). Figure 5 (a) shows time
evolution of the autocorrelation function of the transverse
magnetization,
F̄ (t) =
|F+(r, t)|2
ρ2(r, t)
, (54)
which grows exponentially with the same time constant
as that in Fig. 1, and saturates for t >∼ 100 ms.
Snapshots of |F+(r)| and argF+(r) at t = 100 ms are
shown in Figs. 5 (b) and 5 (c). We see that |F+(r)|
at t >∼ 100 ms contains many holes, around which the
spin direction rotates by 2π. Since this topological spin
structure consists of singly-quantized vortices in the m =
±1 states filled by atoms in the m = 0 state, it is called
the “polar-core vortex.” We can estimate the spin healing
length ξs by equating the kinetic energy h̄
2/(2Mξ2s ) with
0 100 200 300
t [ms]
q=0 q=qc/2
(c) q = q  / 2
(b) q = 0
t = 50 ms t = 100 mst = 75 ms t = 200 ms
t = 50 ms t = 100 mst = 75 ms t = 200 ms
FIG. 5: (Color) (a) Time evolution of the autocorrelation function given in Eq. (54) for the 2D disk geometry. The radius of
the disk is Rw = 100 µm, the atomic density is n = 2.8× 10
14 cm−3, and the number of atoms is N = 107. (b) Profiles of the
magnetization |F+| (upper) and its direction argF+ (lower) for q = 0 and (c) for q = qc/2. The size of each panel is 200 µm
×200 µm.
the energy of magnetization |q − qc|, giving
2M |q − qc|
. (55)
This length scale is ξs ≃ 1.7 µm for q = 0 and ξs ≃ 2.4
µm for q = qc/2, which are in good agreement with the
sizes of the vortex cores in Figs. 5 (b) and 5 (c).
In 2D, the correlation function is defined by
〈Fcorr(δr)〉avg =
drF+(r)F−(r + δr)
drρ(r)ρ(r + δr)
, (56)
which are shown in Figs. 6 (a) and 6 (b). We find that as
in the 1D case the most unstable wave length is reflected
in the shape of the spin correlation function (56), and the
characteristics of these correlation functions in the ra-
1 10 100
R [µm]
q=qc/2
(a) q = 0 (b) q = q  / 2c
-50 500
x [µm]
FIG. 6: (Color) (a) Spin correlation function defined in
Eq. (56) at t = 100 ms for q = 0 and (b) for q = qc/2.
(c) The variance of the winding number along the circumfer-
ence of the circle of radius R. The dashed lines and dotted
lines are proportional to R and R2, respectively. In (a)-(c)
the parameters are the same as those in Fig. 5, and the data
are averages over 1000 runs of simulations for different initial
states produced by random numbers. The error bars in (c)
represent the 95% confidence interval of the χ2 distribution.
dial direction are similar to those in 1D shown in Fig. 2.
For q = 0, the mean distance between spin vortices in
Fig. 5 (b) is not determined by the correlation length
(the whole width of the concentric pattern in Fig. 6 (a))
but by ∼ k−1mu, i.e., the width of the concentric rings in
Fig. 6 (a). On the other hand, for q = qc/2, the density
of spin vortices is determined by the correlation length,
i.e., the size of the blue circle ≃ 30 µm in Fig. 6 (b). The
staggered concentric correlation for q = 0 suggests that
the spin is conserved locally within the region of the cor-
relation length, and domains at a distance larger than the
correlation length grow independently, while conserving
the total spin.
The spin winding number for 2D is defined as
w(R) =
2i|F+|2
(F−∇F+ − F+∇F−) · dr,
where C(R) is a circle with radius R < Rw located at the
center of the system. Figure 6 (c) shows the R depen-
dence of the ensemble average of w2(R), where the radius
of the system is fixed to Rw = 100 µm and the data are
taken at t = 100 ms. It should be noted that 〈w2(R)〉avg
is proportional to R for large R, as expected from the
KZ theory [2], while it is proportional to R2 for small R.
This R2 dependence is due to the fact that the probabil-
ity P for a spin vortex to be in the circle is proportional
to πR2. The variance of the winding number is therefore
0(1−P )+12P/2+(−1)2P/2 ∝ R2, if the probability that
two or more vortices enter the circle is negligible. This
condition is met when the density of spin vortices times
πR2 is much smaller than unity, and hence the radius
R at which the crossover from 〈w2(R)〉avg ∝ R to ∝ R2
occurs is larger for q = qc/2 than for q = 0. As in 1D,
nonzero kmu enhances the winding of magnetization, and
the winding number is larger for q = 0 than for q = qc/2.
Figures 5 (b) and 5 (c) obviously show that the density
of spin vortices is uniform when the size of the system is
large enough. The number of spin vortices in a radius R
is therefore proportional to R2. If the topological charge
of each spin vortex, +1 or −1, was chosen at random,
the net winding number along the circle of radius R, i.e.,
the difference between the numbers of +1 and −1 vor-
tices would be proportional to R. However, from Fig. 6
(c), the winding number is proportional to R1/2 for large
R, consistent with the KZ mechanism. The topological
charge of each spin vortex is thus not at random but
anticorrelated to each other to reduce the net winding
number.
Figure 7 shows the result of the slow quench for 2D,
where q(t) is given by Eq. (30). The winding number fol-
lows the scaling law, 〈w2〉avg ∝ τ−1/3Q , as predicted from
Eq. (38), indicating that Zurek’s discussion is applicable
also to 2D. In order to obtain this scaling law, we must
specify the time at which the winding number is taken,
since the spin winding number decays in time, as shown
in the inset of Fig. 7. From the scaling law in Eq. (37),
we specify the time to take the winding number as
= const., (58)
which is indicated by the arrows in the inset of Fig. 7.
IV. CONCLUSIONS
In this paper, we have studied the dynamics of a spin-1
BEC with a ferromagnetic interaction after quench of the
applied magnetic field in an attempt to investigate spon-
taneous defect formation in the spinor BEC. We have
analyzed the magnetization triggered by quantum fluc-
tuations using the Bogoliubov approximation, and per-
formed numerical simulations of the GP equation with
initial conditions that simulate quantum fluctuations.
We have shown that the correlation functions of the
magnetization have finite correlation lengths (Figs. 2, 6
(a), and 6 (b)), and therefore magnetic domains far from
τQ [s]
slope = -1/3
0 100 200 300 400
t [ms]
τQ=1 2 4 8 [s]
FIG. 7: (Color online) (a) Variance of the spin winding num-
ber versus the quench time τQ for the 2D disk geometry, where
q is varied as in Eq. (30). The inset shows time evolution of
〈w2〉avg. The plots are taken at the times when t/τ
constant is satisfied, which are shown by the arrows in the
inset. The dashed line is proportional to τ
. The radius
of the disk is Rw = 400 µm and the closed path for taking
the winding number is R = 320 µm. The atomic density is
n = 2.8×1014 cm−3 and the number of atoms is N = 1.6×108.
The data are averages over 1000 runs of simulations for dif-
ferent initial states produced by random numbers. The error
bars represent the 95% confidence interval of the χ2 distribu-
tion.
each other grow in random directions. We find that topo-
logical defects — spin vortices — emerge through the KZ
mechanism. We have confirmed that the winding num-
ber along the closed path is proportional to the square
root of the length of the path (Figs. 2 (b) and 6 (c)),
indicating that the topological defects are formed from
domains with random directions of magnetizations.
Even when the total magnetization is conserved for
q = 0, the winding number has the same dependence on
the length of the path (Fig. 2 (b)). This is due to the
fact that domains within the correlation length tend to
be aligned in such a manner as to cancel out local mag-
netization, and consequently the total magnetization is
conserved. Thus, the neighboring domains have local cor-
relation, while domains far from each other are indepen-
dent, which makes the KZ mechanism compatible with
the total spin conservation. The formation of the local
correlation also creates topological defects as well as the
KZ mechanism, and the winding number exhibits the q
dependence as shown in Fig. 3.
When the magnetic field is quenched in finite time τQ
as in Eq. (30), the winding number has been shown to
be proportional to τ
Q (Figs. 4 and 7). This τQ depen-
dence of the winding number can be understood from
Zurek’s simple discussion [2]: the domains are frozen at
which the spin relaxation time becomes the same order
of elapsed time.
In the Berkeley experiment [12], the system is an elon-
gated quasi-2D geometry, and not suitable for testing
the KZ mechanism. The KZ mechanism should apply
to the system in which the size of the system in the x
direction is made much larger. In this case, the har-
monic potential may affect the scaling law, which mer-
its further study. Moreover in the experiment, from the
analysis in Ref. [13], there are some initial atoms in the
m = ±1 components with long-range correlation, which
play a role of seeds for large domains and hinder the ob-
servation of the KZ mechanism. If the residual atoms in
the m = ±1 components is eliminated completely, mag-
netization is triggered by quantum fluctuations as shown
in the present paper. Another way to remove the effect
of the residual atoms may be applying random phases to
the m = ±1 states to erase the initial correlation.
Note added. After our work was completed, the
preprint by Damski and Zurek [20] appeared, which per-
forms 1D simulations of the quench dynamics of a spin-1
Acknowledgments
This work was supported by Grants-in-Aid for Scien-
tific Research (Grant Nos. 17740263 and 17071005) and
by the 21st Century COE programs on “Coherent Op-
tical Science” and “Nanometer-Scale Quantum Physics”
from the Ministry of Education, Culture, Sports, Science
and Technology of Japan. MU acknowledges support by
a CREST program of the JST.
[1] T. W. B. Kibble, J. Phys. A 9, 1387 (1976).
[2] W. H. Zurek, Nature (London) 317, 505 (1985); Phys.
Rep. 276, 177 (1996).
[3] I. Chuang, R. Durrer, N. Turok, and B. Yurke, Science
251, 1336 (1991).
[4] M. J. Bowick, L. Chandar, E. A. Schiff, and A. M. Sri-
vastava, Science 263, 943 (1994).
[5] P. C. Hendry, N. S. Lawson, R. A. M. Lee, P. V. E. Mc-
Clintock, and C. D. H. Williams, Nature (London) 368,
315 (1994); M. E. Dodd, P. C. Hendry, N. S. Lawson,
P. V. E. McClintock, and C. D. H. Williams, Phys. Rev.
Lett. 81, 3703 (1998).
[6] V. M. H. Ruutu, V. B. Eltsov, A. J. Gill, T. W. B. Kibble,
M. Krusius, Yu. G. Makhlin, B. Plaçais, G. E. Volovik,
and W. Xu, Nature (London) 382, 334 (1996); V. M.
Ruutu, V. B. Eltsov, M. Krusius, Yu. G. Makhlin, B.
Plaçais, and G. E. Volovik, Phys. Rev. Lett. 80, 1465
(1998).
[7] C. Bäuerle, Yu. M. Bunkov, S. N. Fisher, H. Godfrin,
and G. R. Pickett, Nature (London) 382, 332 (1996).
[8] S. Ducci, P. L. Ramazza, W. González-Viñas, and F. T.
Arecchi, Phys. Rev. Lett. 83, 5210 (1999).
[9] R. Carmi, E. Polturak, and G. Koren, Phys. Rev. Lett.
84, 4966 (2000).
[10] R. Monaco, J. Mygind, and R. J. Rivers, Phys. Rev. Lett.
89, 080603 (2002).
[11] A. Maniv, E. Polturak, and G. Koren, Phys. Rev. Lett.
91, 197001 (2003).
[12] L. E. Sadler, J. M. Higbie, S. R. Leslie, M. Vengalattore,
and D. M. Stamper-Kurn, Nature (London) 443, 312
(2006).
[13] H. Saito, Y. Kawaguchi and M. Ueda, Phys. Rev. A 75,
013621 (2007).
[14] A. Lamacraft, cond-mat/0611017.
[15] M. Uhlmann, R. Schützhold, and U. R. Fischer,
cond-mat/0612664.
[16] H. Saito, Y. Kawaguchi and M. Ueda, Phys. Rev. Lett.
96, 065302 (2006).
[17] J. Stenger, S. Inouye, D. M. Stamper-Kurn, H. -J. Mies-
ner, A. P. Chikkatur, and W. Ketterle, Nature (London)
396, 345 (1998).
[18] K. Murata, H. Saito, and M. Ueda, Phys. Rev. A 75,
013607 (2007).
[19] T. Kuga, Y. Torii, N. Shiokawa, T. Hirano, Y. Shimizu,
and H. Sasada, Phys. Rev. Lett. 78, 4713 (1997).
[20] B. Damski and W. H. Zurek, arXiv:0704.0440.
http://arxiv.org/abs/cond-mat/0611017
http://arxiv.org/abs/cond-mat/0612664
http://arxiv.org/abs/0704.0440

  The spin vortices are shown to be created through the Kibble-Zurek (KZ)
mechanism in a quantum phase transition of a spin-1 ferromagnetic Bose-Einstein
condensate, when the applied magnetic field is quenched below a critical value.
It is shown that the magnetic correlation functions have finite correlation
lengths, and magnetizations at widely separated positions grow in random
directions, resulting in spin vortices. We numerically confirm the scaling law
that the winding number of spin vortices is proportional to the square root of
the length of the closed path, and for slow quench, proportional to $\tau_{\rm
Q}^{-1/6}$ with $\tau_{\rm Q}$ being the quench time. The relation between the
spin conservation and the KZ mechanism is discussed.

Kibble-Zurek mechanism in a quenched ferromagnetic Bose-Einstein condensate
Hiroki Saito1, Yuki Kawaguchi2, and Masahito Ueda2,3
Department of Applied Physics and Chemistry, The University of Electro-Communications, Tokyo 182-8585, Japan
Department of Physics, Tokyo Institute of Technology, Tokyo 152-8551, Japan
Macroscopic Quantum Control Project, ERATO, JST, Bunkyo-ku, Tokyo 113-8656, Japan
(Dated: October 31, 2018)
The spin vortices are shown to be created through the Kibble-Zurek (KZ) mechanism in a quantum
phase transition of a spin-1 ferromagnetic Bose-Einstein condensate, when the applied magnetic field
is quenched below a critical value. It is shown that the magnetic correlation functions have finite
correlation lengths, and magnetizations at widely separated positions grow in random directions,
resulting in spin vortices. We numerically confirm the scaling law that the winding number of spin
vortices is proportional to the square root of the length of the closed path, and for slow quench,
proportional to τ
with τQ being the quench time. The relation between the spin conservation
and the KZ mechanism is discussed.
PACS numbers: 03.75.Mn, 03.75.Lm, 73.43.Nq, 64.60.Ht
I. INTRODUCTION
Spontaneous symmetry breaking in a phase transition
produces local domains of an order parameter. If do-
mains are separated by such a long distance that they
cannot exchange information, local domains grow ini-
tially with random phases and eventually give rise to
topological defects when they overlap. This scenario
of topological-defect formation in continuous-symmetry
breaking is known as the Kibble-Zurek (KZ) mecha-
nism [1, 2], which originally predicted the cosmic-string
and monopole formation in the early Universe [1], and
has since been applied to a wide variety of systems. Ex-
perimentally, the KZ mechanism has been examined in
liquid crystals [3, 4], superfluid 4He [5] and 3He [6, 7],
an optical Kerr medium [8], Josephson junctions [9, 10],
and superconducting films [11].
Recently, spontaneous magnetization in a spinor Bose-
Einstein condensate (BEC) has attracted much interest
as a new system for studying the KZ mechanism [12, 13,
14, 15]. In the experiment performed by the Berkeley
group [12], a BEC of F = 1 87Rb atoms are prepared in
the m = 0 state, where F is the hyperfine spin and m
is its projection on the direction of the magnetic field.
By quench of the magnetic field, say in the z direction,
magnetization appears in the x-y plane. Since the spinor
Hamiltonian is axisymmetric with respect to the z axis,
the magnetization in the x-y direction breaks the U(1)
symmetry in the spin space. Thus, local domain forma-
tion is expected to lead to topological defects — spin
vortices — through the KZ mechanism.
However, the origin of the spin vortices observed after
the quench in the Berkeley experiment [12] cannot be at-
tributed to the KZ mechanism. In fact, in Ref. [12], the
spin correlation extends over the entire system (at least
in the x direction) and the domains are not independent
with each other. We have shown in Ref. [13] that the ori-
gin of the observed spin vortices is initial spin correlation
due to the residual m = ±1 atoms, which forms domain
structure followed by spin vortex creation [16]. In order
to realize the KZ mechanism in this system, i.e., in order
to ensure that the magnetic domains grow independently,
the size of the system must be much larger than the spin
correlation length and the long-range correlation in the
initial spin state must be absent. The aim of the present
paper is to show that under these conditions spin vortices
are generated through the KZ mechanism.
In the present paper we will consider 1D-ring and 2D-
disk geometries. We will show that in the 1D ring the
average spin winding number after the quench is propor-
tional to the square root of the system size, which is in
agreement with the KZ prediction [2]. In 2D the winding
number along a path with radius R is also proportional
to R1/2 as long as R is much larger than the vortex spac-
ing, while it is proportional to R for small R. When
the magnetic field is quenched slowly, the winding num-
ber is shown to be proportional to τ
Q with τQ being
the quench time. This power law can be understood by
Zurek’s simple discussion [2].
The spinor BEC is different from the other systems in
which the KZ mechanism has been observed, in that the
total spin is conserved when the quadratic Zeeman en-
ergy q is negligible. This fact is seemingly incompatible
with the KZ postulate, since the magnetic domains must
be correlated with each other so that the total magnetiza-
tion vanishes. We will show that for q = 0 small magnetic
domains are aligned to cancel out the local spin averaged
over the correlation length, and that they are indepen-
dent with each other over a greater length scale; the spin
conservation is thus compatible with the KZ mechanism.
The present paper is organized as follows. Section II
analyzes spontaneous magnetization of a spin-1 BEC and
the resultant magnetic correlation functions using the
Bogoliubov approximation. Section III performs numeri-
cal simulations of the dynamics of quenched BECs in 1D
and 2D, and shows that the KZ mechanism does emerge
in the present system. Section IV provides conclusions.
http://arxiv.org/abs/0704.1377v1
II. BOGOLIUBOV ANALYSIS OF A
QUENCHED FERROMAGNETIC
BOSE-EINSTEIN CONDENSATE
A. Hamiltonian for the spin-1 atoms
We consider spin-1 bosonic atoms with mass M con-
fined in a potential Vtrap(r). The noninteracting part of
the Hamiltonian is given by
Ĥ0 =
ψ̂†m(r)
2 + Vtrap(r)
ψ̂m(r),
where ψ̂m(r) annihilates an atom in magnetic sublevel m
of spin at a position r.
The interaction between atoms with s-wave scattering
is described by
Ĥint =
2(r) + c1F̂
:, (2)
where the symbol :: denotes the normal order and
ρ̂(r) =
ψ̂†m(r)ψ̂m(r), (3)
F̂ (r) =
ψ̂†m(r)fmm′ ψ̂m′(r), (4)
with f = (fx, fy, fz) being the spin-1 matrices. The in-
teraction coefficients in Eq. (2) are given by
4πh̄2
a0 + 2a2
, (5a)
4πh̄2
a2 − a0
, (5b)
where aS is the s-wave scattering lengths for two colliding
atoms with total spin S.
When magnetic field B is applied, the linear Zeeman
effect rotates the spin around the direction of B at the
Larmor frequency. Since Ĥ0 and Ĥint are spin-rotation
invariant and we assume the uniform magnetic field, the
linear Zeeman term has only a trivial effect on spin dy-
namics — uniform rotation of spins about B — which is
therefore ignored. The quadratic Zeeman effects for an
F = 1 87Rb atom is described by
Ĥq =
ψ̂†m(r)
(B · f)2
ψ̂m′(r), (6)
where µB is the Bohr magneton and Ehf > 0 is the hy-
perfine splitting energy between F = 1 and F = 2. The
total Hamiltonian is given by the sum of Eqs. (1), (2),
and (6),
Ĥ = Ĥ0 + Ĥq + Ĥint. (7)
B. Time evolution in the Bogoliubov
approximation
In the initial state, all atoms are prepared in them = 0
state. We study the spin dynamics of the system using
the Bogoliubov approximation with respect to this initial
state. For simplicity, we assume Vtrap = 0 in this section.
In the Bogoliubov approximation, the BEC part in the
field operator is replaced by a c-number. In the present
case, we write the m = 0 component of the field operator
ψ̂0(r) = e
−ic0nt/h̄
n+ δψ̂0(r)
, (8)
where n is the atomic density. We expand ψ̂±1(r) as
ψ̂±1(r) = e
−ic0nt/h̄
eik·r â±1,k, (9)
where V is the volume of the system and â±1,k is the
annihilation operator of an atom in the m = ±1 state
with wave vector k. Keeping only up to the second order
of δψ̂0(r) and ψ̂±1(r) in the Hamiltonian, we obtain the
Heisenberg equation of motion for â±1,k as
dâ±1,k(t)
= (εk + q + c1n)â±1,k(t) + c1nâ
∓1,−k(t),
where εk = h̄
2k2/(2M) and q = µ2BB
2/(4Ehf). The mag-
netic field is assumed to be applied in the z direction. The
solution of Eq. (10) is obtained as
â±1,k(t) =
− i εk + q + c1n
â±1,k(0)
∓1,−k(0), (11)
where
(εk + q)(εk + q + 2c1n). (12)
When Ek is imaginary, the corresponding modes are
dynamically unstable and grow exponentially. Since c1 <
0 and q > 0 for F = 1 87Rb atoms, the exponential
growth occurs for
q < 2|c1|n ≡ qc. (13)
This critical value of q agrees with the phase boundary
between the polar phase and the broken-axisymmetry
phase [17, 18]. When q ≤ qc/2, the wave number of
the most unstable mode is
kmu = ±
, (14)
and when qc/2 < q < qc, kmu = 0.
C. Fast quench
We consider the situation in which q is much larger
than the other characteristic energies for t < 0, and q is
suddenly quenched below qc at t = 0. During t < 0, the
time evolution in Eq. (11) is â±1,k(t) ≃ e−iqt/h̄â±1,k(0),
and the m = ±1 state remains in the vacuum state. For
t > 0, we obtain the time evolution of the density of the
m = ±1 component as
±1(r, t)ψ̂±1(r, t)
16|Ek|2
e2|Ek|t/h̄,(15)
where the expectation value is taken with respect to the
vacuum state of the m = ±1 component. In the sec-
ond line of Eq. (15), we have kept the unstable modes
alone with k < kc ≡
2M(qc − q)/h̄ by assuming that
|Ek|t/h̄ ≫ 1. This result indicates that the m = ±1
components grow exponentially after the quench.
Since the operator ψ̂0 in Eq. (4) is replaced by
the Bogoliubov approximation, the magnetization oper-
ator F̂+ = F̂
− = F̂x + iF̂y has the form,
F̂+(r) =
1(r) + ψ̂−1(r)
. (16)
Using Eq. (11), the time evolution of the correlation func-
tion is calculated to be
F̂+(r, t)F̂−(r
′, t)
εk + q
eik·(r−r
(17a)
qc − q − εk
e2|Ek|t/h̄+ik·(r−r
′), (17b)
where in the second line we have kept the unstable modes
alone.
From the exponential factor in Eq. (17b), we see that
the sum is contributed mostly from k around the mode
with maximum |Ek|. The denominator in the summand
of Eq. (17b) is much smoother than the exponential fac-
tor if q is not close to qc, and then we approximate εk
with εmu = h̄
2k2mu/(2M) in the denominator. We expand
2|Ek|t/h̄ around kmu in the exponent as
2|Ek|t
ξ2∆k2 − 1
Ξ4∆k4
+O(∆k6),
where ∆k = k − kmu. It is clear that τ sets the time
scale for the exponential growth. The magnetization is
observed when it sufficiently grows, i.e., t ∼ τ . Replacing
the summation with the Gaussian integral in Eq. (17b),
we find that ξ represents the correlation length. For q <
qc/2, kmu is given by Eq. (14), and
, (19)
qc − 2q
. (20)
For qc/2 < q < qc, kmu = 0 and
q(qc − q)
, (21)
2q − qc
q(qc − q)
. (22)
At q = qc/2, Eqs. (20) and (22) vanish, and the ∆k
4 term
in Eq. (18) becomes important, with
Ξ = 4
2M2q2c
. (23)
We first consider a 1D system with the periodic bound-
ary condition, i.e., the 1D ring geometry. We assume that
the radius of the ring R is much larger than the domain
size, and the curvature of the ring does not affect the
dynamics.
For q < qc/2, the magnetic correlation function is cal-
culated to be
F̂+(θ, t)F̂−(θ
′, t)
cos[kmuR(θ − θ′)]
×et/τ−τR
2(θ−θ′)2/(tξ2), (24)
where τ and ξ are given by Eqs. (19) and (20), and θ and
θ′ are azimuthal angles. For qc/2 < q < qc, we obtain
F̂+(θ, t)F̂−(θ
′, t)
qc − q
et/τ−τR
2(θ−θ′)2/(tξ2)
with Eqs. (21) and (22). At q = qc/2, the correlation
function reads
F̂+(θ, t)F̂−(θ
′, t)
qc − q
τR4(θ − θ′)4
R2(θ − θ′)2
τR4(θ − θ′)4
,(26)
where Γ is the Gamma function and
0F2(a, b, z) =
Γ(a)Γ(b)
Γ(a+ j)Γ(b + j)
is the generalized hypergeometric function. Equa-
tion (26) is shown in Fig. 2(a), where Ξ gives a char-
acteristic width of the correlation function.
Next, we consider the 2D geometry. For qc/2 < q <
qc, and then kmu = 0, the integral can be performed
analytically, giving
F̂+(r, t)F̂−(r
′, t)
2πξ2t
qc − q
et/τ−τ |r−r
′|2/(tξ2),
where τ and ξ are given in Eqs. (21) and (22). For other
q, we can perform only the angular integral as
F̂+(r, t)F̂−(r
′, t)
qc − q − εmu
kJ0(k|r − r′|)e2|Ek|t/h̄dk,
where J0 is the Bessel function. If the exponential factor
is much sharper than the Bessel function around kmu,
the correlation function (29) is approximated to be ∝
J0(kmu|r − r′|) [14, 15].
As shown above, the correlation function (17b) has a
finite correlation length, and the magnetization at posi-
tions widely separated from each other grow with inde-
pendent directions in the x-y plane. Thus, the growth
of the magnetic domains is expected to leave topological
defects through the KZ mechanism.
D. Slow quench
In the previous sections, we have assumed that the
magnetic field is suddenly quenched to the desired value
at t = 0 and q is held constant for t > 0. We assume here
that for t > 0 the magnetic field is gradually quenched
q(t) = qc
. (30)
The magnetic correlation can be estimated to be
F̂+(r, t)F̂−(r
′, t)
dk exp
2|Ek(t)|t
dt+ ik · (r − r′)
. (31)
Since we are interested in the vicinity of the critical
point where correlation starts to grow, we expand |Ek(t)|
around kmu = 0 and keep the terms up to the order of
k2. For the 1D ring, we obtain
F̂+(θ, t)F̂−(θ
′, t)
∝ ef(t)−R
2(θ−θ′)2/ξ2Q , (32)
and for the 2D geometry,
F̂+(r, t)F̂−(r
′, t)
∝ ef(t)−|r−r
′|2/ξ2
Q , (33)
where
f(t) =
tan−1
τQ − t
1− 2t
, (34)
t(τQ − t)
. (35)
For t≪ τQ, f(t) can be expanded as
f(t) =
, (36)
and from f(t) ∼ 1, the time scale for magnetization is
given by
Q . (37)
Substitution of tQ into Eq. (35) yields
Q . (38)
The same power law is obtained in Ref. [14].
It is interesting to note that the results (37) and
(38) are easily obtained also by the simple discussion by
Zurek [2]. Since q(t) depends on time, τ and ξ given in
Eqs. (21) and (22) are time dependent, and hence they
are regarded as the growth time and correlation length
at each instant of time. The local magnetization is de-
veloped after a time tQ has elapsed such that
τ(tQ) ∼ tQ. (39)
Using
τ(t) =
t(τQ − t)
, (40)
we obtain tQ in Eq. (37). Substituting this tQ into
ξ2(t) =
τQ − 2t
t(τQ − t)
yields Eq. (38).
III. NUMERICAL SIMULATIONS AND THE
KIBBLE-ZUREK MECHANISM
A. Gross-Pitaevskii equation with quantum
fluctuations
The multicomponent Gross-Pitaevskii (GP) equation
is obtained by replacing the field operators ψ̂m with the
macroscopic wave function ψm in the Heisenberg equa-
tion of motion:
∇2 + Vtrap + q + c0ρ
F∓ψ0 ± Fzψ±1
, (42a)
∇2 + Vtrap + c0ρ
(F+ψ1 + F−ψ−1) , (42b)
where ρ and F are defined using ψm instead of ψ̂m in
Eqs. (3) and (4). The wave function is normalized as
|ψm|2 = N, (43)
with N being the number of atoms in the condensate.
Suppose that all atoms are initially in the m = 0 state.
It follows then from Eq. (42a) that ψ±1 will remain zero
in the subsequent time evolution. This is because quan-
tum fluctuations in the transverse magnetization that
trigger the growth of magnetization are neglected in the
mean-field approximation. We therefore introduce an ap-
propriate initial noise in ψ±1 so that the mean-field ap-
proximation reproduces the quantum evolution.
Let us write the initial state as
ψ±1(r) =
eik·ra±1,k(0), (44)
where a±1,k are c-numbers. We assume that the c-
number amplitudes a±1,k(0) are stochastic variables
whose average values vanish,
〈a±1,k(0)〉avg = 0, (45)
where by 〈· · · 〉avg we denote the statistical average over
an appropriate probability distribution. The linear ap-
proximation of the GP equation with respect to a±1,k
gives the same time evolution as Eq. (11), in which the
operators are replaced by the c-numbers. We thus obtain
F+(r, t)F−(r
′, t) =
εk + q
e−ik·(r−r
′)|a1,k(0)|2
+eik·(r−r
′)|a−1,−k(0)|2
. (46)
Comparing Eq. (46) with Eq. (17a), we find that they
have the same form if the variance of a±1,k(0) satisfies
〈|a±1,k(0)|2〉avg =
for all k.
0 100 200 300
t [ms]
q=qc/2
0 π 2π
 |F+| /ρ arg F+
q=qc/ 20
0 π 2π
 |F+| /ρ arg F+
(b) (c)
FIG. 1: (Color online) (a) Time evolution of the auto cor-
relation function given in Eq. (49) for the 1D ring geometry.
(b) Magnitude of the normalized magnetization |F+|/ρ (solid
curve, left axis) and direction of the magnetization argF+
(dashed curve, right axis) at t = 70 ms for q = 0 and (c) for
q = qc/2. The radius of the ring is R = 50 µm, the atomic
density is n = 2.8 × 1014 cm−3, and the number of atoms is
N = 106.
In the following, we will perform numerical simulation
of spontaneous magnetization using the GP equation and
show that the ensuing dynamics exhibits defect formation
similar to the KZ mechanism. As the initial state of the
m = ±1 wave functions, we use Eq. (44) with
a±1,k(0) = αrnd + iβrnd, (48)
where αrnd and βrnd are random variables following the
normal distribution p(x) =
2/π exp(−2x2). Equation
(48) then satisfies Eqs. (45) and (47).
B. 1D ring geometry
Let us first investigate the 1D ring system. Experi-
mentally this geometry can be realized, e.g., by an opti-
cal trap using a Laguerre-Gaussian beam [19]. We reduce
the GP equation (42) to 1D by assuming that the wave
function ψm depends only on the azimuthal angle θ. The
average density of atoms is assumed to be n = 2.8× 1014
cm−3. When the radius of the ring is R = 50 µm and
the radius of the small circle is 2 µm, the total number
of atoms is N ≃ 106.
Figure 1 illustrates a single run of time evolution for
an initial state given by Eqs. (44) and (48). Figure 1
(a) shows time evolution of the auto correlation function
defined by
F̄ (t) =
|F+(θ, t)|2
ρ2(θ, t)
. (49)
For both q = 0 and q = qc/2, the transverse magneti-
zation grows exponentially with a time constant ∼ τ =
h̄/qc ≃ 8 ms. Snapshots of the transverse magnetization
at t = 70 ms are shown in Figs. 1 (b) and 1 (c) for q = 0
and q = qc/2, respectively. We define the spin winding
number as
2i|F+|2
, (50)
which represents the number of rotation of the spin vector
in the x-y plane along the circumference, and of course
w is an integer. The spin winding numbers are w = 7 in
Fig. 1 (b) and w = −1 in Fig. 1 (c).
Figure 2 (a) shows the ensemble average of the nor-
malized correlation function,
〈Fcorr(δθ)〉avg =
dθF+(θ)F−(θ + δθ)
dθρ(θ)ρ(θ + δθ)
, (51)
at t = 70 ms. For q = qc/2, the correlation function has
the characteristic width of ∼ Ξ in Eq. (23), indicating
that the ring is filled with magnetic domains with an av-
erage size of ∼ Ξ. According to the KZ theory, the mag-
netic domains with random directions give rise to the spin
winding, which is estimated to be w ∼ (R/Ξ)1/2. This R
dependence of w is clearly seen in Fig. 2 (b). The ensem-
ble average of the winding number, 〈w〉avg , vanishes due
to the random nature of the initial noise, and the square
root of its variance, 〈w2〉1/2avg, should be regarded as a typ-
ical winding number. The variance is expected to obey
the χ2 distribution with 1000 degrees of freedom, and
hence we show the 95% confidence interval to estimate
the statistical errors in Fig. 2. As shown in the inset of
Fig. 2 (b), the typical winding number changes in time,
since the ferromagnetic energy is converted to the kinetic
energy and the system exhibits complicated dynamics.
The situation is different for q = 0, in which the cor-
relation function oscillates with a Gaussian envelope as
shown in Fig. 2 (a). This form of the correlation function
gives us the answer to the question as to how the KZ
mechanism manifests itself in spin conserving systems.
The finite correlation length for q = 0 indicates that the
spin is conserved not only globally but also locally, that
is, the locally integrated spin over the correlation length
|δr|<∼ξ
F (r + δr)dδr, (52)
is held to be zero for any r. This local spin conservation
is due to formation of staggered domain or helical spin
structures whose periodic length is much smaller than
−π/2 0 π/2
q=qc/2
10 100 1000
R [µm]
q=qc/2
0 100
t [ms]
q=qc/2
FIG. 2: (Color online) (a) Numerically obtained correlation
function given in Eq. (51) at t = 70 ms (solid curves), and
theoretical fits (dashed curves) from Eqs. (24) and (26). Other
parameters are the same as those in Fig. 1. (b) R dependence
of the variance of the spin winding number, where the number
of atoms is related to R as N = 106×R [µm] /50. The dashed
lines are semi-log fits to the numerical data. The inset shows
the time dependence of 〈w2〉avg for R = 50 µm. The data
in (a) and (b) are averages over 1000 runs of simulations for
different initial states produced by random numbers. The
error bars in (b) represent the 95% confidence interval of the
χ2 distribution.
ξ. Thus, the neighboring domains tend to have opposite
magnetizations to cancel out the spin locally, and the
domains far from each other grow independently; the spin
conservation and the KZ mechanism are thus compatible.
The oscillation in the correlation function originates
from the fact that the most unstable modes have nonzero
wave numbers ±kmu. Each correlated region of size ∼
ξ = [8h̄2/(Mqc)]
1/2 contains spin waves of eikmuRθ and
e−ikmuRθ. If there is an imbalance between these modes,
the winding number monotonically increases or decreases
in each region of ∼ ξ. This is the reason why 〈w2〉avg is
larger for q = 0 than for q = qc/2 in Fig. 2 (b). It follows
from this consideration that for kmuξ ≫ 1 the winding
0.1 1
1−2q/qc
slope = 3/2
FIG. 3: (Color online) Dependence of the variance of the
spin winding number on q. Except for q, the parameters are
the same as those in Fig. 1. The dashed line is proportional
to (1−2q/qc)
3/2. The plots show the averages over 1000 runs
of simulations for different initial states produced by random
numbers. The error bars represent the 95% confidence inter-
val of the χ2 distribution.
number is proportional to
w ∼ kmuξ
= kmu
1− 2q
, (53)
where Eqs. (14) and (20) are used. Figure 3 shows the
averaged variance of the winding number versus 1−2q/qc.
For small q, 〈w2〉avg is proportional to (1 − 2q/qc)3/2, in
agreement with Eq. (53). When q is close to qc/2, the
spin winding within the correlated region, kmuξ, becomes
small, and then the winding number reduces to the value
shown in Fig. 2 (b), i.e., 〈w2〉avg ≃ 4.
We next discuss the results of simulations of slow
quench as in Eq. (30). Since the winding number for
the slow quench is small compared with the fast quench,
we take a large ring of R = 400 µm. Figure 4 shows
the variance of the winding number as a function of the
quench time. We can clearly see that 〈w2〉avg has a power
law of τ
Q within the statistical error, which is in agree-
ment with ξ−1Q ∼ τ
Q , with ξQ being given in Eq. (38).
Thus, the present system follow the quench-time scaling
of Zurek [2]. We note that, in the slow quench, the wind-
ing number converges to an almost constant value for
varying quench time τQ, as shown in the inset of Fig. 4.
This is because little excess energy other than for exciting
spin vortices is available for the slow quench.
C. 2D disk geometry
When the confinement in the z direction is tight, the
system is effectively 2D. For simplicity, we ignore the den-
sity dependence in the z direction, and assume that the
τQ[s]
slope=−1/3
0 100 200 300 400
t [ms]
0.8 s
1.6 s
3.2 s
6.4 s
FIG. 4: (Color online) Dependence of the variance of the spin
winding number at t = 400 ms on the quench time τQ, where q
is varied as in Eq. (30). The radius of the ring is R = 400 µm,
the atomic density is n = 2.8 × 1014 cm−3, and the number
of atoms is N = 8 × 106. The dashed line is proportional
. The inset shows time evolution of 〈w2〉avg. The
data are averages over 1000 runs of simulations for different
initial states produced by random numbers. The error bars
represent the 95% confidence interval of the χ2 distribution.
2D GP equation has the same form as Eq. (42). We as-
sume that the wave function vanishes at the wall located
at (x2 + y2)1/2 = Rw = 100 µm, and that the potential
is flat inside of the wall. Then the density n = 2.8× 1014
cm−3 is almost constant except within the healing length
{3/[8πn(a0 + 2a2)]}1/2 ≃ 0.16 µm near the wall. When
the thickness in the z direction is ≃ 1 µm, the number of
atoms is N ≃ 107. Such a system will be realized using
an optical sheet and a hollow laser beam.
The initial state of ψ0 is a stationary solution of the GP
equation, and the initial state of ψ±1 is given by Eq. (44)
with random variables (48). Figure 5 (a) shows time
evolution of the autocorrelation function of the transverse
magnetization,
F̄ (t) =
|F+(r, t)|2
ρ2(r, t)
, (54)
which grows exponentially with the same time constant
as that in Fig. 1, and saturates for t >∼ 100 ms.
Snapshots of |F+(r)| and argF+(r) at t = 100 ms are
shown in Figs. 5 (b) and 5 (c). We see that |F+(r)|
at t >∼ 100 ms contains many holes, around which the
spin direction rotates by 2π. Since this topological spin
structure consists of singly-quantized vortices in the m =
±1 states filled by atoms in the m = 0 state, it is called
the “polar-core vortex.” We can estimate the spin healing
length ξs by equating the kinetic energy h̄
2/(2Mξ2s ) with
0 100 200 300
t [ms]
q=0 q=qc/2
(c) q = q  / 2
(b) q = 0
t = 50 ms t = 100 mst = 75 ms t = 200 ms
t = 50 ms t = 100 mst = 75 ms t = 200 ms
FIG. 5: (Color) (a) Time evolution of the autocorrelation function given in Eq. (54) for the 2D disk geometry. The radius of
the disk is Rw = 100 µm, the atomic density is n = 2.8× 10
14 cm−3, and the number of atoms is N = 107. (b) Profiles of the
magnetization |F+| (upper) and its direction argF+ (lower) for q = 0 and (c) for q = qc/2. The size of each panel is 200 µm
×200 µm.
the energy of magnetization |q − qc|, giving
2M |q − qc|
. (55)
This length scale is ξs ≃ 1.7 µm for q = 0 and ξs ≃ 2.4
µm for q = qc/2, which are in good agreement with the
sizes of the vortex cores in Figs. 5 (b) and 5 (c).
In 2D, the correlation function is defined by
〈Fcorr(δr)〉avg =
drF+(r)F−(r + δr)
drρ(r)ρ(r + δr)
, (56)
which are shown in Figs. 6 (a) and 6 (b). We find that as
in the 1D case the most unstable wave length is reflected
in the shape of the spin correlation function (56), and the
characteristics of these correlation functions in the ra-
1 10 100
R [µm]
q=qc/2
(a) q = 0 (b) q = q  / 2c
-50 500
x [µm]
FIG. 6: (Color) (a) Spin correlation function defined in
Eq. (56) at t = 100 ms for q = 0 and (b) for q = qc/2.
(c) The variance of the winding number along the circumfer-
ence of the circle of radius R. The dashed lines and dotted
lines are proportional to R and R2, respectively. In (a)-(c)
the parameters are the same as those in Fig. 5, and the data
are averages over 1000 runs of simulations for different initial
states produced by random numbers. The error bars in (c)
represent the 95% confidence interval of the χ2 distribution.
dial direction are similar to those in 1D shown in Fig. 2.
For q = 0, the mean distance between spin vortices in
Fig. 5 (b) is not determined by the correlation length
(the whole width of the concentric pattern in Fig. 6 (a))
but by ∼ k−1mu, i.e., the width of the concentric rings in
Fig. 6 (a). On the other hand, for q = qc/2, the density
of spin vortices is determined by the correlation length,
i.e., the size of the blue circle ≃ 30 µm in Fig. 6 (b). The
staggered concentric correlation for q = 0 suggests that
the spin is conserved locally within the region of the cor-
relation length, and domains at a distance larger than the
correlation length grow independently, while conserving
the total spin.
The spin winding number for 2D is defined as
w(R) =
2i|F+|2
(F−∇F+ − F+∇F−) · dr,
where C(R) is a circle with radius R < Rw located at the
center of the system. Figure 6 (c) shows the R depen-
dence of the ensemble average of w2(R), where the radius
of the system is fixed to Rw = 100 µm and the data are
taken at t = 100 ms. It should be noted that 〈w2(R)〉avg
is proportional to R for large R, as expected from the
KZ theory [2], while it is proportional to R2 for small R.
This R2 dependence is due to the fact that the probabil-
ity P for a spin vortex to be in the circle is proportional
to πR2. The variance of the winding number is therefore
0(1−P )+12P/2+(−1)2P/2 ∝ R2, if the probability that
two or more vortices enter the circle is negligible. This
condition is met when the density of spin vortices times
πR2 is much smaller than unity, and hence the radius
R at which the crossover from 〈w2(R)〉avg ∝ R to ∝ R2
occurs is larger for q = qc/2 than for q = 0. As in 1D,
nonzero kmu enhances the winding of magnetization, and
the winding number is larger for q = 0 than for q = qc/2.
Figures 5 (b) and 5 (c) obviously show that the density
of spin vortices is uniform when the size of the system is
large enough. The number of spin vortices in a radius R
is therefore proportional to R2. If the topological charge
of each spin vortex, +1 or −1, was chosen at random,
the net winding number along the circle of radius R, i.e.,
the difference between the numbers of +1 and −1 vor-
tices would be proportional to R. However, from Fig. 6
(c), the winding number is proportional to R1/2 for large
R, consistent with the KZ mechanism. The topological
charge of each spin vortex is thus not at random but
anticorrelated to each other to reduce the net winding
number.
Figure 7 shows the result of the slow quench for 2D,
where q(t) is given by Eq. (30). The winding number fol-
lows the scaling law, 〈w2〉avg ∝ τ−1/3Q , as predicted from
Eq. (38), indicating that Zurek’s discussion is applicable
also to 2D. In order to obtain this scaling law, we must
specify the time at which the winding number is taken,
since the spin winding number decays in time, as shown
in the inset of Fig. 7. From the scaling law in Eq. (37),
we specify the time to take the winding number as
= const., (58)
which is indicated by the arrows in the inset of Fig. 7.
IV. CONCLUSIONS
In this paper, we have studied the dynamics of a spin-1
BEC with a ferromagnetic interaction after quench of the
applied magnetic field in an attempt to investigate spon-
taneous defect formation in the spinor BEC. We have
analyzed the magnetization triggered by quantum fluc-
tuations using the Bogoliubov approximation, and per-
formed numerical simulations of the GP equation with
initial conditions that simulate quantum fluctuations.
We have shown that the correlation functions of the
magnetization have finite correlation lengths (Figs. 2, 6
(a), and 6 (b)), and therefore magnetic domains far from
τQ [s]
slope = -1/3
0 100 200 300 400
t [ms]
τQ=1 2 4 8 [s]
FIG. 7: (Color online) (a) Variance of the spin winding num-
ber versus the quench time τQ for the 2D disk geometry, where
q is varied as in Eq. (30). The inset shows time evolution of
〈w2〉avg. The plots are taken at the times when t/τ
constant is satisfied, which are shown by the arrows in the
inset. The dashed line is proportional to τ
. The radius
of the disk is Rw = 400 µm and the closed path for taking
the winding number is R = 320 µm. The atomic density is
n = 2.8×1014 cm−3 and the number of atoms is N = 1.6×108.
The data are averages over 1000 runs of simulations for dif-
ferent initial states produced by random numbers. The error
bars represent the 95% confidence interval of the χ2 distribu-
tion.
each other grow in random directions. We find that topo-
logical defects — spin vortices — emerge through the KZ
mechanism. We have confirmed that the winding num-
ber along the closed path is proportional to the square
root of the length of the path (Figs. 2 (b) and 6 (c)),
indicating that the topological defects are formed from
domains with random directions of magnetizations.
Even when the total magnetization is conserved for
q = 0, the winding number has the same dependence on
the length of the path (Fig. 2 (b)). This is due to the
fact that domains within the correlation length tend to
be aligned in such a manner as to cancel out local mag-
netization, and consequently the total magnetization is
conserved. Thus, the neighboring domains have local cor-
relation, while domains far from each other are indepen-
dent, which makes the KZ mechanism compatible with
the total spin conservation. The formation of the local
correlation also creates topological defects as well as the
KZ mechanism, and the winding number exhibits the q
dependence as shown in Fig. 3.
When the magnetic field is quenched in finite time τQ
as in Eq. (30), the winding number has been shown to
be proportional to τ
Q (Figs. 4 and 7). This τQ depen-
dence of the winding number can be understood from
Zurek’s simple discussion [2]: the domains are frozen at
which the spin relaxation time becomes the same order
of elapsed time.
In the Berkeley experiment [12], the system is an elon-
gated quasi-2D geometry, and not suitable for testing
the KZ mechanism. The KZ mechanism should apply
to the system in which the size of the system in the x
direction is made much larger. In this case, the har-
monic potential may affect the scaling law, which mer-
its further study. Moreover in the experiment, from the
analysis in Ref. [13], there are some initial atoms in the
m = ±1 components with long-range correlation, which
play a role of seeds for large domains and hinder the ob-
servation of the KZ mechanism. If the residual atoms in
the m = ±1 components is eliminated completely, mag-
netization is triggered by quantum fluctuations as shown
in the present paper. Another way to remove the effect
of the residual atoms may be applying random phases to
the m = ±1 states to erase the initial correlation.
Note added. After our work was completed, the
preprint by Damski and Zurek [20] appeared, which per-
forms 1D simulations of the quench dynamics of a spin-1
Acknowledgments
This work was supported by Grants-in-Aid for Scien-
tific Research (Grant Nos. 17740263 and 17071005) and
by the 21st Century COE programs on “Coherent Op-
tical Science” and “Nanometer-Scale Quantum Physics”
from the Ministry of Education, Culture, Sports, Science
and Technology of Japan. MU acknowledges support by
a CREST program of the JST.
[1] T. W. B. Kibble, J. Phys. A 9, 1387 (1976).
[2] W. H. Zurek, Nature (London) 317, 505 (1985); Phys.
Rep. 276, 177 (1996).
[3] I. Chuang, R. Durrer, N. Turok, and B. Yurke, Science
251, 1336 (1991).
[4] M. J. Bowick, L. Chandar, E. A. Schiff, and A. M. Sri-
vastava, Science 263, 943 (1994).
[5] P. C. Hendry, N. S. Lawson, R. A. M. Lee, P. V. E. Mc-
Clintock, and C. D. H. Williams, Nature (London) 368,
315 (1994); M. E. Dodd, P. C. Hendry, N. S. Lawson,
P. V. E. McClintock, and C. D. H. Williams, Phys. Rev.
Lett. 81, 3703 (1998).
[6] V. M. H. Ruutu, V. B. Eltsov, A. J. Gill, T. W. B. Kibble,
M. Krusius, Yu. G. Makhlin, B. Plaçais, G. E. Volovik,
and W. Xu, Nature (London) 382, 334 (1996); V. M.
Ruutu, V. B. Eltsov, M. Krusius, Yu. G. Makhlin, B.
Plaçais, and G. E. Volovik, Phys. Rev. Lett. 80, 1465
(1998).
[7] C. Bäuerle, Yu. M. Bunkov, S. N. Fisher, H. Godfrin,
and G. R. Pickett, Nature (London) 382, 332 (1996).
[8] S. Ducci, P. L. Ramazza, W. González-Viñas, and F. T.
Arecchi, Phys. Rev. Lett. 83, 5210 (1999).
[9] R. Carmi, E. Polturak, and G. Koren, Phys. Rev. Lett.
84, 4966 (2000).
[10] R. Monaco, J. Mygind, and R. J. Rivers, Phys. Rev. Lett.
89, 080603 (2002).
[11] A. Maniv, E. Polturak, and G. Koren, Phys. Rev. Lett.
91, 197001 (2003).
[12] L. E. Sadler, J. M. Higbie, S. R. Leslie, M. Vengalattore,
and D. M. Stamper-Kurn, Nature (London) 443, 312
(2006).
[13] H. Saito, Y. Kawaguchi and M. Ueda, Phys. Rev. A 75,
013621 (2007).
[14] A. Lamacraft, cond-mat/0611017.
[15] M. Uhlmann, R. Schützhold, and U. R. Fischer,
cond-mat/0612664.
[16] H. Saito, Y. Kawaguchi and M. Ueda, Phys. Rev. Lett.
96, 065302 (2006).
[17] J. Stenger, S. Inouye, D. M. Stamper-Kurn, H. -J. Mies-
ner, A. P. Chikkatur, and W. Ketterle, Nature (London)
396, 345 (1998).
[18] K. Murata, H. Saito, and M. Ueda, Phys. Rev. A 75,
013607 (2007).
[19] T. Kuga, Y. Torii, N. Shiokawa, T. Hirano, Y. Shimizu,
and H. Sasada, Phys. Rev. Lett. 78, 4713 (1997).
[20] B. Damski and W. H. Zurek, arXiv:0704.0440.
http://arxiv.org/abs/cond-mat/0611017
http://arxiv.org/abs/cond-mat/0612664
http://arxiv.org/abs/0704.0440
