Spin fragmentation of Bose-Einstein condensates with antiferromagnetic interactions

We study spin fragmentation of an antiferromagnetic spin 1 condensate in the presence of a quadratic Zeeman (QZ) effect breaking spin rotational symmetry. We describe how the QZ effect turns a fragmented spin state, with large fluctuations of the Zeemans populations, into a regular polar condensate, where atoms all condense in the $m=0$ state along the field direction. We calculate the average value and variance of the Zeeman state $m=0$ to illustrate clearly the crossover from a fragmented to an unfragmented state. The typical width of this crossover is $q \sim k_B T/N$, where $q$ is the QZ energy, $T$ the spin temperature and $N$ the atom number. This shows that spin fluctuations are a mesoscopic effect that will not survive in the thermodynamic limit $N\rightarrow \infty$, but are observable for sufficiently small atom number.


Introduction
The natural behavior of bosons at low enough temperatures is to form a Bose-Einstein condensate, i.e. a many-body state where one single-particle state becomes macroscopically occupied [1]. There are, however, situations where bosons can condense simultaneously in several single-particle states, forming a so-called fragmented condensate. Several examples are known, where fragmentation occurs due to orbital (Bose gases in optical lattices or in fast rotation) or to internal degeneracies (pseudospin 1/2 or spin 1 Bose gases). These examples have been reviewed in [2,3].
The spin 1 Bose gas, first studied by Nozières and Saint James [4], is a striking example where fragmentation occurs due to rotational symmetry in spin space. For antiferromagnetic interactions of the form V 12 = g s s 1 · s 2 between two atoms with spins s 1 and s 2 (g s > 0), the many-body ground state is expected to be a spin singlet state [2,5]. In such a state the three Zeeman sublevels are occupied, leading to three macroscopic eigenvalues of the single-particle density matrix (instead of just one for a regular condensate). As pointed out in [5,6,2], the signature of fragmentation is then the occurrence of anomalously large fluctuations of the populations N m in the Zeeman states m = 0, ±1 (see also [7], where a similar behavior is predicted in a pseudo spin 1/2 system). In the singlet state for instance, the expectation value and variance of N 0 are N 0 = N/3, and ∆N 2 0 ≈ 4N 2 /45, respectively (N is the total number of particles). Such super-Poissonian fluctuations (∆N 2 0 ∝ N 0 2 ) deviate strongly from the value expected for a single condensate or any ensemble without correlations where ∆N 2 0 ∝ N 0 §. It was pointed out by Ho and Yip [6] that such state was likely not realized in typical experiments, due to its fragility towards any perturbation breaking spin rotational symmetry (see also [8,9,10,11,12,13]). In the thermodynamic limit N → ∞, an arbitrary small symmetry-breaking perturbation is enough to favor a regular condensed state, where almost all atoms occupy the same (spinor) condensate wave function and ∆N 0 N . In this article, we give a detailed analysis of the phenomenon of spin fragmentation for spin 1 bosons. Our analysis assumes the conservation of the total magnetization m z . The fact that magnetization is an (almost) conserved quantity follows from the rotational invariance of the microscopic spin exchange interaction, and from the isolation of atomic quantum gases from their environment. A key consequence is that in an external magnetic field B, the linear Zeeman effect only acts as an energy offset and does not play a role in determining the equilibrium state. The dominant effect of an applied magnetic field is a second-order (or quadratic) Zeeman energy, of the form q(m 2 − 1) for a single atom in the Zeeman state with magnetic quantum number m . § Note that the problem we discuss here is unrelated to the anomalous fluctuations of the total condensate number found for ideal gases in the grand canonical ensemble [1]. In this work, we assume implicitly the canonical ensemble, and study the fluctuations of the populations of individual Zeeman states discarding quantum and thermal depletion of the condensate.
This second-order shift originates from the hyperfine coupling between electronic and nuclear spins, and corresponds to the second order term in an expansion of the well-known Breit-Rabi formula for The QZ energy breaks the spin rotational symmetry, and favors a condensed state with m = 0 along the field direction. In [9,10,11,12], the evolution of the ground state with the QZ energy q was studied theoretically. Since experiments are likely to operate far from the ground state, it is important to understand quantitatively how the system behaves at finite temperatures. This is the main topic we address in this paper.
Our focus in this article will be to calculate the first two moments (average value and variance) of N 0 . These moments illustrate clearly the evolution of the system from fragmented to unfragmented and thus constitute the main experimental signature of fragmentation. The main findings are summarized in Figure 1, where we plot the standard deviation of n 0 = N 0 /N in a q − T plane. Large fluctuations and depletion of the m = 0 state are observed for small q. We can distinguish three different regimes. For low q U s /N 2 and low temperatures k B T U s /N (U s ∝ g s is the spin interaction energy per atom), the system is close to the ground state in a regime we call "quantum spin fragmented" [5,6,2,11]. We also observe a thermal regime for k B T N q, U s /N dominated by thermally populated excited states. We call this second regime "thermal spin fragmented". Finally, for q large enough and temperature low enough, the bosons condense into the single-particle state m = 0, forming a so-called "polar" condensate [15,16]. In this limit, N 0 ≈ N and ∆N 0 N . We indicate this third regime as "BECm = 0" in Figure 1.
The evolution from the fragmented, singlet condensate to an unfragmented condensate with increasing QZ energy q is similar to a well-known example in the literature on quantum magnetism, the Lieb-Matthis model of lattice Heisenberg antiferromagnets [17]. This model describes collective spin fluctuations of an Heisenberg antiferromagnet on a bipartite lattice. It constitutes a popular toy model for demonstrating the appearance of broken symmetry ground states in condensed matter [18,19,20,21,22]. The ground state of such system (in principle also a spin singlet) was found theoretically to evolve to a Néel state in the thermodynamic limit in the presence of an arbitrarily small staggered magnetic field (whose sign alternates from one site to the next). The underlying theory is close to the one presented here. An essential difference is that the present model of antiferromagnetic spin 1 BECs is expected to accurately describe actual experimental systems [23,24]. In the antiferromagnet case, the staggered magnetic field is a theoretical object that cannot be produced in the laboratory for real solids. In contrast, the QZ energy is easily controllable in spin 1 BEC experiments. Another important difference is that experiments with ultracold quantum gases are typically done with relatively small atom numbers, from N ∼ 10 2 to N ∼ 10 6 , so that conclusions that hold in the thermodynamic limit do not necessarily apply and spin fragmentation can be observed experimentally.
The article is organized as follows. In Section 2, we present the basic model that describes an ensemble of spin 1 bosons with antiferromagnetic interactions condensing in the same orbital wave function irrespective of the internal state (single-mode alkalis (sees, e.g., [14]). approximation, or SMA). In Section 3, we use the basis of total spin eigenstates (exact in the absence of an applied magnetic field, q = 0). We derive approximate solutions for the spectrum and eigenstates for q > 0 in section 3.1, and discuss how they evolve with increasing QZ energy. Using these results, we compute in section 4 the average value and variance of N 0 at finite temperatures, and compare the approximate solution to numerical diagonalization of the Hamiltonian. We finally present in Section 5 an alternative approach, where the fragmented condensate is described as a statistical mixture of mean-field (symmetry broken) states. We find excellent agreement with the exact diagonalization of the Hamiltonian.  We mark three different regimes in the q − T plane. "Spin fragmentation" refers to a fragmented spin state with large population fluctuations, where ∆n 0 ∼ 1. In the quantum regime (N q/U s 1/N and k B T /U s 1/N ), this is due to quantum fluctuations: the system is then close to the singlet ground state. In the thermal regime(k B T N q, U s /N ), on the other hand, thermal fluctuations dominate over the quantum one and over the effect of the QZ energy (QZ energy). Conversely, "BEC m=0" refers to atoms forming a regular polar condensate with almost all atoms in m = 0, and ∆n 0 1. The plot was drawn by numerically diagonalizing the Hamiltonian (1) and computing thermodynamic averages from the spectrum and eigenstates, using N = 300 particles. Note the logarithmic scales on both the horizontal and the vertical axis.

Single-mode description of spin 1 condensates
We consider a gas of ultracold spin 1 bosons in a trap with Zeeman components m = −1, 0,or +1. We discuss the case of antiferromagnetic interactions and assume the validity of the SMA, i.e. that all bosons condense in the same spatial orbital irrespective of their internal state [25]. The Hamiltonian is [26] where U s > 0 is the spin interaction energy per atom ¶, q > 0 is the QZ energy,Ŝ is the total spin operator, andN α is the number operator in the Zeeman state α. We assume that the number of atoms N is even for simplicity. Odd values of N could be treated in a similar way, without modifying the final results to order 1/N . Typical experimental values for the parameters of the SMA model are N = 10 3 − 10 5 , U s /k B ∼ 2 − 5 nK, while q can be varied from zero to values much larger than U s by changing the magnetic field [23,24].
In the absence of an external magnetic field (q = 0), the Hamiltonian reduces to a quantum rotor with moment of inertia N/U s [5,10]. The energy eigenstates are thus simply the total spin eigenstates |N, S, M , with S the spin quantum number and M its projection on the z axis. The corresponding eigenvalues are E(S) = (U s /2N )S(S + 1), with a degeneracy 2S + 1. The wave functions for these states are known explicitly in the Fock basis [5,6,2] (see also Appendix A).
When q = 0, since [Ŝ z ,N 0 ] = 0, the magnetic quantum number M (eigenvalue of S z ) remains a good quantum number. One can diagonalizeĤ by block in each M sector. For each M , the energy eigenstates can be expressed in the angular momentum basis, To express the Hamiltonian in (1) with E the energy eigenvalue and where the coefficients h M S,S are easily obtained from the expressions given in Appendix A.

Spectrum and eigenstates for M = 0
A first approach for finding the spectrum and eigenstates is to diagonalize numerically the matrix h M in (3). Our goal this Section is to propose an analytical approximation to understand better the structure of the spectrum and eigenstates. The discussion allows one to describe how the ground state evolves with q, and will also be useful to understand qualitatively the behavior of the systems at finite temperatures later in this paper. For simplicity, we focus in this Section on the M = 0 sector. The conclusions we obtain remain qualitatively correct for M = 0 provided its value is not too large (|M | N ).

Continuum approximation for large q
We make the assumption that the thermodynamic behavior is dominated by states, such that the dominant coefficients in the |N, S, M basis obey 1 S N . As we will see later in this paper, this assumption is justified for large enough q at T = 0, and for any q at finite temperatures k B T U s /N . In this limit, the matrix elements h S,S , h S,S±2 can be simplified. We obtain to lowest order in 1/S, S/N (see Appendix B), where we have set x = S/N , = 2/N , c(x) = c S,0 . This equation maps the spin problem to a tight-binding model for a particle hopping on a lattice, with an additional harmonic potential keeping the particle near x = 0. The model is characterized by an inhomogeneous tunneling parameter J(x) = N q(1 − x 2 /2)/4 and a harmonic potential strength N U s . Boundary conditions confine the particle to 0 ≤ x ≤ 1.
If c(x) changes smoothly as a function of x, the tight-binding model can be further simplified in a continuum approximation. We show in Appendix B that the tight-binding equation reduces to the one for a fictitious one-dimensional harmonic oscillator, The boundary condition c(0) = 0 selects eigenstates of the standard harmonic oscillator with odd parity. The mass m and oscillation frequency ω of the fictitious oscillator are found from 2 /2m ≡ q/N and mω 2 ≡ N (q + 2U s )/2. The oscillator frequency is thus This collective spectrum was also obtained by the Bogoliubov approach of [9,11].

Ground state
In this Section, we use the results established previously to examine the evolution of the ground state with increasing q. Our results reproduce the ones from [11] obtained using a different method. The ground state of the truncated fictitious harmonic oscillator (with boundary condition c 0 (0) = 0) is given by with the quantum harmonic oscillator size The continuum approximation is valid only if c(x) varies smoothly on the scale of the discretization step , or equivalently when σ 1/N . This gives the validity criterion for this approximation, For q < U s /N 2 , the ground state is very close to the singlet state, with a width σ 1/N . Here spin fragmentation occurs purely due to quantum spin fluctuations (related to antiferromagnetic interactions) of a polar BEC . We indicate this state in Figure 1 as "quantum spin fragmented".
For q U s /N 2 , the continuum approximation is valid. We see from (8) that as q increases, the QZ energy mixes an increasing number of S states. Asymptotically, for q U s , the true ground state is a superposition of ∼ N σ ≈ √ 2N total spin eigenstates. In this regime, we can compute the moments of N 0 by expressing the depletion operator N −N 0 in terms of the ladder operatorsb andb † associated with the fictitious harmonic oscillator. We find For U s /N 2 q U s , the depletion N − N 0 and variance ∆N 2 0 are larger than unity but small compared to N, N 2 , respectively, while for q U s , they become less than one particle : in the latter case, the ground state approaches the Fock state â † 0 N |vac expected from mean field theory. We indicate both regimes as "BEC m=0" in Figure 1, without marking the distinction.

Excited states for M = 0
We now turn to the description of excited states, still limiting ourselves to the case M = 0 for simplicity. The tight-binding model (4) is characterized by a tunneling parameter J = N q(1 − x 2 /2)/4 and a harmonic potential strength κ = N U s . Let us examine two limiting cases. For q = 0 (no hopping), the energy eigenstates coincide with "position" eigenstates with energy E(S) ≈ (U s /2N )S 2 for S 1. Conversely, when U s = 0 the energy eigenstates are delocalized states, which form an allowed energy band of width ∼ 4J ∼ N q. The weak inhomogeneity of the tunneling parameter does not play a large role since these states are confined near x = 0 by the harmonic potential.
For the general case where J, κ = 0, the eigenstates can be divided in two groups [27,28], • low-energy states with energy E < 4J, which are extended "Bloch-like" states modified by the harmonic potential; the continuum approximation introduced earlier corresponds to an effective mass approximation, valid for low-energy states with E 4J = N q (the requirement q U s /N 2 found before still holds).
• high-energy states with E 4J, that would be in the band gap in the absence of the potential energy term (and thus forbidden).They are better viewed as localized states, peaked around x(E) ≈ 2E/N U s with a width ∼ 1/N . As a result they are very similar to the angular momentum eigenstates |N, S, 0 for the corresponding value of S. For these states, the continuum approximation does not hold.
We illustrate this classification in Figure 2, where we show the probability densities | c(S) | 2 as a function of energy. One can see the change from a "delocalized" regime at small energies to a "localized" regime at large energies. The wave functions were calculated exactly by diagonalizing the Hamiltonian 1 for N = 1000. We also show the corresponding energy spectrum in Figure 3, showing the same crossover from delocalized states at low energies to localized states at high energies. For low energies, the spectrum is given by the harmonic oscillator model, n ≈ ω(2n + 3/2) with n integer. For high energies, the energy eigenstates are localized around x n = n/N , with a spectrum given by n ≈ U s n 2 /2N with n integer. Both expressions agree well with the numerical result in their respective domains of validity.

Spin fragmentation at finite temperatures
We have seen in Section 3.2 that for a system in its ground state, the depletion and fluctuations of the M = 0 state were rapidly collapsing as q was increased above U s /N 2 , and the system turned from a fragmented to a single condensate with all atoms in the Zeeman state m = 0. The energy gap to the first excited state is 3U s /N near q = 0. For typical experimental values [23,24], this corresponds to a few pK, vastly smaller than realistic temperatures for a typical experiment (a few tens of nK) due to the 1/N scaling. Therefore, it is natural to ask how the crossover from a fragmented to a single condensate is modified at finite temperatures. In the remainder of the paper, we thus consider the high temperature case k B T U s /N . We will compute the first two moments of N 0 at finite temperatures, N 0 T and (∆N 2 0 ) T = N 2 0 T − N 0 2 T , and use these quantities to study the fragmented to single condensate crossover.

Spin fragmentation for q = 0
Let us first consider the case q = 0. An important remark is that super-Poissonian fluctuations are not unique to the ground state, but also occur for low-energy eigenstates with S N . This is best seen by considering values of S such that 1 S N . In this limit, we find where N p 0 SM = N, S, M | N p 0 | N, S, M . Hence, we find super-Poissonian fluctuations for M S N , which eventually vanish as S (resp. M ) increases to its maximum value N (resp. S).
We calculate now the thermally averaged n 0 T and (∆n 2 0 ) T in the canonical ensemble. The average population in m = 0 is given by Spin fragmentation of Bose-Einstein condensates with antiferromagnetic interactions10 The second moment N 2 0 T and the variance (∆N 2 0 ) T are given by similar expressions. Here Z is the partition function and β = U s /2N k B T . Assuming that the temperature is large compared to the level spacing (k B T U s /N ), the thermodynamic sums over energy levels is dominated by states with large S 1. There are two regimes to consider.
At intermediate temperatures, states with 1 S N dominate the thermodynamics. To calculate the thermal average over all S in this regime, we replace the discrete sums by integrals and send the upper bound N of the integral to infinity. A simple estimate of the mean value of S, S ∼ (N k B T /U s ) 1/2 , shows that the condition 1 S N corresponds to the boundaries In this regime, we find We note that to leading order in 1/N , the moments of N 0 are identical for those found in the singlet state. The second regime arises when the temperature becomes very large (k B T /N U s > 1), where one expects the sum to saturate due to the finite number of states. In this limit, the upper bound of the integral cannot be taken to infinity, and one must take the restriction S ≤ N into account. On the other hand, the Boltzmann factor can be replaced by unity, and the sums can then be calculated analytically. One finds To summarize (see Figure 1), for q = 0 we always find large depletion and super-Poissonian fluctuations (∆N 2 0 ∼ N 0 2 ). The average population is always N/3 as expected from the isotropy of the Hamiltonian. The relative standard deviation remains approximately constant (to order N ) at the value ∆N 0 /N ≈ 2/3 √ 5 ≈ 0.298 for k B T N U s , and changes to 1/3 √ 2 ≈ 0.236 for very large temperatures k B T > N U s where all states are occupied with equal probability.

Bogoliubov approximation for q = 0
For large q > 0 (and S z constrained to vanish only in average), we expect that the system will form a condensate in the m = 0 Zeeman state, with small fluctuations. Such a system can be described in the Bogoliubov approximation (as described in the Appendix of [11]), which extends to any M the harmonic oscillator approximation made earlier for the M = 0 sector. One setsâ 0 ≈ √ N 0 , and expresses the fluctuationsâ ±1 in terms of new operatorsα ± , Spin fragmentation of Bose-Einstein condensates with antiferromagnetic interactions11 Here the Bogoliubov amplitudes u, v defined by are chosen to put the Hamiltonian in diagonal form, The energy ω of the Bogoliubov mode is identical to the one previously found in the harmonic oscillator approximation for M = 0 [Eq. (6)]. Note that we have now two such modes (instead of only one in the case M = 0) + . In the Bogoliubov approximation, the moments of N 0 can be obtained analytically. The quantum (T = 0) depletion of N 0 is smaller than one atom. The thermal part of the depletion and variance of n 0 = N 0 /N read for k B T ω The prefactors take values of order unity, and both the depletion 1 − n 0 and standard deviation ∆n 0 scale as k B T /N q. The above expressions are valid provided they describe small corrections to a regular polar condensate where almost all atoms accumulate in m = 0 ( n 0 = 1), or in other words for temperatures

Comparison between the different approximations
We compare in Figure 4 the predictions for the moments of N 0 obtained from the various approximations discussed the paper, Bogoliubov approximation, and q = 0 limit. These approximations are compared to the results obtained by diagonalization of the original Hamiltonian (1) and computing thermodynamic averages using the exact spectrum and eigenstates. When N q/k B T 1, the localized states of Section 3.3, which are dominated by their potential energy, will be populated. Because these localized states are close to the angular momentum eigenstates found in the q = 0 limit, to a good approximation the formula derived in Section 4.1 [see (17,18) and the continous blue line in Figure 4]. On the other hand, for N q/k B T 1, thermal states mostly populate states with E ∼ qN , i.e. "delocalized" states within the low-energy "Bloch band" of width ∼ N q. Those states correspond to small depletion and fluctuations, and they are well described by + We expect in general three modes of excitations for a spin 1 system. When the constraint of constant particle number is taken into account, this reduces the number of modes to two. The suppressed mode would correspond to density fluctuations in an extended system, and is explicitly ruled out by the SMA. When a further constraint M = 0 is imposed, another mode is cancelled -corresponding to magnetization fluctuations which are explicitly forbidden, thus leaving only one excitation mode.
the Bogoliubov approximation presented in Section 4.2 [see (23,24) and the red dashed line in Figure 4]. The numerical solution of the original model (3) interpolates between the two well-defined asymptotic limits, either a thermal mixture of total spin eigenstates for q k B T /N or a thermal state of Bogoliubov-like excitations for q k B T /N . We note to conclude this section that in the regime U s /N k B T N U s , the tightbinding model defined Eq. (3) has a quasi-universal form at finite temperatures, in the sense that the model is entirely specified by two dimensionless parameters, for instance k B T /U s and N q/U s . We found that the physical quantities N 0 , (∆N 0 ) depend only on their ratio N q/k B T , to a very good approximation. This quasi-universality, which can be explored by experiments, will be easily justified in the broken symmetry approach presented in the next Section.  (17,18), and the red dashed line shows the Bogoliubov approximation. Deviations are observed for q/T ∼ N , which is expected from our approximation: This regime corresponds to a depletion of one atom or less, and corrections ∝ 1/N that we neglect become important.

Comparison with the broken-symmetry picture
So far, we have treated the problem by the most natural method, by looking for the eigenspectrum of the Hamiltonian. Another approach [2,3] to the problem of spin 1 bosons with antiferromagnetic interactions relies on the set of so-called polar or spinnematic states, defined as where the vector Ω reads in the standard basis Spin fragmentation of Bose-Einstein condensates with antiferromagnetic interactions13 For a single particle, the states |Ω = i=0,±1 Ω i |m = i form a continuous family of spin 1 wavefunctions with vanishing average spin. In fact, |Ω is the eigenvector with zero eigenvalue of the operator Ω ·ŝ, withŝ the spin 1 operator. The states |N : Ω correspond to a many-body wave function where all particle occupy the single-particle state |Ω . As a result, one has N : Ω|Ŝ|N : Ω = 0.

Zero temperature
It is interesting to connect the spin nematic states to the angular momentum eigenstates.
where Y SM denotes the usual spherical harmonics and where dΩ = sin(θ)dθdφ. In particular, the singlet ground state |N, 0, 0 appears to be a coherent superposition with equal weights of the nematic states. Consider now the average value in the singlet state Ô k singlet = N, 0, 0|Ô k |N, 0, 0 of a k−body operatorÔ k , As shown in [2], for few-body operators with k N this expectation value can be approximated to order 1/N by the much simpler expression This approximation shows that the system can equally well be described by a statistical mixture of spin-nematic states described by the density matrix [2,3] ρ BS = 1 4π dΩ|N : Ω N : Ω|.
At zero temperature and zero field, there is no preferred direction for the vector Ω so that each state can appear with equal probability. This approach is known as a "broken symmetry" point of view, where one can imagine that the atoms condense in the same spin state for each realization of the experiment, but this spin state fluctuates arbitrarily from one realization to the next. The important point is that the overlap integral N : Ω |N : Ω between two spin-nematic states vanishes very quickly with the distance |Ω − Ω |. This allows one to use the approximation N : Ω |N : Ω ≈ δ(Ω − Ω ) + O(1/N ), which leads to This result can be written as a general statement concerning average values of few-body observables with k N [2]: to leading order in 1/N , the exact and broken symmetry approaches will give the same results after averaging over the ensemble. The differences between the two approaches are subtle and vanish in the thermodynamic limit as 1/N .
The variance in the ensemble is thus super-Poissonian, and differs from the result in the exact ground state only by the sub-leading term ∝ N . This is in agreement with the general statement made above.

Moments of N 0 at finite temperatures
We now extend the broken symmetry approach summarized above to finite temperatures. The density matrix should include a weight factor proportional to the energy of the states |N : Ω . To leading order in 1/N , these states have zero interaction energy * and a mean QZ energy given by −N q cos 2 (θ). In the spirit of the mean-field approximation, we replace the Boltzmann factor by its mean value and write the density matrix asρ with β = 1/k B T . The partition function can then be expressed as Here we introduced the family of functions which are related to the lower incomplete gamma functions. In a similar way, we can compute the moments of n 0 = N 0 /N to leading order in N as From this result, one can easily deduce the average and variance of n 0 . This calculation provides an explicit proof of the numerical evidence that, to leading order in N , the moments of N 0 obey a universal curve depending only on N q/k B T and not on q/U s or T /U s separately. From the properties of the functions F α , we recover the results established in the previous Section. When x → 0, F α (x) ∼ x α /α. Using this result we recover for q = 0 the previous results, i.e. n 0 = 1/3, n 2 0 = 1/5 and ∆n 2 0 = 4/45. When x → ∞, . This leads to the asymptotic behavior n m 0 ∼ 1 − m/(N βq) + (m 2 − 3m/2)/(N βq) 2 + · · · when N βq 1, which reproduces the Bogoliubov results (23,24) for q U s . We finally compare in Figure 5 the results from the broken symmetry approach to the results obtained by diagonalizing the Hamiltonian 1. We find excellent agreement between the two in the regime of thermal fragmentation, supporting the picture of mean-field states with random orientation fluctuating from one realization to the next. We note that the ansatz (33) for the density matrix is by no means obvious, and the good agreement with the numerical results is obtained only because the set of polar states is a good description for sufficiently low temperatures : Although these states are not true eigenstates of the Hamiltonian (1), the action ofĤ yields off-diagonal matrix elements scaling as 1/N [29], and thus vanishing in the thermodynamic limit. At high temperatures (k B T ∼ N U s ), where all high energy states are populated the broken-symmetry ansatz is no longer adequate.

Conclusion
We have studied the properties of an ensemble of antiferromagnetic spin 1 bosons with QZ energy breaking the spin rotational symmetry. The system evolves with increasing QZ energy from a super fragmented condensate with large fluctuations to a regular polar condensate where atoms condense in m = 0. We focused in particular on the behavior of a thermal mixture of excited states, and discussed the evolution of the moments of N 0 with increasing q. Two approaches were explored, one relying on diagonalization of the Hamiltonian (either exactly or approximately in certain parameter regimes), and the other relying on a broken symmetry picture where the system is described as a statistical mixture of degenerate polar condensates. Both approaches were found in remarkable agreement. In this article, we focused on equilibrium properties and assumed thermal equilibrium from the start. An interesting question is how the physical system (i.e. also including the dynamics of non-condensed modes not described in the SMA) can reach such an equilibrium state, e.g. following a quench in q [30]. This problem, which can be linked to the more general question of thermalization of closed quantum systems [31] provides an interesting direction for future work.

Appendix B. Continuum approximation
We expand the matrix elements h S,S , h S,S±2 to first order in 1/S, S/N, M 2 /S 2 , and obtain For M = 0, we obtain where we have set x = S/N and = 2/N . We now take the continuum limit, where 1 is taken as a discretization step and c s becomes a continuous function c(S). We write N 2 4 (c S+2 + c S−2 ) ≈ ∆c(s) + N 2 2 c(s).

(B.4)
Substituting in (B.4) and neglecting a term ∝ (qx 2 /N )∆c, we arrive at (5). This derivation is valid as long as the relevant states are well localized around x = 0. This is always the case in the ground state, which has a width at most ∼ 1/ √ N for q U s . For thermal states, the width is ∼ k B T / [N (2U s + q)], which gives the condition k B T N (U s + q). Finally, the cross-term ∝ (qx 2 /N )∆c is of order 2E p E c c/ [N (2U s + q)] in terms of the kinetic and potential energies E c , E p of the harmonic oscillator. In the thermal regime, a typical order of magnitude for this term is thus (k B T ) 2 /[N (2U s + q)], small compared to the energy k B T typical for the other terms we kept in the equation provided the condition above is fulfilled.
We acknowledge discussions with members of the LKB, in particular Yvan Castin. This work was supported by IFRAF, by Ville de Paris (Emergences project) and by DARPA (OLE project).