Dynamics of a quantum phase transition in a ferromagnetic Bose-Einstein
  condensate
Dynamics of a quantum phase transition in a ferromagnetic Bose-Einstein condensate
Bogdan Damski and Wojciech H. Zurek
Theory Division, Los Alamos National Laboratory, MS-B213, Los Alamos, NM 87545, USA
We discuss dynamics of a slow quantum phase transition in a spin-1 Bose-Einstein condensate.
We determine analytically the scaling properties of the system magnetization and verify them with
numerical simulations in a one dimensional model.
Studies of phase transitions have traditionally focused
on equilibrium scalings of various properties near the crit-
ical point. Dynamics of the phase transition presents new
challenges and there is a strong motivation for analyzing
it. Nonequilibrium phase transitions may play a role in
the evolution of the early Universe [1]. Their analogues
can be studied in the condensed matter experiments. The
latter observation led to development of the theory based
on the universality of critical behavior [2], which in turn
resulted in a series of beautiful experiments [3]. The re-
cent progress in the cold atom experiments allows for
time dependent realizations of different models under-
going a quantum phase transition (QPT) [4, 5]. These
experimental developments are only a proverbial tip of
the iceberg, but they call for an in-depth theoretical un-
derstanding of the QPT dynamics.
A QPT is a fundamental change in ground state (GS)
of the system as a result of small variations of an exter-
nal parameter, e.g., a magnetic field [6]. It takes place
ideally at zero absolute temperature, which is in strik-
ing contrast to thermodynamical phase transitions. The
most complete description of the QPT dynamics has been
obtained so far in spin models [7, 8] that are exactly solv-
able. In these systems the gap in the excitation spectrum
goes to zero at the critical point, which precludes the adi-
abatic evolution across the phase boundary. It leads to
creation of excitations whose density and scaling with a
quench rate follow from a quantum version [7, 9] of the
Kibble-Zurek (KZ) theory [1, 2].
We study dynamics of a ferromagnetic condensate of
spin-1 particles [10]. For simplicity, we consider 1D ho-
mogeneous (untrapped) clouds: atoms in a box as in the
experiment [11] with spinless bosons. We adopt the pa-
rameters for our 1D model such that the length and time
scales are comparable to experimental ones [12]. Assum-
ing that the system is placed in a magnetic field B aligned
in the z direction, one gets the following dimensionless
mean-field energy functional [12]
E[Ψ] =
+Q〈Ψ|F 2z |Ψ〉
〈Ψ|Fα|Ψ〉2 (1)
where ΨT = (ψ1, ψ0, ψ−1) describes the m = 0,±1 con-
densate components,
m |ψm|2 = 1, and Fx,y,z are
spin-1 matrices [13]. The first term in (1) is the kinetic
energy, the second and the fourth term describe spin-
independent and spin-dependent atom interactions re-
spectively, the third term is a quadratic Zeeman shift
coming from atom interactions with a magnetic field.
For 87Rb atoms considered here c1 < 0, which results
in an interesting phase diagram due to the competition
between the last two terms in (1). Restricting analysis
to zero longitudinal magnetization case, and introducing
q = Q/(n|c1|), n = Ψ†Ψ
one finds a polar phase for q > 2, described by ΨTP ∼
(0, 1, 0), and the broken-symmetry phase where
ΨTB ∼ (
4− 2qeiχ1 ,
8 + 4qei(χ1+χ−1)/2,
4− 2qeiχ−1)
for 0 ≤ q < 2. The freedom of choosing the χ±1 results
in rotational symmetry of the transverse magnetization
on the (x, y) plane. The transition between these phases
can be driven by the change of the magnetic field B im-
posed on the atom cloud, q ∼ Q ∼ B2 [14], which was
experimentally done in [5].
The dynamics of a QPT depends on the rate of quench
driving the system across the phase boundary. For very
fast “impulse” transition, the system has no time to ad-
just to the changes of the Hamiltonian and arrives in a re-
gion where a new phase is expected with the “old” wave-
function untouched during the evolution. Slow transi-
tions are different: the system has time to “probe” vari-
ous broken symmetry “vacua” in the neighborhood of the
critical point where it gets excited. We are interested in
evolutions that are fast enough to produce macroscopic
excitations of the system, but slow enough to reflect scal-
ings of the critical region. By comparing analytical find-
ings to numerical simulations for experimentally relevant
parameters we provide the first complete description of
QPT dynamics in a ferromagnetic condensate.
Fast transitions were realized in the Berkeley experi-
ment [5]. The 3D numerical simulations closely following
this experiment were reported in [14]. Analytical studies
of the evolution after “impulse” quench were presented in
[15, 16]. The paper of Lamacraft [15] also discusses dy-
namics of non instantaneous transitions in 2D spinor con-
densates focusing on analytical predictions on the growth
of the transverse magnetization correlation functions.
We start with a qualitative discussion. Consider-
ing small perturbations around the GS of the broken-
symmetry phase one finds three Bogolubov modes as in
http://arxiv.org/abs/0704.0440v2
[13] where quantum fluctuations are studied. In the long
wavelength limit (important for slow transitions) there
is only one nonzero eigenvalue mode: the gapped mode
having eigenenergy ∆ ∼
4− q2. Suppose now that we
drive the system from polar to broken-symmetry phase.
The system reaction time to Hamiltonian changes in the
broken-symmetry phase is given by the inverse of the gap:
τ ∼ 1
[7, 9]. For example, when τ is small enough the
evolution becomes adiabatic so the system adjusts fast
to parameter changes. Right after entering the broken-
symmetry phase, the reaction time is large with respect
to the transition time, ∆/ d∆
, and so the system under-
goes the “impulse” evolution where its state is “frozen”.
The gapped mode starts to be populated around the in-
stant t̂ after entering the broken-symmetry phase: the
system leaves the “impulse” regime to catch up with in-
stantaneous GS solution. This happens when the two
time scales become comparable: 1/∆(t̂) ∼ ∆/ d∆
|t=t̂. We
consider here transitions driven by
q(t) = 2− t/τQ, (2)
where τQ is the quench time inversely proportional to the
speed of driving the system through the phase transition.
For slow transitions of interest here, τQ ≫ 1, we obtain
t̂ ∼ τ1/3Q . (3)
In the following we analyze dynamics induced by a
linear decrease of q(t) (2). The evolution starts from
t < 0, i.e., in the polar phase, and ends at t = 2τQ (q =
0). Such q(t) dependence is achieved by ramping down
the magnetic field as ∼
2− t/τQ. The initial state
is chosen as a slightly perturbed GS in the polar phase,
ΨT ∼ (δψ1, 1/
L+δψ0, δψ−1), where |δψm| ≪ 1/
L are
random. We generate the real and imaginary part of δΨm
at different grid points with the probability distribution
p(x) = exp(−x2/2σ2)/
2πσ. We take σ = 10−4 to start
evolution closely to the polar phase GS.
To find the full numerical solution within the mean-
field approximation, we integrate three coupled nonlin-
ear Schrödinger equations for the ψm condensates that
can be easily obtained by the variation of (1). During
evolution we look at the magnetization of the sample
fα = 〈Ψ|Fα|Ψ〉, α = x, y, z.
The transverse magnetization. A total transverse (to the
magnetic field in the z direction) magnetization reads
MT (t) =
dz[f2x(z, t) + f
y (z, t)] =
dz mT , (4)
and is experimentally measurable. It disappears in the
GS of the polar phase and equals (1 − q2/4)/L in the
broken-symmetry GS. Its typical evolution is depicted in
Fig. 1. We see there that nothing happens in the po-
lar phase. The system starts nontrivial evolution in the
0 0.5 1 1.5 2 2.5 3
q(t)=2-t/τ
1.5 1.75 2
0.4 1 4 16 64
FIG. 1: (color) Main plot: numerical solution (black solid
line) vs. static prediction (red dashed line). The arrow depicts
direction of evolution. Inset (a): the same as in the main plot
plus a numerically obtained solution of the linearized problem
(green divergent line). Inset (b): numerical data vs. fit to
τQ ≥ 10 data only (see text for details). In the main plot
τQ = 10 (see [12] for units).
FIG. 2: The vectors represent (fx(z), fy(z)) × 103. Plot (a):
snapshot at q(t = 2.81) = 1.72, i.e., at the first peak in MTL
(see Fig. 1). Plot (b): snapshot by the end of time evolution:
q(t = 20) = 0. The results come from the same numerical
simulation as in Fig. 1 (see [12] for units).
broken-symmetry phase at a distance t̂/τQ after the criti-
cal point was passed. The magnetization grows fast from
that point until it exceeds the static prediction and starts
oscillations with the amplitude decreasing in time. We
consider slow transitions. Therefore, by the end of time
evolution, when q = 0, the system is in the slightly per-
turbed ferromagnetic GS: globallyMTL ≈ 1 (Fig. 1) and
locally L2mT (z) ≈ 1 (Fig. 3). We can now ask: Does the
scaling (3) hold? To find out we define arbitrarily t̂ as the
instant when MTL intersects 1%. A fit to numerics for
τQ ≥ 10 yields ln t̂ = (0.056±0.01)+(0.332±0.002) lnτQ
which confirms prediction (3). This fit is presented in Fig.
1a, where the gradual departure of the numerical data for
τQ < 10 from t̂ ∼ τ1/3Q indicates that τQ ≫ 1 or 37ms has
to be taken for the observation of 1/3 exponent: quench
has to be slow enough to reflect the critical dynamics.
60 65 70 75
FIG. 3: (color) Magnetization of the system at t = 2τQ
(q = 0). The dashed lines facilitate observation of extrema
coincidences. Results come from the same simulation as in
Figs. 1, 2; see [12] for units.
In the GS configuration of the broken-symmetry phase
the vector (fx, fy) can have arbitrary orientation, so in
the dynamical problem considered here it is interesting to
find out how is this rotational symmetry broken. When
unstable evolution starts, spatial correlations in magne-
tization appear (Fig. 2a). In the subsequent evolution
these correlations evolve such that the correlation length
increases: see Fig. 2b obtained by the end of time evolu-
tion. This is a generic picture though the details depend
on the quench time τQ and initial state of the system.
This behavior suggests creation of spin textures [17, 18].
In our case, topological textures are spin configurations
where the magnetization direction varies in space so that
the kinetic energy term in (1) is not minimized, but mag-
netization magnitude follows closely a GS result. Such
structures appear in 1D when the first homotopy group
of the vacuum manifold M is nontrivial, which happens
here: π1(M) = Z [19]. These textures are characterized
by the winding number, 1
Arg(fx + ify), which
is not conserved. Indeed, it reads +1 in Fig. 2a, while
by the end of that evolution (Fig. 2b) it equals 0.
Are different stages of this evolution experimentally
observable? Let’s look at τQ = 10 case presented in Figs.
1-3. The evolution from the phase boundary to the first
peak in magnetizationMT (the q = 0 point) takes 2.81×
37ms ∼= 104ms (2τQ = 740ms). Both these time scales
are well within the reach of the experiment [5].
The longitudinal magnetization. Initially, fz(z) ≈ 0 so
dzfz ≈ 0. The conservation of the latter allows
only for creation of a network of magnetic domains (non-
topological structures with fixed fz sign) having opposite
polarizations. The domains appear by the time when the
system enters unstable evolution and the maxima of |fz|
tend to move towards the minima of mT (Fig. 3). More
quantitatively, we performed Nr evolutions starting from
different initial conditions, but fixed σ. As in the exper-
iment [5], we average over these runs to wash out shot-
to-shot fluctuations. In Fig. 4 we plot the mean domain
size: ξ =
i ξz(i)/Nr, where i = 1, ..., Nr and ξz(i) is
the mean domain size in the i-th run. As shown in Fig.
4a, for t . t̂ we observe ξ ≈ f(t/τ1/3Q ) as for MT (t). The
domains are formed on a time scale of ∼ t̂. A simple anal-
ysis based on KZ theory [1, 2] suggests that their charac-
teristic post-transition size, ξ̂, should be roughly given by
∫ ∼t̂
dtvs(t), where vs(t) is a sound velocity. There are two
sound modes in the broken-symmetry phase that prop-
agate both spin and density fluctuations [13]: the faster
(slower) one has velocity ∼ √c0 (∼
q|c1|). Putting any
of these as vs into the integral, and assuming τQ ≫ 1 for
the slower mode, we get
ξ̂ ∼ τ1/3Q . (5)
This result correctly predicts the scaling property of
the size of post-transition “defects” as is evident from
the overlap of different curves in Fig. 4, which shows
up for τQ ≥ 25 or 0.9s. Quantitatively, we define ξ̂
as the value of ξ averaged over q ∈ [1/2, 1] to wash
out post-transition fluctuations. A fit got us ln ξ̂ =
(−0.38 ± 0.03) + (0.30 ± 0.01) ln τQ, in good agreement
with (5). The fit was done to τQ ≥ 30 data and is pre-
sented in Fig. 4b which illustrates that smaller τQ data
gradually departs from 1/3 scaling law.
Now we focus on the analytical calculations provid-
ing predictions about early stages of time-evolution.
We assume that the wave-function stays close to
the polar phase GS, ΨT = (δψ1(t), 1/
δψ0(t), δψ−1(t)) exp(−iµt), where the chemical potential
is µ = c0/L, |δψm| ≪ 1/
L, and
dz(δΨ0 + δΨ
0) ≡ 0
to keep
dzΨ†Ψ = 1 + O(δΨ2). Linearizing the cou-
pled nonlinear-Schrödinger equations that describe the
system we get fχ = ReGχ, where χ = x, y, Gx =√
2(δΨ1 + δΨ−1)/
L, Gy = i
2(δΨ1 − δΨ−1)/
L, and
Gχ = −
qGχ −
(Gχ +G
where α = 2|c1|/L. To solve this equation we go to
momentum space, aχ(k) =
dzfχ exp(ikz) and bχ(k) =
dzImGχ exp(ikz), getting
0 k2 + αq
2α− k2 − αq 0
. (6)
Diagonalizing the matrix from Eq. (6) we see that there
is instability for k2/α < 2 − q as in the Bogolubov spec-
trum of this model [13]. Thus, the system is stable in
the polar phase, and so small initial perturbations do
not grow during the evolution towards broken-symmetry
phase. The instability for q < 2 is responsible for the
magnetization jump in Fig. 1 and the subsequent break-
down of the linear approach (Fig. 1b).
To solve Eq. (6) with q(t) given by (2) we derive the
equation for d2aχ(t)/dt
2, keep leading order terms in the
slow transition (τQ ≫ 1) and long-wavelength (k2/α ≪
2) limits, and get that
aχ(k, t) = αkχAi(s)+βkχBi(s),
, (7)
0 0.5 1 1.5 2
q(t)=2-t/τ
0.75 1.5
4 16 64
FIG. 4: (color) Dynamics of magnetic domains in fz. Black
line (τQ = 30), red line (τQ = 50), green line (τQ = 80). The
arrow show direction of evolution on the main plot. Inset (a):
early stages of ξ(t) evolution. Inset (b): dependence of the
typical post-transition domain size, ξ̂, on quench time. The
fit was done to τQ ≥ 30 data (see text for details). The figure
shows results averaged over Nr = 44 runs; see [12] for units.
where κ = (α2/2)1/3, αkχ and βkχ are constants given
by initial conditions, while Ai and Bi are Airy func-
tions. From (7) we see that the instability arises from
unbounded increase of the Bi(s) function happening for
s > 0, i.e., k2/α < 2−q(t), which is a dynamical manifes-
tation of the static result for unstable modes. This solu-
tion works till t ∼ t̂ ∼ τ1/3Q when a significant increase of
fχ invalidates the linearized theory: this calculation rig-
orously derives scaling (3). Additionally, the solution (7)
can be reliably used as long as τQ ≫ 1 or 37ms, which
is also supported by numerics (Fig. 1a). The quench
time scale in the experiment [5] is much smaller than
this bound. Finally, these results hold for any initial
state spread over the k modes.
The (re)scalings t/τ
Q and ξ̂ ∼ t̂ ∼ τ
Q derived above
in a 1D system were also found by different means in a 2D
spinor condensate [15]. A trivial extension of our mean-
field analytical calculations to 2D and 3D systems shows
that they hold for any number of spatial dimensions.
To summarize, we have developed a theory of the dy-
namics of symmetry-breaking in the quantum phase tran-
sition inspired by the experiment [5], but for the range of
quench rates that are sufficiently slow so that the critical
scalings can determine phase transition dynamics. This
regime should be accessible by a “slower” version of the
quench [5]. Our analysis points to a Kibble-Zurek-like
scenario, where the state of the system departs from the
old symmetric vacuum with a delay ∼ t̂ after the critical
point was crossed. This sets up an initial post-transition
state with a characteristic length scale ξ̂ (5). This scale
should determine the initial density of topological fea-
tures. In our 1D simulations textures appear, but we
predict that in real 3D experiments other topological de-
fects are created (as they were in [5]), and the distance
between them should be initially ∼ ξ̂. Such topological
defects are more stable than textures so measurement of
their density should be possible and would be a good test
of the theory we have presented.
[1] T.W.B. Kibble, J. Phys. A 9, 1387 (1976); Phys. Rep.
67, 183 (1980).
[2] W.H. Zurek, Nature (London) 317, 505 (1985); Acta
Phys. Pol. B 24, 1301 (1993); Phys. Rep. 276, 177 (1996).
[3] I. Chuang et al., Science 251, 1336 (1991); M.J. Bowick
et al., ibid. 263, 943 (1994); C. Bauerle et al., Nature
(London) 382, 332 (1996); V.M.H. Ruutu et al., ibid.
382, 334 (1996); S. Ducci et al., ibid. 83, 5210 (1999); A.
Maniv, E. Polturak, and G. Koren, Phys. Rev. Lett. 91,
197001 (2003); R. Monaco et al., ibid. 96, 180604 (2006).
[4] M. Lewenstein et al., Adv. Phys. 56, 243 (2007).
[5] L.E. Sadler, J.M. Higbie, S.R. Leslie, M. Vengalattore,
and D. M. Stamper-Kurn, et al., Nature (London) 443,
312 (2006).
[6] S. Sachdev, Quantum Phase Transitions (Cambridge
University Press, Cambridge UK, 2001).
[7] W.H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett.
95, 105701 (2005).
[8] J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005); R.W.
Cherng and L.S. Levitov, Phys. Rev. A 73, 043614
(2006).
[9] B. Damski, Phys. Rev. Lett. 95, 035701 (2005). B.
Damski and W.H. Zurek, Phys. Rev. A 73, 063405
(2006).
[10] T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998); T. Ohmi and
K. Machida, J. Phys. Soc. Jpn. 67, 1822 (1998).
[11] T.P. Meyrath et al., Phys. Rev. A 71, 041604(R) (2005).
[12] The experiment [5] is done in the harmonic potential
mω2z(λ
2 + λ2yy
2 + z2)/2, where m is the 87Rb mass,
ωz = 2π×4.3Hz, λx = 13, and λy = 81.4. We use the os-
cillator units through the paper for time (1/ωz = 37ms),
length (
~/mωz = 5.1µm), and energy (~ωz), and as-
sume that the system stays in the harmonic oscillator
ground states in x and y directions. Then, after skip-
ping constant terms, the 3D energy functional [10] re-
duces to dimensionless (1) with the additional Ψ†Ψz2/2
term, and ci = 2Nαi
mωz/~
λxλy/3 (N = 2× 106 is
the atom number, α0 = 16nm, and α1 = −α0/216.1).
To approximate such a system by a box we assume that
the total density of atoms in the harmonic trap center
is the same as in the box of size L. Neglecting in the
Thomas-Fermi limit the first and the last term of (1) we
get L = 2(
c02/3)
2/3 ≈ 78, i.e., 78× 5.1µm ≈ 0.4mm.
[13] K. Murata, H. Saito, and M. Ueda, Phys. Rev. A 75,
013607 (2007).
[14] H. Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. A 75,
013621 (2007).
[15] A. Lamacraft, Phys. Rev. Lett. 98, 160404 (2007).
[16] M. Uhlmann, R. Schützhold, and U.R. Fischer,
cond-mat/0612664.
[17] G.J. Stephens, Phys. Rev. D 61, 085002 (2000).
[18] A. Vilenkin and E.P.S. Shellard, Cosmic strings and
other topological defects (Cambridge University Press,
Cambridge, 1994).
[19] H. Mäkelä, Y. Zhang, and K.-A. Suominen, J. Phys. A
36, 8555 (2003).
http://arxiv.org/abs/cond-mat/0612664
