Dynamics of F=1 87Rb condensates at finite temperatures

We investigate the dynamics of a F=1 spinor Bose-Einstein condensate of 87Rb atoms confined in a quasi-one-dimensional trap both at zero and at finite temperature. At zero temperature, we observe coherent oscillations between populations of the various spin components and the formation of multiple domains in the condensate. We study also finite temperature effects in the spin dynamics taking into account the phase fluctuations in the Bogoliubov-de Gennes framework. At finite T, despite complex multidomain formation in the condensate, population equipartition occurs. The length scale of these spin domains seems to be determined intrinsically by nonlinear interactions.


INTRODUCTION
The seminal theory papers of T.-L. Ho [1] and T. Ohmi and K. Machida [2], as well as the experiments performed by the MIT group on optically trapped sodium Bose-Einstein condensates (BEC) [3] have stimulated the development of a new interesting area of research in the field of multi-component quantum gases: BEC with spin degrees of freedom, i.e. the so called spinor condensates. The bosonic quantum field operator in such systems is no longer a scalar, but a vector, and depending on various parameters the system can be found in a coherent superposition or in an incoherent mixture of condensates with different spin components, that can exchange population depending on the intrinsic nonlinear coupling. The couplings between the different spin components, similarly as in the standard scalar BEC can be described by zero range potentials that however, acquire a matrix form and describe couplings that assure total spin conservation in elastic binary collisions [4].
It is worth mentioning that spinor condensates are closely related to the effective spin-1/2 systems realized by radio-frequency coupling of the two hyperfine states in 87 Rb [5]; in such systems spin-waves above the critical temperature [6] and decoherence effects were observed [7]. From the theoretical point of view, collective modes in two-component BEC's have been studied in [8], while an analysis on the formation of spin domains can be found e.g. in [9].
Spinor condensates are suitable systems to study various quantum phenomena, that do not occur in single component BEC's. Equilibrium states of spinor condensates in an optical trap can exhibit magnetic ordering of various kinds. For instance, for sodium (F = 1), the free-dom of spin orientation leads to the formation of spin domains in an external magnetic field, which can be either miscible or immiscible with one another [3]. Excitations in such systems may manifest the different spin character; stability of topological excitations and textures, such as ordinary vortices, coreless vortices (Skyrmions) [1,10], and t'Hooft-Polyakov monopoles [11], depends in a complex manner on the parameters of the system. The complexity of these systems becomes even more transparent in the limit of strongly correlated systems: ultracold Bose spin gases in optical lattices exhibit fascinating properties, including new types of quantum phases, such as a polar condensate phase, a condensate of singlet pairs, a crystal spin nematic phase, and a spin singlet crystal phase [12].
The ground state, its magnetic properties and the low temperature thermodynamics of spinor condensates have been studied in several experiments: in 23 Na [13] (F = 1) which has an anti-ferromagnetic ground state, in 87 Rb in the F = 1 spin state [14,15,16,17] which is ferromagnetic, as well as in the F = 2 spin state, which presents a rich ground state behavior [14,17]. Recent experiments involving non-destructive imaging of magnetization of a spin F = 1 Bose gas of 87 Rb atoms [18] open a new route to study magnetism of spinor condensates. Also, the recent experiment with chromium [19] condensates opens the possibility to study even higher spin states (F = 3) with long range (dipolar) interactions and, presumably a much more complex phase diagram.
In the last two years the focus of investigation has shifted towards two other very important aspects of the spinor BEC physics: the equilibrium properties at finite T , and dynamics of spinor BEC's. The thermodynamic properties of an ultra-cold spin Bose gas have been experimentally investigated in Ref. [20]. Recently, finite tem-perature effects to describe the properties of the equilibrium density distribution have been considered within the Hartree-Fock-Popov theory [21]. Spin dynamics has been experimentally studied in 87 Rb condensates for F = 1 in [15], and more exhaustively for F = 2 in [14,15,17]. Coherent collisional spin dynamics in 87 Rb for F = 2 has been recently observed also in optical lattices [22] -these experiments pave way towards the efficient creation of entangled atom pairs in optical lattices.
Here we focus on the dynamics of F = 1 elongated spinor condensates at zero and finite temperatures. One of our major motivation to study the thermal effects on spin dynamics is to learn something about decoherence in multi-component systems. In the case F = 1 the internal coupling of the spin components depends only on 2 coefficients, while for F = 2 it depends on 3. Due to this simplicity, the case F = 1 allows for a better understanding of the interplay between nonlinear interactions, spin couplings and temperature effects. We consider here the full coupled dynamical equations of the spin components, obtained within a mean-field framework, without any further approximation such as the Single Mode Approximation (SMA) [23,24], or variational ansatz [25], that would mask some aspects of the complex dynamical evolution. We restrict our analysis to the quasi-1D case, when the condensates are kinematically frozen in the transverse directions. We investigate then the influence of thermal effects on the spin dynamics. To this aim we use the Bogoliubov-de Gennes description of phonon modes and treat them as classical random fields, similarly as in Ref. [26]. For quasi-1D condensates at finite T the main contribution to phonon fluctuations comes from the fluctuations of the condensate phase [27], and this contribution is fully accounted in our simulations.
The paper is organized as follows. First, in section II we introduce the model for T = 0. We describe some of the details of the numerical method in section III. In section IV we describe our results for T = 0. Section V we present a discussion of finite temperature effects. We conclude in section VI.

DESCRIPTION OF THE SYSTEM
The many-body Hamiltonian describing a F=1 spinor condensate in absence of an external magnetic field is given by [1] is the field operator that annihilates (creates) an atom in the m-th hyperfine state |F = 1, m at point r (m = 1, 0, −1). The trapping potential V ext (r) is assumed harmonic and spin independent. The terms with coefficients c 0 and c 2 describe binary elastic collisions of spin-1 atoms in the combined symmetric channel of total spin 0 and 2, and are expressed in terms of the s-wave scattering lengths a 0 and a 2 : 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 [1,21]. The system is ferromagnetic if c 2 < 0 ( 87 Rb), and anti-ferromagnetic if c 2 > 0 ( 23 Na).
In the mean field approach, a condensate order parameter is introduced for each magnetic sublevel ψ m (r) = Ψ m (r) , and by neglecting quantum fluctuations it yields the following energy functional According to i ∂ψ m /∂t = δE/δψ * m , the coupled dynamical equations for the spin components are obtained [23,28] n m (r) = |ψ m | 2 is the density of the m-th component and n = |ψ 1 | 2 + |ψ 0 | 2 + |ψ −1 | 2 is the total density normalized to the total number of atoms N . The population of the hyperfine state |1, m is N m = dr|ψ m | 2 such that , and the effective potentials that will determine the spatial dynamics of each component Analogously to a spin-polarized condensate, the multicomponent Gross-Pitaevskii equations (3) can be rewritten in the form of continuity equations [29]: where j m = (ψ * m ∇ψ m − ψ m ∇ψ * m )/2iM is the current, and δṅ m (r, t) = −(2c 2 / )Im(T m ψ m ) is the rate of transfer of populations between spin components. Since the total number of atoms and magnetization are conserved, and the dynamical equations for m = ±1 are symmetric, it is verified that . In our calculation, we assume an axially-symmetric harmonic confinement V ext = M (ω 2 ⊥ r 2 ⊥ + ω 2 z z 2 )/2. In the limit of highly elongated traps (ω ⊥ ≫ ω z ) the tight confinement ensures that no excited states are available in the transverse direction and thus the dynamics takes place along the axial direction. Factorizing ψ m (r) into a longitudinal and a transverse function, and approximating the transverse part as the two-dimensional ground state of the transverse oscillator [23,25], the equations of motion (3) become one-dimensional for the longitudinal wave functions ψ m (z), and the coupling constants c 0 and c 2 are accordingly rescaled by a factor 1/(2πa 2 ⊥ ), with a ⊥ = /mω ⊥ the transverse oscillator length [30].

NUMERICAL METHOD
The dynamics of the spin components is obtained by numerically integrating the coupled nonlinear differential equations (3). For the time evolution our numerical procedure combines the split operator method with the Fast Fourier Transform to treat the kinetic terms and a fourth-order Runge-Kutta method for the remaining terms of the dynamical equations. We have compared our combined numerical method with an evolution using a pure fourth-order Runge-Kutta as the one used in Ref. [23], and we have found that our method allows larger time steps making the computation notably faster.
In this paper, we consider N = 20000 atoms of spin-1 87 Rb trapped in a quasi-1D harmonic trap with ω ⊥ = 2π × 891Hz and ω z = 2π × 21 Hz. The coupling constants are a 0 = 101.8a B , and a 2 = 100.4a B [31], with a B the Bohr radius. Since c 2 < 0 the atomic interactions for 87 Rb atoms are ferromagnetic. The initial wave functions ψ m (z, t = 0) have the same spatial profile as the ground state of the scalar condensate with coupling constant c 0 /(2πa 2 ⊥ ). Each initial wave function component is however normalized to the corresponding initial population in spin-m, N m . Since ψ m (z, t) are complex functions, ψ m (z, 0) = |ψ m (z, 0)| exp (iθ m ), the initial phases have also to be fixed. However, for initial symmetric configurations (N 1 = N −1 ) the dynamical evolution depends only on one initial relative phase θ = 2θ 0 − θ 1 − θ −1 [23], and thus only one initial phase has to be fixed. Experimentally the initial phases are not well determined and the results are averaged over several repeated measurements [15]. Therefore we will average the dynamics over different initial relative phases randomly distributed over (0, 2π) in order to obtain results easy to compare with experiments.

SPINOR DYNAMICS AT T = 0
We consider that initially a quasi-pure condensate in the m = 0 spin component is populated. Spin component mixing requires, at least, a small seed of atoms populating the other components but keeping the total magnetization equal to zero. In Figure 1 we plot the population of each spin component as a function of time at zero temperature, for the initial populations (N 1 /N, N 0 /N, N −1 /N ) = (0.5%, 99%, 0.5%). In absence of an external magnetic field gradient the dynamical evolutions of the spin m = ±1 components are symmetric. Dashed lines correspond to the dynamical evolution with a given initial relative phase θ = 0, and solid lines are the average over 20 random initial relative phases. The oscillations of the populations with a given initial phase, are smeared out by the average over different runs. It is interesting to point out the damping of the oscillations obtained in our numerical calculations which is a finite size effect related to collapse phenomena characteristic of discrete anharmonic spectra [32], and to dephasing of Josephson oscillations [33]. During the time evolution the magnetization is conserved, as it has been experimentally observed [15]. The populations oscillate around the ground state configuration of the system with M = 0 that in absence of applied external magnetic field is (25%, 50%, 25%) [28]. Our numerical results are in qualitative agreement with the experimental measurements of Ref. [15] obtained in a strongly anisotropic diskshaped spin-1 condensate, where the relaxation to the steady state is also not monotonic but damped.
In Fig.2 we plot the density profiles of the spin components for the configuration θ = 0 (see Fig.1). At the initial stages of the evolution, the population of the spin-0 component decreases due to the spin exchange interaction, and thus the ±1 spin components start to be populated by the same amount, keeping the symmetry of the initial state. The total magnetization is conserved along the time evolution. Initially (t < 100 ms), the conversion of atoms from 0 to ± 1 states mainly occurs at the central part of the condensate, where the density is higher and thus the coupling between different spin components is more effective, see Eq. (4).
Then, the ± 1 spin components swing back to the 0 component, and vice versa. The oscillations between the populations of the m = 0 and m = ±1 states are not regular and present a dynamical instability around t ∼ 100 ms, when the large amplitude oscillations become small amplitude oscillations [34,35,36]. At this moment the condensate starts the multi-domain formation process into small dynamical spin domains. A simple estimation of the time scale for the appearance of the instability t dom ∼ 2π /(c 2 n) has been provided by studying the normal excitation modes of the system [35,36]. In our case, taking n ∼ 2 × 10 14 cm −3 at the center of the trap, t dom ∼ 140 ms, which is in agreement with our results. Notice that, in the present study, the analysis is performed directly from the spatial and temporal evolution of the population of the spin components.
The transfer of population remains coherent along all the time evolution, and it does not become chaotic even for large times. The number of small spin domains does not grow indefinitely, but it is limited by a characteristic size of the spin domains (l dom ) that depends on the in-ternal coupling between different spin components, and more weakly on temperature. The formation of multiple spin domains, also found recently in Ref. [35,36], is a result of the ferromagnetic character of the spin-1 87 Rb condensate: c 2 < 0 favors the spatial separation of ±1 atoms from 0-atoms, as obtained in our numerical results. Since in our simulations the initial condition for m = ±1 are identical, the density profiles of m = ±1 components are equal along the time evolution. Apparently, the value of l dom ∼ 5 µm, being much smaller than the harmonic oscillator length, does not depend on the external trapping potential, but it is an intrinsic characteristic of the spin coupling strength [36].
The dynamical evolution between spin components can be easily interpreted in terms of the effective potential V ef f m (z, t) felt by each spin component ψ m (z, t), and from the continuity equations. For the same initial configuration as in Fig.2, we plot in Fig. 3 the effective potentials V ef f 0 and V ef f 1 = V ef f −1 (top panels), and the local transfer of population δṅ 0 (z, t) and δṅ 1 (z, t) = δṅ −1 (z, t) (bottom panels) for t = 0, 40 and 80 ms. For a fixed time t, δṅ m (z, t) represents the local exchange of atoms between spin components. If δṅ m (z, t) < 0 at that position the population of the spin-m component is decreasing and thus the atoms convert to the other spin components. For example at the first stages of the time evolution (t ≤ 60 ms) δṅ 1 (z, t) = δṅ −1 (z, t) > 0, and consequently δṅ 0 (z, t) < 0: the spin-exchange interaction is favoring the conversion from 0 to ±1 states as it is shown in Fig.1. At t = 0 this is also true, but it cannot be appreciated from the scale of the figure. Moreover, since the minimum and maximum of δṅ 0 (z, t) and δṅ ±1 (z, t), respectively are at the center of the condensate, the spinexchange mainly occurs at the central region, as we have already commented in Fig.2. At t = 80 ms δṅ ±1 is positive at the boundaries of the condensate, and negative at the central region, whereas δṅ 0 has the opposite behavior. Therefore, the population with m = ±1 is decreasing at the center and increases at the boundaries, and the ±1 atoms convert to 0 state mainly at the center, where the m = 0-condensate develops a new central peak.
We have also performed simulations starting from other initial conditions. In particular, we observe that multi-domain formation is inhibited if the starting configuration coincides with the ground state composition (25%, 50%, 25%). Moreover, we find a good agreement with the very recent experimental results of Chang et al. [37], starting from (0, 3/4, 1/4) and converging to (1/5, 2/5, 2/5).

SPINOR DYNAMICS AT FINITE T
The spinor dynamics of these multicomponent quantum gases is rather complex due to the internal coupling between different spin components. We consider now thermal effects in the spinor dynamics. At low temperature, thermal excitations can be described within the Bogoliubov-de Gennes theory. Recently, finite temperature effects in the equilibrium density distribution of the condensed and non-condensed components of spin-1 trapped atoms has been investigated within Hartree-Fock-Popov theory [21]. We will investigate thermal effects in the spinor dynamics using a Bogoliubov-de Gennes description of the thermal cloud. For a highly elongated trap, the condensate is quasi-1D, and total density fluctuations are strongly suppressed at small temperatures, whereas phase fluctuations are relevant and in the Thomas-Fermi regime can be described analytically in terms of Legendre polynomials P j (z) [27,38] Analogously, we assume that spin fluctuations can be disregarded in a first approximation.
In Refs. [7,20] it has been shown that at finite temperature each spin component has its own thermal cloud. Thus, we assume that initially the condensate corresponding to each spin-m component can be described by an order parameter ψ m (z) = n m (z) exp(iφ m ), which has a random fluctuating phase [39] where n 0 (z) is the equilibrium total density profile of the initial condensate, µ and R TF are its chemical potential and Thomas-Fermi radius, ǫ j = ω z j(j + 1)/2 is the spectrum of low lying axial excitations [40], and a m j (a m * j ) are complex amplitudes that replace the quasiparticle annihilation (creation) operators in the meanfield approach. In the numerical calculation in order to reproduce the quantum statistical properties of the phase [26], a m j and a m * j are sampled as random variables with a zero mean value and |a m j | 2 = N m j , where N m j = [exp(ǫ j /k B T ) − 1] −1 is the occupation number for the quasi-particle mode j [41,42].
In Fig.4 we plot the dynamical evolution of the spin populations at T = 0.2T c for the same initial populations as in Fig.1. The critical temperature of Bose-Einstein condensation, T c = ω z N/ ln(2N ), corresponds here to a single component Bose gas in a harmonic trap of frequency ω z . Solid lines correspond to the numerical results averaged over 20 random initial relative phases, and dashed lines to a single run with initial phase θ = 0. The interaction of the condensate atoms with the thermal clouds smears out the oscillations present at T = 0 ( Fig.1) and leads to an asymptotic configuration with all m components equally populated (equipartition).
Thermal effects in the density profiles of the spin components are shown in Fig.5. The occurrence of phase fluctuations due to thermal excitations at t = 0 transforms to modulations of the density during the time evo-lution as in a spin-polarized single component elongated condensate [26]. In a spinor condensate this leads to multi-domain formation at much earlier times than at T = 0 (see Fig.2), and the density profiles are no longer symmetric at finite T . Moreover, the existence of different thermal clouds for each spin component [20] breaks also the symmetry of the m = ±1 spin components of the initial state of system. Therefore the m = ±1 density profiles are different during the time evolution, and the three components separate in different spin domains. The local magnetization is not longer conserved as at T = 0 but, as expected, the total magnetization is still a conserved quantity along all the time evolution. Similar results are obtained at lower temperatures T /T c = 0.01 and 0.001.

CONCLUSIONS
In this paper, we have studied the spinor dynamics of a spin-1 87 Rb condensate in a highly elongated trap. We have solved the full three coupled dynamical equations for the spin components within Gross-Pitaevskii framework without any further approximations. This is in fact necessary, since approximated approaches frequently mask some of the aspects of the dynamics. We have also considered here finite temperature effects, using the approach of Ref. [41]. We have found that the spinor dynamics towards the steady state is not monotonic but rather slowly damped, involving a coherent transfer of population between different spin components. At finite temperature the coherent oscillations of populations are almost smeared out. The internal coupling of the spin components leads to the formation of multiple spin domains of a small, but finite characteristic length l dom , which does not decrease with time, and seems thus to be determined intrinsically by the nonlinear interactions. This scale is evidently larger than the condensate healing length, l heal = 2π / √ 2M c 0 n, which for scattering lengths a of order of 5 nm is of order of 10-100 nm. In fact, l dom is of the order of the spin healing length ξ s = 2π / 2M |c 2 |n. The presence of different thermal clouds for each spin component breaks the symmetry of the m = ±1 components, and therefore separates them in different spin domains. For a condensate with initially zero magnetization, the spin populations oscillate around the ground state configuration (25%, 50%, 25%) at T = 0, whereas at finite temperature the interaction of the condensate atoms with their corresponding thermal clouds leads to equipartition in populations, i.e. (1/3, 1/3, 1/3).
Our results shed, in our opinion, also some light on the question of decoherence. Our simulations, and in particular the finding of multi-domain formation suggest that decoherence undergoes enhancement with the number of components in the system. Of course, there are many open questions connected to this, e.g. does the multi-domain formation go along with a loss of phase relations, and give rise to some enhanced (generalized) phase fluctuations? These questions go beyond the present study, and will be discussed in detail in a future publication, in which we will compare the F = 1 and F = 2 case. and n0(z, t) (green). The initial configuration corresponds to (0.5%, 99%, 0.5%), and θ = 0 (dashed line in Fig.1).