Predicting spinor condensate dynamics from simple principles

We study the spin dynamics of quasi-one-dimensional F=1 condensates both at zero and finite temperatures for arbitrary initial spin configurations. The rich dynamical evolution exhibited by these non-linear systems is explained by surprisingly simple principles: minimization of energy at zero temperature, and maximization of entropy at high temperature. Our analytical results for the homogeneous case are corroborated by numerical simulations for confined condensates in a wide variety of initial conditions. These predictions compare qualitatively well with recent experimental observations and can, therefore, serve as a guidance for on-going experiments.

We study the spin dynamics of quasi-one-dimensional F = 1 condensates both at zero and finite temperatures for arbitrary initial spin configurations. The rich dynamical evolution exhibited by these non-linear systems is explained by surprisingly simple principles: minimization of energy at zero temperature, and maximization of entropy at high temperature. Our analytical results for the homogeneous case are corroborated by numerical simulations for confined condensates in a wide variety of initial conditions. These predictions compare qualitatively well with recent experimental observations and can, therefore, serve as a guidance for on-going experiments. Spinor condensates realized by optically trapped ultracold atoms in a hyperfine Zeeman manifold allow to address a broad scope of problems, which are related to magnetic ordering. These include quantum phase transitions, exotic topological defects and spin domains either within a mean-field regime or within the stronglycorrelated one [1]. Already in the mean-field regime the dynamics of spinor condensates shows some intriguing features, with an apparent randomness regarding the evolution from a given initial state towards a steady state, accompanied by formation of spin domains. Even for the simplest case of an F = 1 spinor condensate, the interplay between spin-spin interactions, non-linear terms and temperature effects rends the analysis of the dynamics rather complex. Previous studies of the spinor dynamics show that, at early stages of the evolution, there is a coherent population transfer between the different hyperfine coupled sublevels [2,3]. Inclusion of temperature (T ) not only smears out the population transfer [2], but also leads to a different distribution of population among the different hyperfine levels.
In this contribution we analyze the dynamics of spinor F = 1 condensates both at zero and finite T from a new perspective. As we shall show, the complex dynamics displayed by these systems can be understood in terms of oscillations in phase space around a steady state. The configuration of this state can be approximately determined by analyzing the trajectories of constant energy in the homogeneous case. To a good approximation, the populations that characterize this state are rather close to those that minimize the energy associated to spin-spin interactions for a given magnetization. In contrast, at finite T , the system is able to exchange energy with the thermal clouds. At large enough temperature, the populations of the steady state can be simply determined by maximizing the entropy of the homogeneous condensate. This leads to a different long-time configuration of the populations than the one attained starting from the same initial conditions at T = 0. Our claims are supported by analytical results for an homogeneous condensate and by numerical investigations for trapped condensates. They provide a good qualitative agreement with recent experimental data on the dynamics of spinor condensates [4,5]. Furthermore, our results can be straightforwardly used to analyze and predict experimental outcomes.
In the mean-field approach, the spinor condensate is described by a vector order parameter Ψ whose components ψ m correspond to the order parameter of each magnetic sublevel |F = 1, m with m = 1, 0, −1. In absence of an external magnetic field and at T = 0, the properties of the condensate are determined by the spin-dependent energy functional [2,6,7]: where repeated indices are assumed to be summed. The trapping potential is taken to be harmonic, axiallysymmetric and spin independent. The coefficients c 0 and c 2 describe the spin-independent and the spin-dependent contributions of the binary elastic collisions of two identical spin-1 bosons. They are expressed in terms of the s-wave scattering lengths a f of the combined symmetric channels of total spin f = 0, 2 as c 0 = 4π 2 (a 0 + 2a 2 )/(3M ) and c 2 = 4π 2 (a 2 − a 0 )/(3M ), where M is the atomic mass. F are the spin-1 matrices [6,7]. In this Letter, we focus on 87 Rb spinor condensates [4,5,8,9], which have ferromagnetic character (c 2 < 0).
According to i ∂ψ m /∂t = δE/δψ * m , a set of coupled dynamical equations for the spin components are ob-tained: with H s = − 2 /(2M ) ∇ 2 + V ext + c 0 n being the spinindependent part of the Hamiltonian. The density of the m-th component is given by |ψ m (r)| 2 , while n(r) = m |ψ m (r)| 2 is the total density normalized to the total number of atoms N . The population of each hyperfine state is N m = dr|ψ m (r)| 2 . Defining the relative populations λ m = N m /N , the magnetization of the system is M = λ 1 − λ −1 . The total number of atoms and the magnetization are both conserved quantities [2]. Then, the conditions λ 1 + λ 0 + λ −1 = 1 and M = λ 1 − λ −1 , together with the positiveness of λ m , restrict the allowed configurations of the system in the (λ 0 , M) phase space to the region below the dotted line in Fig. 1.
To understand the physics governing the spin dynamics, it is instructive to analyze first the homogeneous system at T = 0. Translational invariance leads to the following ansatz: Notice that, since the system is homogeneous, the value of n is not relevant to determine the relative populations that minimize the energy. Thus, the only relevant contributions to the energy are those associated to c 2 : For a given value of M, the ground state configuration is obtained by minimizing Eq. (3) with respect to λ 0 and θ. For ferromagnetic systems, this yields θ = 0 and which defines the solid line drawn in Fig. 1. For instance, in the symmetric case (M = 0) the minimization leads to a ground state with λ 0 = 0.5 and θ = 0.
We study now the dynamical evolution of the homogeneous system by solving numerically the set of Eqs. (2), after removal of the terms − 2 /(2M ) ∇ 2 + V ext . We consider three independent phases θ m in the time evolution, although the energy depends only on the relative phase θ [Eq. (3)]. As expected, we find that starting from an equilibrium configuration (e.g. λ 0 = 0.5 and θ = 0 for M = 0) there is no evolution at all. However, if the initial configuration does not correspond to the equilibrium state, then the system evolves exchanging populations between spin components, but conserving energy and magnetization.
Therefore, in the homogeneous system the dynamical evolution of the spin populations is characterized by trajectories of constant energy and magnetization in the λ 0 − θ plane, which are defined by the starting values of M, λ 0 and θ. In most cases, i.e. if the initial conditions are not very far from the equilibrium configuration, defined by λ eq 0 and θ = 0, the trajectories in the λ 0 − θ plane are closed orbits (see also [11]). Then, the evolution of the populations -in particular of λ 0 -shows up as oscillations of a well defined frequency, around the geometrical center of the orbit. An analysis of the orbits defined by Eq. (3), shows that these centers lie always between λ eq 0 and λ min 0 = (1 − |M|)/(2 − |M|), represented by the dashed curve in Fig. 1. With this information in mind we move now to the confined systems.
To this aim, we solve the equations of motion for 20000 atoms of 87 Rb, in a highly elongated trap ω ⊥ = 2π × 891 Hz and ω z = 2π × 21 Hz, i.e. ω ⊥ ≫ ω z , so that Eqs. (2) become quasi-one-dimensional. The corresponding scattering lengths are a 0 = 101.8a B and a 2 = 100.4a B , and the coupling constants c 0 and c 2 are accordingly rescaled by a factor 1/(2πa 2 ⊥ ), being a ⊥ = /mω ⊥ the transverse oscillator length [10]. The validity of this approximation for the case of a spin-1 atomic condensate has been proved for a very wide range of conditions in Ref. [11]. We take as initial wave functions ψ m (z, t = 0) the wave function of the ground state of the scalar condensate, obtained from H s with coupling constant c 0 /(2πa 2 ⊥ ), and normalized to their corresponding initial populations N m . This corresponds to a Single Mode Approximation (SMA) ansatz for the initial conditions. SMA provides a good description of spin systems in certain conditions [12,13]. However, in solving the time evolution equations we keep the full dependence of ψ m on position and therefore we go beyond the SMA. Actually, this is a necessary condition to study the formation of spin domains.
The numerical solutions of the dynamical equations show oscillations of the integrated spin populations around a steady state, which coincides with that in the homogeneous case [Eq. (4)]. Indeed, if we start with a configuration lying on the curve defined by λ eq 0 with θ = 0, the system does not evolve at all. However, any initial configuration with (λ 0 ,M) on the curve, but θ = 0, yields dynamical evolution of the spin populations around the corresponding steady configuration of the homogeneous case.
Our results are summarized in Fig. 1. The dots correspond to the time average value of λ 0 (t) obtained by solving the equations of motion for different initial outof-equilibrium values of M, λ 0 , and θ = 0. All cases studied would correspond in the homogeneous system to closed orbits in the λ 0 − θ plane. The magnetization is constant throughout the evolution, while λ 0 (t) presents slightly damped oscillations around the same λ 0 as in the homogeneous case, which is restricted to be in the narrow region limited by the solid and the dashed curves described above. Thus, the curve λ eq 0 provides a reliable estimation of the configuration of the steady state. Notice that in all reported cases, the dots lie very close to the curve λ eq 0 as we have taken, in the selection of the starting conditions for M = 0, configurations and phases not far from the equilibrium conditions.
As an example of the dynamical evolution of the system at T = 0, we display in Fig. 2(a) the spin dynamics for an initial symmetric spin configuration corresponding to λ 1 = λ −1 = 5% and θ = 0. The oscillations between m = 0 and m = ±1 populations are not regular and present a dynamical instability around t dom , when the condensate starts the multidomain formation process [2,14,15,16]. An estimate of the time scale for the appearance of the instability is t dom ∼ 2π /(c 2 n 0 ) [15,16], where n 0 is the density at the center of the trap. In our simulation, n 0 ∼ 2 × 10 14 cm −3 , leading to t dom ∼ 140 ms, which is in agreement with the dynamics of the system. Around this time, the condensate starts the formation of small dynamical spatial spin structures, also observed experimentally [4,9], while the relative populations oscillate around the equilibrium state (λ eq 0 = 50%, λ eq ±1 = 25%). Due to this inhomogeneity of the phases and density profiles [2], it is not possible to associate the dynamics to a well-defined orbit in the λ 0 − θ plane.
It is interesting to note the formation of spin domains despite the total density profile remaining almost unchanged [2,15]. Hence, spin domains are linked to the kinetic and spin-spin energy terms, since all the other contributions to the energy functional depend only on the total density. The evolution of the kinetic energy and E c2 is shown in Fig. 2(c)   first 500 ms. Clearly the appearence of domains is linked to a rapid variation of kinetic and E c2 energies. Moreover, since the system has no dissipation at T = 0, the total energy has to be conserved, and therefore, the energy difference between the starting and the equilibrium configurations provides an upper bound to the energy associated with the formation of spin domains. Actually this energy is always rather small, as can be seen in the figure. In the case of Fig. 2, the excess energy amounts to 0.28 ω z which represents a 0.3% of the total energy.
We incorporate now finite-temperature effects by means of a Bogoliubov-de Gennes description of the thermal clouds. Further we assume that each spin component has its own cloud [17], and solve the corresponding timedependent equations under these conditions [2]. With our choice of parameters, the condensation in an elongated trap occurs at k B T 3D [18]. In Fig. 2(b) we show for T = 0.2T 3D c the time evolution of the relative populations, starting from the same initial parameters as in the T = 0 case. It is worth noticing that even if T ω ⊥ , the quasione-dimensional description of the thermal effects is still valid [19].
A simple comparison between Fig. 2(a) and 2(b) shows that at finite T the oscillations in the populations are clearly damped and that the steady state configuration that is reached is different from the one at T = 0. In particular, in Fig. 2(b) we observe that starting from the same initial configuration (5%, 90%, 5%), the long-time population distribution changes from (25%, 50%, 25%) at T = 0 to a steady configuration with all m-states almost equally populated (equipartition). Moreover, even though the local magnetization is no longer conserved as it was at T = 0 [2], the total magnetization is still a conserved quantity along all the time evolution. In general, temperature effects reflect in a faster damping and in an earlier domain formation in comparison to T = 0. A way to understand this result is to assume that the steady configuration is attained by maximizing the entropy, subjected to the conservation of total population and magnetization: where we have introduced the Lagrange multipliers Λ and µ associated with the magnetization and the normalization, respectively. Imposing that ∂S/∂λ m = 0, together with the normalization and magnetization constraints, leads to the following condition for the steady state configurations: which corresponds to the dot-dashed curve plotted in Fig. 1. Notice that the curve thus obtained is independent of temperature. However, the hypothesis underlying its derivation is appropriate for a high temperature regime. Therefore, one expects that this curve together with the one for the T = 0 case [Eq. (4)] limit the zone of the configurations around which the system will oscillate depending on temperature. In fact, the results of our numerical simulations for different temperatures and initial configurations, shown in Fig. 1, are spread between the two curves. We present results for T /T 3D c = 0.05 (rhombi), 0.1 (squares), 0.2 (triangles) and 0.4 (inverted triangles). For a given initial configuration, when the temperature increases, the steady state gets closer to the curve of maximum entropy. This is illustrated for M = 0.2 and 0.4. In both cases, the initial configuration is indicated by open squares. Notice that for M = 0.2, the initial configuration is on the equilibrium curve (T = 0) and its evolution is a direct consequence of the temperature. Also one observes a larger dependence of the final steady state configuration on the initial conditions than at T = 0. In the figure, this is illustrated for M = 0, and 0.1 by the presence of more than one symbol belonging to the same temperature that correspond to different initial conditions. Summarizing, we have shown that at zero temperature and in absence of external magnetic fields the spin populations of F = 1 condensates oscillate around a steady state which is well determined by analyzing the simpler homogeneous case. This steady state is close to the one described by minimizing the E c2 energy at constant magnetization for the homogeneous case [Eq. (4)]. In confined systems, spin domain formation is associated to the excess of spin-spin coupling energy, which converts into kinetic energy back and forward dynamically.