How to fix a broken symmetry: Quantum dynamics of symmetry restoration in a ferromagnetic Bose-Einstein condensate

We discuss the dynamics of a quantum phase transition in a spin-1 Bose-Einstein condensate when it is driven from the magnetized broken-symmetry phase to the unmagnetized ``symmetric'' polar phase. We determine where the condensate goes out of equilibrium as it approaches the critical point, and compute the condensate magnetization at the critical point. This is done within a quantum Kibble-Zurek scheme traditionally employed in the context of symmetry-breaking quantum phase transitions. Then we study the influence of the nonequilibrium dynamics near a critical point on the condensate magnetization. In particular, when the quench stops at the critical point, nonlinear oscillations of magnetization occur. They are characterized by a period and an amplitude that are inversely proportional. If we keep driving the condensate far away from the critical point through the unmagnetized ``symmetric'' polar phase, the amplitude of magnetization oscillations slowly decreases reaching a non-zero asymptotic value. That process is described by the equation that can be mapped onto the classical mechanical problem of a particle moving under the influence of harmonic and ``anti-friction'' forces whose interplay leads to surprisingly simple fixed-amplitude oscillations. We obtain several scaling results relating the condensate magnetization to the quench rate, and verify numerically all analytical predictions.


I. INTRODUCTION
Any symmetry breaking phase transition can be traversed in two opposite directions. The usual way is to start in the symmetric phase, and move into the phase with a broken symmetry. When this process happens on a finite timescale, critical slowing down will generally lead to the random local choices of the broken-symmetry vacua. This in turn can result in the formation of topological defects [1,2]. They will appear with the density set by the size of the regions that choose the same broken-symmetry vacuum. That size -and, hence, the density of the resulting non-equilibrium structures -can be deduced from the critical scalings of the relaxation time and the healing length.
The obvious question that has not been addressed to date concerns traversing the transition from the phase in which the symmetry is broken to where it is restored. Such transition can be also induced on a finite timescale, which will again bring into the discussion the critical scalings. Now, however, they will determine remaining excitation of the system rather than the size of broken-symmetry regions.
We study here the dynamics of quantum phase transitions [3,4,5,6] rather than their classical (thermodynamical) counterparts that have received extensive theoretical [1,2] and experimental [7,8] attention (see [9] for a recent review). More precisely, we are interested in the dynamics of a ferromagnetic Bose-Einstein condensate (BEC) composed of spin-1 atoms [10]. Such a condensate is studied experimentally by several groups [11,12]. Similar issues may also arise in an even more complex case of higher spin condensates [13]. The research on symmetry breaking in the spin-1 BEC has been the subject of the seminal paper of the Berkeley group [11]. There the condensate was rapidly driven from the polar (symmetric) to the ferromagnetic (broken symmetry) phase by the decrease of a magnetic field. Formation of topological defects (e.g., vortices) was observed. This work has triggered several theoretical investigations [14,15,16,17,18,19]. In particular, vortex density after a slow quench was found through the extensive numerical simulations by Saito, Kawaguchi and Ueda [18] to follow the scaling result obtained from the quantum Kibble-Zurek mechanism [17] and the studies of the magnetization correlation functions [15,18].
Below we consider dynamics of a ferromagnetic condensate driven through a transition "in the opposite direction": from the broken-symmetry phase to the polar phase. As in [17], we are interested in slow transitions that exhibit adiabatic-impulse behavior so that the condensate goes out of equilibrium close to the critical point. In this case the condensate excitation in the polar phase reflects the scalings of the critical regime. This scenario is supported by the following general discussion. Suppose that we change some dimensionless parameter q to drive the system towards the critical point at q c . Close to the phase boundary, the gap ∆ in the excitation spectrum behaves as ∆ 0 (q c − q) zν , where z and ν are critical exponents [20], q < q c , and ∆ 0 is a constant with dimension of energy. As long as the system is driven far enough from the critical point, the gap is large and the system excitation does not occur. The relevant "reaction time" for the system is given by The system undergoes adiabatic evolution when τ is small compared to the timescale on which changes occur in the Hamiltonian. That timescale, in turn, is characterized by how fast the "instantaneous" excitation gap changes, or in other words by how fast the system is driven. It is simply given by having a proper dimension of time. Comparison of the timescales (1) and (2), gives the location of the border between adiabatic -non-adiabatic behavior at some q = q c −q. Assuming, that d dt q(t) = τ −1 Q near a critical point (τ −1 Q is the speed of transition between the phases) we get that ( The qualitative picture of system dynamics is as follows. For q c − q(t) q the system undergoes adiabatic evolution following instantaneous changes to its Hamiltonian. When q c − q(t) q the evolution ceases to be adiabatic and to a first approximation one can expect the impulse stage, where the state of the system does not change -remains "frozen". As a result, system properties near the critical point are determined by its ground state at q = q c −q, whereq is computed from the quench time τ Q and the product of the critical exponents: Eq. (3). This approach was already successfully used to study the Landau-Zener dynamics [5], the dynamics of the quantum Ising model [3], and the ferromagnetic condensate dynamics during the symmetry-breaking transitions from the polar to the broken-symmetry phase [17]. We begin our study in the next section, where we discuss the basics of a quantum phase transition in the mean-field model that we consider. Section III presents how the condensate gets excited when driven through the brokensymmetry phase. Section IV is devoted to studies of magnetization oscillations after quench. We point out the observables that contain the information about condensate excitation at the phase boundary. Section V determines the range of applicability of the Single Mode Approximation (SMA) used in Sections III and IV. Finally, Section VI provides the summary of the paper.

II. THE MODEL
The energy of the ferromagnetic condensate placed in an external, homogeneous, magnetic field aligned in the z direction is given in the mean-field approximation as where F α is the standard spin-1 matrix, and c 0 > 0 and c 1 < 0 provide the strength of spin independent and spin-dependent interactions, respectively: with a S being the s-wave scattering length in the total spin S channel, M is the atom mass, while Q (proportional to square of magnetic field imposed on the condensate) is the prefactor of the quadratic Zeeman shift resulting from atom-magnetic field interactions. The linear Zeeman term is skipped in Eq. (4) because it has a trivial effect on condensate dynamics. Indeed, it only rotates the condensate magnetization around the z axis with the Larmor frequency. This rotation, of course, can be removed from theoretical studies by going to a reference frame rotating with the Larmor frequency (Sec. II.A of Ref. [14]). The wave function has three condensate components, ψ m , corresponding to m = 0, ±1 projections of spin-1 onto the magnetic field: where N is the total number of atoms. In subsequent calculations we will use the phases χ m defined as: ψ m = |ψ m | exp(iχ m ). The simplest meaningful approach to the dynamics of a ferromagnetic Bose-Einstein condensate is provided by the Single Mode Approximation where the translationally invariant mean-field approach is used. This simplification, while obviously inapplicable for symmetry-breaking phase transitions where the orientation of magnetization is explicitly position-dependent [17,18], can be successfully used for the transitions considered here where the system starts from a uniform broken-symmetry ground state (GS). The influence of inhomogeneities will be illustrated in Sec. V where we will compare simulations done with the SMA to full mean-field calculations. In the SMA both the gradient term and the density-density interaction term ∼ c 0 are irrelevant and so are skipped in Secs. II-IV. These terms, however, reappear in the numerical calculations in Sec. V.
As the ferromagnetic condensate is described in the SMA by translationally invariant spinor (5), it is instructive to parametrize the coupling constant Q from (4) by where n = Ψ † Ψ is the atom density, and q is a dimensionless parameter proportional to the square of the magnetic field. This parameter will be changed to drive the condensate from one quantum phase to another.
The system longitudinal (f z ) and transverse (f x and f y ) magnetizations read We consider the f z = 0 case, i.e., |ψ 1 | = |ψ −1 |. This condition is dynamically conserved during all the evolutions considered in Secs. III and IV because d dt drf z (r, t) = 0, which in the SMA reduces to d dt f z (t) = 0. The f x,y components are conveniently combined to a complex transverse magnetization The quantum phases of the ferromagnetic condensate were recently discussed in [21]. The condensate is in the polar phase when q > 2. There the ground state is given by so that f x = f y = 0. For 0 ≤ q < 2 the system is in the broken-symmetry phase where the GS spinor reads and The system dynamics in the SMA is governed by [22] i where d/dt is denoted with a "dot" and we have used the identity |ψ 1 | = |ψ −1 | to arrive at Eqs. (9). The initial condition for time evolutions in Secs. III and IV is provided by the spinor state (7) taken at q = 0. The last equation of (9) states that the orientation of the transverse magnetization on the (x, y) plane is conserved throughout the evolution and so it is fixed by the initial conditions. This constraint prohibits full restoration of the symmetry in the state evolved from the broken-symmetry phase to the symmetric polar phase. In particular, We will discuss below the dynamics driven by the variation of q that can be achieved with a proper manipulation of the magnetic field imposed on the condensate. A linear increase of the magnetic field strength in time is well within the experimental capabilities [11]. Due to the quadratic Zeeman coupling, it results in This is the time dependence we will assume in all subsequent calculations. Close to the phase transition point, where t = 4τ Q − δt, we have For the transitions considered here the system will cease to adiabatically follow the ramp up of q near a critical point where the above linearization holds: the familiar time-dependence from the Kibble-Zurek theory emerges near the phase boundary with τ −1 Q providing the quench rate [23], or in other words, the speed in the parameter space of driving the system through the critical point.

III. SYSTEM EXCITATION
For small perturbations around the GS of the broken-symmetry phase one finds three Bogolubov modes as in [21]: two gapless modes and one gapped. We do not consider analytically the dynamics induced by the excitation of the gapless modes. Instead, we focus on the gapped mode and show that its excitation is responsible for changes of the magnitude of the system transverse magnetization during non-equilibrium dynamics. This mode in the long wavelength limit, obviously relevant for the SMA, has the energy gap [21] As we drive the system to the critical point changing q(t) we expect qualitatively the following dynamics. First, for q(t) small enough the evolution is adiabatic with respect to the gapped mode -far enough from the critical point the energy cost of exciting the gapped mode is too large. However, as the critical point is approached, the gap becomes too small to allow for adiabatic evolution and the system starts to populate the gapped mode: the non-adiabatic dynamics starts. Below we assume that the separation of the two regimes takes place at q = 2 −q, and determine the location ofq as a function of the quench rate. This can be done as in Sec. I. Comparison of the timescales (1) and (2) brings us to the (approximate) equation This equation is illustrated in Fig. 1. Its solution in the slow transition limit,q ≪ 1, giveŝ This, of course, corresponds to (3) with zν = 1/2 and ∆ 0 = 2n|c 1 |. A simple estimation of the constant τ 0 by the peak density in the Berkeley experiment, 2.8 × 10 14 cm −3 , gives τ 0 ≈ 16ms. In this section and Sec. IV all the results can be obtained without assuming any particular value of τ 0 . Qualitatively, prediction (13) is in agreement with what we expect: the slower the quench is (the larger τ Q is), the closer to the critical point the condensate can be driven adiabatically.
Quantitatively, solution (13) turns out to be remarkably accurate in the wide range of quench times τ Q as depicted in Fig. 2. In that plot we define the instantq to be the distance from a critical point when the departure of system transverse magnetization from a static (ground state) prediction starts exceeding some threshold, e.g., 5% of the instantaneous GS value. From the fit there we get that  |fT | 2 static refers to the squared modulus of (8). We arbitrarily define that when these departures exceed 5%, the system leaves the adiabatic regime and starts the non-equilibrium dynamics (other reasonable thresholds < 10% give the same scaling result). The inset shows the result of power law fit pointing toq = 1.05 × (τ0/τQ) 0.651±0.001 , which compares well to theoretical prediction (13). The fit is presented as a solid line, while numerics comes as pluses. Note the log-log scale on the inset. The fit was done for τQ/τ0 = 30 · · · 200. The larger τQ's are taken to the fit, the closer scaling exponent approaches 2/3. On the main plot we have a simulation for τQ/τ0 = 30.
in excellent agreement with (13) with respect to the scaling exponent. The prefactor, which obviously depends on the threshold used for the numerical determination ofq is of the same order as in (13). Additional details are seen in Fig. 2. The error in the determination of the scaling exponent is one standard deviation coming from a linear fit on a log-log plot. All fitting errors in this paper are determined in this way.
Let's look at the system magnetization at the critical point. To proceed further it is convenient to define the condensate energy density as which is equal in the SMA to the contribution of the last two terms in Eq. (4) to the condensate energy per volume (the first two terms in Eq. (4) are unimportant in the SMA as the gradient term equals zero while the density-density interaction term ∼ c 0 is a constant). In Fig. 3 we observe that |f T | 2 ∼ τ The upper data presents |fT | 2 /n 2 , where fT is the transverse magnetization, while the lower data is for E/n 2 |c1|, where E (energy density) is given by (14). The pluses come from numerics, while the solid lines present the curves ∼ (τ0/τQ) γ with γ = 2/3 (4/3) for the upper (lower) data. This confirms that at the critical point |fT scalings into Eq. (14) we see that at the critical point in the leading, i.e., τ −2/3 Q , order in the quench rate. It is now interesting to ask whether there is any relation between the condensate magnetization at the point where the system goes out of equilibrium, q = 2 −q, and the condensate magnetization at the critical point, q = 2. The condensate magnetization at q = 2 −q equals (in the slow transition limit whereq ≪ 1) scales with the quench rate in the same way as |f T (q = 2)| 2 . This resembles the adiabatic-impulse simplification of dynamics [3,5], where it is assumed that the evolution of a system undergoing a phase transition is either adiabatic or impulse, and the impulse stage implies no changes in the state of the system. Here this simple picture gives a correct scaling of the system magnetization at the critical point from the value ofq -the location where the system leaves the adiabatic regime. This is a good approximation despite the fact that the condensate does not enter an ideal "impulse" regime. Indeed, the wave-function is not "frozen" and does change from q = 2 −q to q = 2. Remarkably, however, these changes are such that instead of getting magnetization equal to |f T (q = 2 −q)| 2 at the critical point (ideal impulse dynamics) we get const · |f T (q = 2 −q)| 2 , where the constant is τ Q -independent. This character of condensate evolution toward a critical point is universal with respect to quench time. It is depicted in Fig. 4. Now we would like to address how these critical scalings can be experimentally observed. One obvious way is to watch the condensate dynamics at different times as was done in the Berkeley experiment [11]. Another approach involving interesting physics would be to stop a change of q and watch subsequent magnetization dynamics. This approach will be discussed in the next section.

IV. DYNAMICS OF MAGNETIZATION AFTER THE QUENCH
Suppose that one stops driving the system at some q ≥ 2 and lets it evolve freely, which leads to magnetization oscillations. We would like to investigate here whether the information about the system excitation at the critical point can be extracted from the amplitude and/or period of these oscillations. In other words, we want to find out what are the signatures of the nonequilibrium dynamics in the broken-symmetry phase that may be observable in the polar phase. This problem is of both experimental and theoretical interest. The critical point is reached at t/τ0 = 4τQ/τ0 = 40, so q(t ≤ 40τ0) is given by (10) while q(t > 40τ0) = 2. The system excitation due to approaching the critical point is visible from the nonzero magnetization at t/τ0 ≥ 40: if the evolution would be adiabatic, the condensate magnetization at the critical point would be zero and no magnetization oscillations would be present on the plot.

A. Quench stops at the critical point
Suppose the condensate is driven from q = 0 to the critical point with (10). The critical point is reached at t = 4τ Q and then q(t ≥ 4τ Q ) = 2, i.e., free (without any driving) evolution takes place. The evolution of magnetization for such a problem is presented in Fig. 5, where periodic oscillations are easily seen. Our aim is to describe them and show that both their amplitude and their period contain information on how the condensate was excited at the critical point.
In this case one can derive the following equation for the evolution of the magnetization: where E is the system energy density (14) at the critical point (in fact, also at any other time since lack of parameter changes results in energy conservation). This equation can be obtained straightforwardly from Eqs. (6) and (9). It is an interesting equation: in the limit of slow transition, when the amplitude of oscillations becomes very small, the nonlinear term dominates the physics rather than being negligible compared to the ∼ E term promoting harmonic oscillations. Indeed, from Sec. III we see that |c 1 |E/|c 1 | 2 |f T | 2 ∼ (τ 0 /τ Q ) 2/3 ≪ 1 for slow enough quench.
The exact solution reads where cn is the Jacobi cosine [24], and A and θ are given by initial conditions. They both might be found from a fit to experimental data. In particular, |A| 2 is given by the amplitude of free |f T | 2 oscillations, because the Jacobi cosine (16) oscillates periodically between ±1.
The Jacobi cosine was formerly used in the spin-1 condensate context in [12,25]. It is periodic -cn(x, κ) = cn(x + 4K(κ), κ) -where K(κ) is the complete elliptic integral of the first kind. In the time domain this periodicity translates into To simplify this expression we note that from the fact that the Jacobi cosine in (16) is bounded by ±1, f T and A are of the same order and so we expect we get that in the limit of large τ Q (slow transitions) α ∼ (τ Q /τ 0 ) 2/3 ≫ 1 and so κ → 1/ √ 2. This last result allows us to write This finding is quite interesting: the periodicity is amplitude-dependent unlike in the harmonic oscillator case. This simple result predicting the inverse proportionality of the period and the amplitude of oscillations is accurate to about 1% for experimentally relevant quenches as depicted in inset (a) of Fig. 6. The numerical results on the oscillation period ∆t when the quench stops at q = 2 are presented in Fig. 6. From the fit we see that ∆t scales as τ 1/3 Q in accordance with expression (18) supplemented with the above observation that |A| ∼ τ −1/3 Q (proven numerically in the inset (b) of Fig. 6).

B. Quench stops after passing the critical point
According to our numerical simulations there are three stages of such evolution depicted in Fig. 7: (i) the system is driven to the critical point through the entire broken symmetry phase and its magnetization at the critical point scales as |f T | 2 ∼ τ −2/3 Q (only barely noticeable magnetization oscillations are present -see Ref. [23]); (ii) the condensate is driven by a change of q, Eq. (10), in the polar phase and we observe there damped magnetization oscillations; (iii) the free evolution after we stop ramping up q(t) takes place and periodic magnetization oscillations appear. These oscillations are described byf where againĖ = 0 and Eqs. (6) and (9) have been employed in derivation of (19). The exact solution in terms of the Jacobi cosine function exists and has form (16) with α = |Ac 1 | 2 q(q − 2) 2 /2τ 2 0 + 2E|c 1 | . Now, in the limit of fixed q > 2 and large τ Q (slow transition) we have a different behavior of the Jacobi cosine than above. Namely, now α ∼ (τ 0 /τ Q ) 2/3 ≪ 1 which results in κ → 0 and cn(•, κ) → cos(•). In this limit the last two  (18) and (17) divided by (17), i.e., relative departures of the approximate formula for oscillation period from the exact result. The approximate formula (18)  terms in (19) can be neglected and the Jacobi cosine turns into a normal cosine function: harmonic oscillations show up with a repetition period Naturally, the transition between Jacobi cosine oscillations and typical harmonic dynamics is gradual and can be traced quantitatively with the exact solution (16). We can compare (20) to the numerics. For the evolutions where the increase of q stops at q = 3, as is the case in Fig. 7, we get from (20) that ∆t = 2 × 1.814τ 0 , while from numerics we obtain 2 × 1.804τ 0 for τ Q /τ 0 = 30 and 2 × 1.812τ 0 for τ Q /τ 0 = 200 (notice that |f T | 2 oscillates with 1/2 of the period (20) and so our numerical results extracted from |f T | 2 oscillations were multiplied by a factor of two). As we see, the larger τ Q , the better the agreement because the nonlinear corrections become smaller. This is in fact bad news: as (20) is τ Q -independent, we can no longer use repetition period ∆t to investigate the dynamics of the quantum phase transition. Fortunately, however, behavior of the amplitude of free magnetizations oscillations provides an easily visible signature of the non-equilibrium dynamics. As shown in the inset of Fig. 7, the amplitude of |f T | 2 scales as τ −1 Q . This scaling results from two observations. First, the driving in the polar phase starting from q = 2 and proceeding to q 3 damps the amplitude of |f T | 2 by a factor of τ −1/3 Q . We have verified this numerically by simulating different evolutions that begin from a fixed, τ Q independent, state. Second, the scaling of |f T | 2 at the phase boundary is τ −2/3 Q . Combining these two results one easily justifies τ −1 Q scaling from the numerics. Finally, we comment on what happens when the system is continuously driven through the polar phase. Combining Eqs. (6) and (9) one easily arrives aẗ that has to be solved with initial conditions given at instant Two remarks are in order. First, the same initial conditions apply to (15) and (19). Second, equation (21) is no longer self-consistent: an additional equation for ψ 1 (t) dynamics has to be solved simultaneously unless the term ∼ |ψ 1 (t)| 2 is negligible.
As we have observed in Fig. 7 the magnetization oscillations are damped when the system is driven in the polar phase. It turns out, however, that when we continue driving the condensate, the amplitude of magnetization oscillations reaches some nonzero asymptotic value, scaling as τ −1 Q , instead of decreasing to zero. This is supported by  magnetization, |fT | 2 , during an uninterrupted driving through the polar phase with q(t) given by (10). The solid line comes from numerics and provides the upper envelope for oscillations of |fT (t)| 2 . The data is for one run with τQ/τ0 = 30. The dashed line is the asymptotic (large q) value for the solid line. It is obtained from a linear fit (Maxima of |fT | 2 )/n 2 = a + b 1 q for q ∼ 100 data. The fit gives a ≈ 0.0074 plotted as a dashed line. Additionally, we have verified numerically, that for different τQ's the following scaling holds: a ∼ τ −1 Q , which can be explained in the same way as the fitting result from the inset of Fig. 7. 8 and the following analysis. First, we assume that q(t) ≫ 2 such that it is safe to approximate q(t)(q(t) − 2) by q(t) 2 . Second, we neglect the second term in (21) as it is small compared to the first one in the slow transition limit. After that we end up with the following equation that can be solved exactlÿ This is precisely a driven harmonic oscillator with anti-friction sinceq(t)/q(t) is positive in our case. In a classical mechanical picture we may imagine that (22) describes a particle in a harmonic potential under the influence of a force that acts along its velocity. Such a force fights the tendency to squeeze particle motion due to increase of the frequency of the harmonic oscillator. Interestingly, it turns out that the particle perfectly maintains the amplitude of its motion because the exact solution of (22), constrained byφ = 0 from (9), is Additionally, we stress the fact that this solution is valid for any smooth q(t) dependence and not just (10). In particular, in the opposite limit when q(t) decreases in time the "anti-friction" term works as a special friction term that slows the particle down so that it oscillates with constant amplitude instead of spreading out. This phenomenon can be qualitatively explained as follows. As it happens at q(t) ≫ 2 we can safely assume that the magnon modes of the system do not become excited any more during driving. This is so because they have energy gap that becomes large far away from the critical point. These modes are responsible for the scattering of atoms between m = 0 and m = ±1 condensate components [21]. In the absence of that process the populations of m = 0, ±1 sublevels should be constant, which implies that will undergo fixed amplitude oscillations induced by rotation of the condensate phases.

V. DYNAMICS IN A SLIGHTLY INHOMOGENEOUS SYSTEM
In this section we would like to compare the results obtained from a Single Mode Approximation enforcing the translational symmetry onto wave-function of the system to the more experimentally relevant inhomogeneous ("disordered") problem. To this aim we assume that the atom cloud is placed in a box-like trap and impose density and phase fluctuations onto it. These may come from experimental imperfections and quantum fluctuations.
The results presented here come from numerics done in a one-dimensional configuration and so we present an explicitly 1D description below unless stated otherwise. We consider untrapped system: atoms in the box as in the experiment [26] done with spinless bosons. We estimate the parameters for these 1D simulations as in [17] taking into account the experimental setup of [11]. Below L ≈ 400µm is the size of the box keeping N = 2 × 10 6 atoms confined along the z direction.
We introduce disorder by modifying the initial wave function at q(t = 0) = 0 in the following way: where the mode amplitude is proportional to: There δ gives the strength of the disorder, η (m) n , χ ±1 ∈ [0, 2π) are random phases, r (m) n ∈ [0, 1] are random amplitudes, and z (m) n ∈ [0, L] are random positions of gaussian disturbances. All these parameters are generated with uniform probability density. We have chosen σ to be equal to the spin healing length ξ s and assumed that the average spacing between Gaussian perturbations is 2ξ s . The spin healing length in the experiment [11] is 2.4µm, so that κ, the number of Gaussians in (23), is about 80 (L/2ξ s ). The meaning of δ is that it provides a characteristic relative-to-background size of density variations. A typical disordered pattern in (23) is depicted in Fig. 9. Now our goal is to answer what is the range of the δ parameter that allows for observation of qualitatively the same physics as in the homogeneous system described in Secs. III and IV.
To look at the dynamics in the broken-symmetry phase, we have done calculations for different δ ∈ (0, 0.2). For each δ we randomly generate many-times the pattern (23) and evolve the system to the critical point at different quench times τ Q . The patterns are independently generated for every evolution, which shall resemble the experimental situation where the structure of "imperfections" fluctuates from run to run. Then for each pair (δ, τ Q ) we average the mean system magnetization over N r = 30 runs, where f (i) T (z, t) is the condensate transverse magnetization in the i-th run. We look at two quantities: the distanceq(τ Q ) from the critical point where the system starts the non-adiabatic evolution due to approaching the critical point, and the scaling of the condensate magnetization at the critical point.  . The values and errors (one standard deviation) come from a linear fit on a log-log plot to data in the range of τQ = 74ms · · · 1.1s. The Single Mode Approximation (δ = 0) points to λ, σ ≈ 2/3. The definition ofq is the same as in Fig. 2 Both are defined in the same way as in Sec. III except that now we use |f T | 2 instead of |f T | 2 . Our results suggest that the scaling exponents in the disordered system can match the SMA predictions within a few percent accuracy for δ 2%. Quantitatively, as we see in Table I, the scaling exponent λ given byq ∼ τ −λ Q stays reasonably close to the 2/3 value for the homogeneous system for δ as large as 3%. For larger δ, e.g., δ = 5% the numerical data clearly departs from the power law. The scaling exponent σ, defined at the critical point as |f T | 2 ∼ τ −σ Q , deviates noticeably from the SMA prediction already for δ = 2% (Table I). For larger δ, e.g., δ = 3% numerics does not follow the power law anymore.
There are at least two reasons for discrepancies between disordered and homogeneous results. First and most importantly, the presence of disorder perturbs magnetization affecting determination of both exponents from Table  I. Even when these departures from a GS are relatively small compared to the GS magnetization when the system starts time evolution at q = 0, they can be significant when compared to the condensate magnetization close to or at the critical point where scaling exponents are determined. Second, our prediction of the scaling exponents is based on translationally invariant theory, i.e., momentum k = 0 problem, while the introduction of disorder leads to population of k = 0 modes. A more involved description with the gap ∆ (11) being a function of k should be used when the population of k = 0 modes becomes significant. Now we would like to look at the system dynamics in the polar phase. In Sec. IV we have observed that the transverse magnetization undergoes periodic oscillations in a homogeneous problem after stop of driving. Now we see dephasing dynamics spoiling this picture after a time inversely proportional to the disorder magnitude and the quench time τ Q . The latter observation results from the fact that the slower we go at fixed δ, the more significant the disorder is compared to the system magnetization in the polar phase. Fig. 10 presents the evolution with τ Q = 74ms and δ = 1%. The change of q stops at the critical point there. For these parameters the scaling exponents are very close to the Single Mode Approximation result (see Table I). As is illustrated in Fig. 10, however, dephasing takes place after about 7 − 9 oscillations. To observe more magnetization oscillations closely following the SMA result, for twice that time, δ has to be about one order of magnitude smaller. Therefore, the system dynamics is quite sensitive to inhomogeneities when the condensate is left at the critical point. This should result from the coupling between different k modes due to nonlinear terms in the evolution equation. Fig. 11 presents the evolution where driving stops at q = 3 (similarly as in Fig. 7). The parameters and the initial state taken for this simulation are the same as in Fig. 10. This time the undriven dynamics takes place in the regime  The change of q(t) stops at 3 in the polar phase, q(t ≥ √ 24τQ ≈ 0.36s) = 3, and δ = 1%. The insets show enlarged data from the main plot corresponding to the regions depicted with boxes. The period of oscillations is so small that individual oscillations are not resolved in the main plot and the upper inset. The lower inset does not show the SMA result as it is practically indistinguishable from the inhomogeneous calculation. Arrows indicate which inset belongs to what box. The initial state for that evolution is the same as in Fig. 10.
where nonlinear couplings are small. As a result, the dephasing occurs after hundreds of oscillations even though the disorder amplitude, δ = 1%, is the same as in Fig. 10.
To explain dephasing when the evolution stops away from a critical point, we can linearize the coupled Gross-Pitaevskii equations. Below we do it in an explicitly 3D configuration to use the same notation as in Secs. III and IV (the 1D result is qualitatively the same). We assume that where the chemical potential is µ = nc 0 , |δψ m | ≪ 1, and dr(δψ 0 + δψ * 0 ) ≡ 0 to keep drΨ † Ψ = N + O(δψ 2 m ). In leading order f T = √ 2n(δψ * 1 + δψ −1 ), and its dynamics at fixed q is given bÿ where the last term is known from (19). That can be solved by the substitution f T ∼ cos(kr − ω(k)t), giving ω(k) = ǫ 2 k + ǫ k 2(q − 1)n|c 1 | + n 2 |c 1 | 2 q(q − 2), where ǫ k = 2 k 2 /2M . When disorder is present, different k modes will be occupied and they will oscillate at different frequencies ω(k), which will result in the dephasing observed in Fig. 11. The expression for ω(k) shall not be used at q = 2, where linearized theory fails as was shown within the SMA in Sec. IV.

VI. SUMMARY
We have analyzed transitions from the broken-symmetry phase to the polar phase occurring on the finite timescale τ Q . Our focus was on slow transitions, when the condensate goes out of equilibrium near a critical point. We have presented a simple theory predicting that the evolution will cease to be adiabatic before reaching the polar phase at the distance ∼ τ −2/3 Q from the critical point. This result was found to be in excellent agreement with numerics. Applying the basic assumptions of the adiabatic-impulse approach [3,5], which originates from the Kibble-Zurek theory of nonequilibrium dynamics of classical phase transitions [1,2], we have explained why the transverse magnetization of the condensate driven to the critical point scales as τ −2/3 Q as well. Subsequently we have analyzed how the latter result can be experimentally extracted from the observation of magnetization oscillations. Three cases were considered. First, we have studied what happens when the quench stops at the critical point and the system follows free evolution. It was shown that the periodic oscillations appear with the repetition period coupled to the amplitude of transverse magnetization oscillations. The repetition period was found to scale as τ 1/3 Q , where the exponent was directly related to the scaling of the transverse magnetization at the critical point. That nonlinear dynamics was exactly described within the mean-field approach in terms of Jacobi elliptic functions. Second, we have studied what happens when the system is driven into the polar phase, and then undergoes a free evolution. In that case we have shown that its excitation at the critical point can be easily extracted from the scaling of the amplitude of the free magnetization oscillations given by τ −1 Q . Third, we have studied what happens when the condensate is driven without any interruptions not only in the broken-symmetry phase, but also in the polar phase. In this case the amplitude of driven magnetization oscillations reaches a non-zero asymptotic value: the condensate magnetization is described as an unusual harmonic oscillator with anti-friction term that perfectly cancels the squeezing induced by the increase of the harmonic oscillator frequency.
All the above results were obtained within a translationally invariant Single Mode Approximation whose range of applicability was determined by introducing a controlled disorder into an initial wave function. We have found out that the amount of disorder has to be quite small to have the homogeneous system predictions applicable.
Future extensions of this work will include studies of quantum phase transitions dynamics in a harmonically trapped ferromagnetic condensate and investigations of the role of quantum fluctuations on the condensate dynamics.

VII. ACKNOWLEDGMENTS
We gratefully acknowledge the support of the U.S. Department of Energy through the LANL/LDRD Program for this work.