Quantum phase transition in space in a ferromagnetic spin-1 Bose-Einstein condensate

A quantum phase transition between the symmetric (polar) phase and the phase with broken symmetry can be induced in a ferromagnetic spin-1 Bose-Einstein condensate in space (rather than in time). We consider such a phase transition and show that the transition region in the vicinity of the critical point exhibits scalings that reflect a compromise between the rate at which the transition is imposed (i.e., the gradient of the control parameter) and the scaling of the divergent healing length in the critical region. Our results suggest a method for the direct measurement of the scaling exponent nu.


I. INTRODUCTION
Studies of phase transitions have traditionally focused on equilibrium scalings of various properties near the critical point of a homogeneous system. The dynamics of phase transitions presents new interdisciplinary challenges. Nonequilibrium phase transitions may play a role in the evolution of the early Universe [1]. Their analogues can be studied in condensed matter systems [2]. The latter observation led to the series of beautiful experiments [3] (see [4] for an up-to-date review) and to the development of the theory based on the universality of critical behavior [2]. The recent progress in cold atom experiments allows for the temporal and spatial control of different systems undergoing a quantum phase transition (QPT) [5,6]. These experimental developments call for an in-depth understanding of non-equilibrium QPTs.
A QPT is a fundamental change in the ground state of the system as a result of small variations of an external parameter, e.g., a magnetic field [7]. In contrast to thermodynamic phase transitions, it takes place ideally at temperature of absolute zero.
The problem of how a quantum system undergoes a transition from one quantum phase to another due to timedependent (temporal) driving has attracted lots of attention lately [8,9,10,11,12,13,14,15]. Basic insights into the QPT dynamics can be obtained through the quantum version [8,16] of the Kibble-Zurek mechanism (KZM) [1,2]. The KZM recognizes that the time evolution of the quantum system is adiabatic far away from the critical point where the gap in the excitation spectrum is large. The system adjusts to all changes imposed on its Hamiltonian. Near the critical point, however, the gap closes precluding the adiabatic evolution and resulting in the non-equilibrium dynamics. This switch of behavior happens when the system reaction time /∆ (∆ is excitation gap) becomes comparable to the time scale on which the critical point is approached, ε/|dε/dt| (ε = |q − q c |, q is the parameter driving the transition, q c is the location of the critical point). This brings us to To solve (1) we define d dt q(t) = τ −1 Q . The solution of (1) gives us timet, left to reaching the critical point, when the non-equilibrium dynamics startst ( Above we have set ∆ = ∆ 0 |q c − q| zν near the critical point (z and ν are critical exponents). The quantum KZM proposes that the non-equilibrium evolution in the neighborhood of the critical point is to a first approximation impulse (diabatic): the state of the system does not change there. Suppose that the system evolves towards the critical point from the ground state far away from it. Thus its properties near the critical point are determined by its ground state properties at a distancet/τ Q from the critical point (the adiabatic regime was abandoned there). As have been shown in spin models (e.g., [8,9]) and ferromagnetic spin-1 Bose-Einstein condensates [11,12,13], this simplification allows for a correct analytical computation of the density of excitations resulting from the non-equilibrium quench.
In this paper we study a quantum phase transition induced by spacial rather than temporal driving. It means that the driving parameter, q, depends on position and is independent of time. Such a transition was recently studied in the quantum Ising model [17,18] and the mean-field Ginzburg-Landau theory [17]. Here we present the theory of the spatial quench in a ferromagnetic spin-1 Bose-Einstein condensate, which is one of the most flexible physical systems for studies of QPTs.
We focus on the ground state. In the simplest approximation, the local homogeneous approximation (LHA), the system is locally characterized by its homogeneous ground state properties perfectly tracking spatial variations of q.
This approximation is good far enough from the critical point, where the healing length [19] is small compared to the imposed scale of spatial driving. Using the time quench analogy, the system evolves perfectly adiabatically in space away from the spatial critical point, r c , where q(r c ) = q c .
Around the critical point, however, the healing length diverges and so the LHA breaks down. The system enters "non-equilibrium" regime similarly as during a time quench. The spatial coupling term in the Hamiltonian (a gradient term in mean-field theories, a spin-spin interaction term in spin models) smoothes out the sharp spatial boundary between the phases predicted by the LHA.
Even more interestingly, these similarities are not only qualitative [18]. The distance from the spatial critical point where the switch between the two regimes takes place,x, can be determined by looking at the length scales relevant for the spatial quench. We compare the healing length to the length scale on which the spatial critical point is approached: ε/|dε/dx| assuming that ε = |q − q c | changes in x-direction only. This brings us to the spatial analog of (1) To solve (4) we define d dx q(x) = λ −1 Q and getx = ξ The results (2) and (5) are analogous. Indeed, apart from a slight difference in the scaling exponents -due to the absence of z in (5) -t maps ontox when time scales /∆ 0 and τ Q are replaced by length scales ξ 0 and λ Q , respectively. This shows striking parallels between quenches in time and space. It is also instructive to note that the same general result, Eq. (5), can be obtained in a different way from the scaling theory [17].

II. MODEL
In the following, we will study a QPT in space in a ferromagnetic spin-1 Bose-Einstein condensate. We consider untrapped clouds: atoms in a box. Such a system can be realized in an optical box trap [20]. Assuming that the condensate is placed in the magnetic field B(r) aligned in the z direction, its mean-field energy functional reads [21] where F α=x,y,z is the spin-1 matrix. The strength of spin independent interactions is c 0 = 4π 2 (a 0 + 2a 2 )/3M > 0, while the strength of spin dependent interactions reads c 1 = 4π 2 (a 2 − a 0 )/3M < 0. The constant a S is the s-wave scattering length in the total spin S channel (a 0 = 101.8a B , a 2 = 100.4a B ), and M is the atom mass. The prefactors of linear and quadratic Zeeman shifts are given by P = µ B B(r)/2 and Q = µ 2 B B 2 (r)/4E hf , respectively. There E hf is the hyperfine splitting energy.
The condensate magnetization reads It is convenient to define a dimensionless parameter where n ≈ N/V is the condensate density (V is the system volume). We consider below setups where n fluctuates negligibly in space. The critical point corresponds to q c = 2.  (7)]. This numerical simulation is for λQ/λ0 = 3/2 -see Appendix B for numerical details.

III. QUENCH IN A FERROMAGNETIC SPIN-1 BOSE-EINSTEIN CONDENSATE
Similarly as in our work on time-dependent quench [11,12], we assume that drf z = 0. Experimentally, such a constraint can be achieved by putting all atoms into an equal mixture of m = ±1 magnetic sublevels and letting the sample relax to equilibrium [21]. The spin conservation ensures that the final and initial states will have the same total f z magnetization: drf z = 0. Additionally, we consider ψ ±1,0 ≥ 0 in the ground state: a condition that can be . This sets f y (r) = 0. Moreover, our numerical simulations indicate that f z (r) ≈ 0 in the spatial quench considered here. Thus we investigate the f x (r) condensate magnetization only.
For simplicity, we assume that and forgetting about the gradient term in (6) -i.e., using our LHA -the ground state phase diagram can be sketched. The part of the condensate exposed to 0 ≤ q(r) < 2 is in the ground state of the broken-symmetry phase. There ψ GS ±1 (r) = √ n 1/4 − q(r)/8 and ψ GS 0 (r) = √ n 1/2 + q(r)/4. The condensate is magnetized: f x (r) = n 1 − q(r) 2 /4, while f y,z (r) = 0. This ground state breaks rotational symmetry on the (x, y) plane present in energy functional (6). Parts of the condensate exposed to q(r) > 2 are in the polar phase ground state with ψ GS ±1 (r) = 0 and ψ GS 0 (r) = √ n, which implies f x,y,z (r) = 0. The dependence of the condensate magnetization on x and q in the LHA is depicted with the red dashed line on Fig. 1. There one condensate accommodates the two phases. The inclusion of the gradient term in the energy functional (6) prohibits singularities in the first derivative of the wave function. The smooth crossover region replaces the sharp phase boundary predicted by the LHA (Fig. 1). Its size is proportional tox (5), as will be carefully discussed below.
In the simplest setup q(r) dependence (7) is achieved by using the static, inhomogeneous, magnetic field to drive the system between the two phases. That would correspond to the following spatial variation of the Zeeman coefficients across a typical 87 Rb condensate: ∆Q ∼ n|c 1 | and ∆P ∼ 10 4 n|c 1 |. The latter result is obtained by using P = h × 70 × 10 4 × B HzG −1 , Q = h × 70B 2 HzG −2 [6], and n|c 1 |/h = 9.77 (evaluated at a peak density of the Berkeley experiment [6]). As can be found numerically, such an abrupt P -variation makes the LHA phase diagram qualitatively incorrect.
To avoid these complications we propose to expose the condensate to the inhomogeneous magnetic field B pointing in the z-direction and changing in time faster than the characteristic time scales for the condensate dynamics. When B(x, t) t = 0, where · · · t denotes time average, the linear Zeeman shift will be wiped out (which we assume from now on), while the quadratic term will be position-dependent only: Q(x) ∼ B(x, t) 2 t . Thus, we are still dealing with a phase transition in space.
To make use of the general predictions worked out in (5), we have to determine the value of the critical exponent ν and the prefactor ξ 0 from the generic expression for the healing length (3). We do it by comparing (3) we also see from (A3) that ξ 0 equals ξ s in the polar phase and ξ s / 2 − 2|c 1 |/c 0 in the broken-symmetry phase. Coming back to (5) we get that the size of the crossover region on either the polar or the broken-symmetry side iŝ This expression can be verified by looking at experimentally measurable quantity: the ground state magnetization of the condensate [6]. We start our discussion on the polar side and linearize the wave-function around polar phase ground state: This leads to f x = √ 2n(δψ 1 +δψ −1 )+O(δψ 2 m ). Using position-dependent q from (7) and linearized version of (A1) we arrive at which is solved by There f crit x is the condensate magnetization at the spatial critical point (it depends on λ Q only), ∆x = x − 2λ Q ≥ 0 is the distance from the spatial critical point, and Ai is an Airy function -the only nondivergent solution of (10). This solution rigorously shows that decay of the magnetization, f x (r)/f crit x , takes place on a length scale (λ Q ξ 2 s ) 1/3 in full agreement with the strikingly simple scaling result (9). Therefore, in the polar phase, the magnetization approaches the LHA value, f x = 0, at a distancex from the critical point.
On the broken-symmetry side of the QPT, we refer to numerics (Appendix B) to verify accuracy of (9). In the two extreme limits, very large and very small λ Q the theory does not work well. Indeed, using (7) we see that the spatial extent of the broken-symmetry phase is 2λ Q (Fig. 1). Thus, we expect thatx ≪ 2λ Q , i.e., λ Q ≫ ξ 0 /2 3/2 , to see well the crossover region on the broken-symmetry side. In the other limit, large λ Q , we have to take into account that the system size in numerical simulations is finite, say l. There we need to have λ Q ≪ l/2 so that there will be a spatial point in our system where q(x) = 2, and the finite size corrections coming from proximity of the system boundary to the spatial critical point will be negligible. These two estimates tell us that λ Q in our system shall be much larger than 0.1λ 0 and much smaller than 39λ 0 (see Appendix B for numerical value of l, ξ 0 ≈ ξ s / √ 2 and the unit of length λ 0 relevant for our simulations). Increasing (decreasing) the lower (upper) limit on λ Q by a factor of five we fit a power law to λ Q /λ 0 ∈ (1/2, 8). We get from a fit that lnx/λ 0 = −1.7984 ± 0.0005 + (0.3082 ± 0.0004) ln λ Q /λ 0 -see Fig. 2. The numerical data forx clearly departs from the fitted line for λ Q /λ 0 10. Departures for λ Q /λ 0 1/2 are much less pronounced. Therefore, we get approximately thatx ∼ λ 0.31 Q , which is in good qualitative agreement with the predicted scaling, i.e., λ 1/3 Q . We expect that better agreement can be obtained for large systems where we can explore larger λ Q 's for which crossover region is closer to the critical point (x/λ Q ∼ λ −1/3 Q ). More precisely, by taking the limit of N, l, λ Q → ∞ at ξ 0 = const (i.e., N/l = const) we expect that the scaling (9) will be fully recovered from numerical simulations. Now let's focus on the magnetization at the critical point. From (11) we know that it scales in the same way as the condensate magnetization at the border of the crossover region in the polar phase (x + = 2λ Q +x): Assuming that similar scaling relation will also hold between f crit x and f x at x − = 2λ Q −x, i.e., at the broken-symmetry phase border of the crossover region where we get that f crit Fig. 2, we have indeed a robust f crit scaling for λ Q /λ 0 in the large range of 10 −1 to 10. It is quite remarkable that the above intuitive prediction of the f crit x scaling matches numerics so well. It suggests that magnetization at the spatial critical point might be a robust observable for studies of spatial quenches. Similar observation was made for condensate magnetization at the critical point during a temporal quench from the broken-symmetry to the polar phase [12].
All these results are in qualitative agreement with our expectations: in the limit of large λ Q -very slow spatial driving -the condensate ground state approaches the LHA. Indeed, both the size of the crossover region in the q-parameter space,x/λ Q , and the condensate magnetization at the critical point, f crit x , go to zero as λ Q increases. The measurements of the dependence of the size of crossover region on λ Q should allow for the first experimental determination of the scaling exponent ν in a ferromagnetic spin-1 Bose-Einstein condensate. Thus, our considerations lead to the new way of investigating the critical region of phase transitions. Another scheme that can be used to determine the exponent ν was explored experimentally by Esslinger's group in the context of the classical phase transition from a normal gas to a Bose-Einstein condensate [22]. The setup involved in this experiment is quite complicated as it requires presence of a high finesse optical cavity for detection purposes. We expect that the approach proposed by us should be easier to implement in the context of quantum phase transitions. Moreover, according to (5) our scheme shall be also ready for experimental exploration in the whole zoo of other physical systems undergoing a QPT. Finally, we would like to stress that applicability of Eq. (5) is not limited to mean-field theories: any experimental departures from the λ 1/3 Q scaling (apart from finite size corrections) should indicate that mean-field value of the critical exponent ν does not describe the condensate properly. In such a case the correct value of the ν exponent should be easily extracted from experimental data with the help of Eq. (5).

IV. SUMMARY
We have explored physics of the QPT in space in a ferromagnetic spin-1 Bose-Einstein condensate [6]: the most flexible system for studies of QPTs in space and time. Our scaling results explain how singularities of the critical point affect the ground state magnetization of the condensate. They also suggest a new way for the measurement of the critical exponent ν. These findings are generally applicable to all systems undergoing a second order QPT. The quest for the full understanding of a spatial quench opens up a prospect of interdisciplinary studies similarly as time quenches have done to date [4].

V. ACKNOWLEDGMENTS
We acknowledge the support of the U.S. Department of Energy through the LANL/LDRD Program.

APPENDIX A: HEALING LENGTH OF A FERROMAGNETIC SPIN-1 CONDENSATE
The key ingredient of our theory is the divergent healing length (3). Below we derive it from mean-field equations that come from minimization of E − µ dr m |ψ m | 2 without the linear Zeeman term (6) -see Sec. III for the explanation of why we remove it from our considerations. The healing length is a typical length scale over which a local perturbation of the wave-function gets forgotten (here we will have three characteristic length scales, but we will identify the leading one relevant for a long distance healing process). To find it, we assume the constant parameter q = Q/n|c 1 | across the condensate, and linearize mean-field equations (A1). We write the wave-function as ψ m = ψ GS m + √ nδψ m with m = 0, ±1 and ψ GS m being the ground state solution. In the broken-symmetry phase we have ψ GS ±1 (r) = √ n 1/4 − q(r)/8 and ψ GS 0 (r) = √ n 1/2 + q(r)/4, while in the polar phase, ψ GS ±1 (r) = 0 and ψ GS 0 (r) = √ n. We use here the freedom to work with ψ 0,±1 ≥ 0, which is explained in Sec. III.
Linearized equations take the form where ξ s = 2 /2M n|c 1 | and S stands for in the broken-symmetry phase, and for in the polar phase. The factor u ± above is very close to unity as |c 1 |/c 0 ≈ 1/216.1 ≪ 1 for 87 Rb atoms considered here.
To proceed we diagonalize matrix S by the transformation S = CDC −1 , where D = diag(Γ 1 , Γ 2 , Γ 3 ) is a diagonal matrix made of eigenvalues of matrix S. In the polar phase we find Γ m = q − 2, q, 2c 0 /|c 1 |, while in the broken-symmetry phase we get After defining δψ m = n C −1 mn δψ n and some elementary algebra we get where the form of the function f depends on system dimensionality (1D, 2D, or 3D), while r is "attached" to the point in space where the perturbation is imposed on the wave-function. This solution implies that the eigenvalue vanishing at the critical point, Γ min , provides the leading contribution to the long-distance healing process near the critical point (q ≈ 2). Indeed, a simple calculation shows that the wave function will forget about the perturbation imposed on it (δψ m , δψ m ≈ 0) at a distance scaling as ξ s / √ Γ min . Thus, the divergent healing length equals ξ s / √ Γ min . Around the critical point we find it to be ξ broken = ξ s 2 − 2|c 1 |/c 0 |q − 2| 1/2 , ξ polar = ξ s |q − 2| 1/2 (A3) on the broken-symmetry (q < 2) and polar (q > 2) sides of the QPT, respectively. The expansion around q = 2 was applied to obtain ξ broken , i.e., we have approximated the eigenvalue c 0 /|c 1 | − (c 0 /c 1 ) 2 − (4 − q 2 )(c 0 /|c 1 | − 1) by (2 − 2|c 1 |/c 0 )(2 − q) near the critical point. These expressions fix the value of the mean-field critical exponent ν for our model: ν = 1/2 (compare (A3) to (3)). We conclude by showing another simplification of the full expression for the divergent healing length on the broken symmetry side. Taking the limit of |c 1 |/c 0 → 0, well justified for 87 Rb, we get without restricting ourselves to q ≈ 2.

APPENDIX B: DETAILS OF NUMERICAL SIMULATIONS
Our analytical predictions are valid in any dimensional system, but we test them numerically in a 1D model. Our numerics is based on conjugate gradient minimization of the 1D energy functional obtained through dimensional reduction of the full 3D energy functional (6) without the linear Zeeman term -see Sec.
III for the explanation of why we remove it from our considerations. Moreover, all the components of the above energy functional are dimensionless, and chosen to match qualitatively the experiment of the Berkeley group [6]. Derivation of (B1) is explained in detail in [11]: please account for different normalization of the wave-function between this paper and [11].