# Spin fragmentation of Bose-Einstein condensates with antiferromagnetic interactions

Luigi De Sarlo, Lingxuan Shao, Vincent Corre, Tilman Zibold, David Jacob,

Jean Dalibard, Fabrice Gerbier

► To cite this version:

Luigi De Sarlo, Lingxuan Shao, Vincent Corre, Tilman Zibold, David Jacob, et al. Spin fragmentation of Bose-Einstein condensates with antiferromagnetic interactions. New Journal of Physics, 2013, 15, pp.113039. 10.1088/1367-2630/15/11/113039. hal-01587034

HAL Id: hal-01587034

https://hal.sorbonne-universite.fr/hal-01587034

Submitted on 13 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Spin fragmentation of Bose-Einstein condensates
with antiferromagnetic interactions

To cite this article: Luigi De Sarlo et al 2013 *New J. Phys.* **15** 113039

View the [article online](#) for updates and enhancements.

## Related content

- [PhD Tutorial: Quantum field theory for trapped Bose-Einstein condensates](#)
J Rogel-Salazar, S Choi, G H C New et al.

- [Spin-1 condensates at thermal equilibrium: A SU(3) coherent state approach](#)
V. Corre, T. Zibold, C. Frapolli et al.

- [Fragmented condensation in Bose-Hubbard trimers with tunable tunnelling](#)
A Gallemi, M Guilleumas, J Martorell et al.

## Recent citations

- [Shortcut to adiabaticity in spinor condensates](#)
Arnau Sala et al

- [Competition between Bose-Einstein Condensation and Spin Dynamics](#)
B. Naylor et al

- [Spatiotemporal coherence of non-equilibrium multimode photon condensates](#)
Jakov Marelic et al
# Spin fragmentation of Bose–Einstein condensates with antiferromagnetic interactions

Luigi De Sarlo<fnref n="1" />, Lingxuan Shao, Vincent Corre, Tilman Zibold, David Jacob, Jean Dalibard and Fabrice Gerbier<fnref n="2" />

Laboratoire Kastler Brossel, CNRS, ENS, UPMC, 24 rue Lhomond, F-75005 Paris, France
E-mail: fabrice.gerbier@lkb.ens.fr

New Journal of Physics 15 (2013) 113039 (18pp)

Received 2 July 2013

Published 19 November 2013

Online at http://www.njp.org/

doi:10.1088/1367-2630/15/11/113039

**Abstract.** 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 Zeeman populations, into a regular polar condensate, where the atoms all condense in the $m=0$ state along the field direction. We calculate the average value and variance of the population 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 the spin fluctuations are a mesoscopic effect that will not survive in the thermodynamic limit $N \to \infty$, but are observable for a sufficiently small atom number.

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹ Current address: SYRTE, Observatoire de Paris, CNRS, UPMC; 61 Avenue de l'Observatoire, 75014 Paris, France.</footnote>

<footnote n="2">² Author to whom any correspondence should be addressed.</footnote>
Contents

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Single-mode description of spin-1 condensates</td><td>5</td></tr><tr><td>3. Spectrum and eigenstates for M = 0</td><td>5</td></tr><tr><td>3.1. Continuum approximation for large q</td><td>6</td></tr><tr><td>3.2. Ground state</td><td>6</td></tr><tr><td>3.3. Excited states for M = 0</td><td>7</td></tr><tr><td>4. Spin fragmentation at finite temperatures</td><td>8</td></tr><tr><td>4.1. Spin fragmentation for q = 0</td><td>9</td></tr><tr><td>4.2. Bogoliubov approximation for q ≠ 0</td><td>10</td></tr><tr><td>4.3. Comparison between the different approximations</td><td>11</td></tr><tr><td>5. Comparison with the broken-symmetry picture</td><td>12</td></tr><tr><td>5.1. Zero temperature</td><td>12</td></tr><tr><td>5.2. Moments of N₀ at finite temperatures</td><td>14</td></tr><tr><td>6. Conclusion</td><td>15</td></tr><tr><td>Acknowledgments</td><td>15</td></tr><tr><td>Appendix A. Total spin eigenstates</td><td>16</td></tr><tr><td>Appendix B. Continuum approximation</td><td>17</td></tr><tr><td>References</td><td>17</td></tr></table>

1. Introduction

The natural behavior of bosons at low enough temperatures is to form a Bose-Einstein condensate (BEC), 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 (pseudo-spin-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 \cdot 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 [2, 5, 6], the signature of fragmentation is then the occurrence of anomalously large fluctuations of the populations $N_m$ in the Zeeman states $m = 0, \pm 1$ (see also [7, 8]), where 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 $\langle N_0 \rangle = N/3$ and $\Delta N_0^2 \approx 4N^2/45$, respectively ($N$ is the total number of particles). Such super-Poissonian fluctuations ($\Delta N_0^2 \propto \langle N_0 \rangle^2$) deviate strongly from the value expected for a single condensate or any ensemble without correlations where
$\Delta N_0^2 \propto \langle N_0 \rangle$.<fnref n="3" /> It was pointed out by Ho and Yip [6] that such a state was probably not realized in typical experiments, due to its fragility toward any perturbation breaking spin rotational symmetry (see also [9–14]). In the thermodynamic limit $N \to \infty$, an arbitrary small symmetry-breaking perturbation is enough to favor a regular condensed state, where almost all the atoms occupy the same (spinor) condensate wave function and $\Delta N_0 \ll N$.

In this paper, we give a detailed analysis of the phenomenon of spin fragmentation for spin-1 bosons. Our analysis assumes the conservation of 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$.<fnref n="4" /> The quadratic Zeeman (QZ) energy breaks the spin rotational symmetry, and favors a condensed state with $m=0$ along the field direction. In [10–13], 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 paper 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 of the $m=0$ state are observed for small $q$. We can distinguish three different regimes. For low $q \ll U_s/N^2$ and low temperatures $k_B T \ll U_s/N$ ($U_s \propto 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’ [2, 5, 6, 12]. We also observe a thermal regime for $k_B T \gg Nq$, $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 [16, 17]. In this limit, $\langle N_0 \rangle \approx N$ and $\Delta N_0 \ll N$. We indicate this third regime as ‘BEC $m=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 [18]. 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 [19–24]. 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

<footnote n="3">³ 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.</footnote>

<footnote n="4">⁴ 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 alkalis (see e.g. [15]).</footnote>
<figure><img src="image_3.png" /></figure>

**Figure 1.** Standard deviation $\Delta n_0 = \Delta N_0/N$ of the population $N_0$ of the Zeeman state $m=0$, normalized to the total atom number $N$. We mark three different regimes in the $q-T$ plane. ‘Spin fragmentation’ refers to a fragmented spin state with large population fluctuations, where $\Delta n_0 \sim 1$. In the quantum regime ($Nq/U_s \ll 1/N$ and $k_B T/U_s \ll 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 \gg Nq$, $U_s/N$), on the other hand, thermal fluctuations dominate over the quantum one and over the effect of the QZ energy. Conversely, ‘BEC $m=0$’ refers to the atoms forming a regular polar condensate with almost all the atoms in $m=0$, and $\Delta n_0 \ll 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.

is that the present model of antiferromagnetic spin-1 BECs is expected to accurately describe actual experimental systems [25, 26]. 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 conducted with relatively small atom numbers, from $N \sim 10^2$ to $10^6$, so that the conclusions that hold in the thermodynamic limit do not necessarily apply and spin fragmentation can be observed experimentally.

The paper 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 approximation (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. By 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.

## 2. 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 [27]. The Hamiltonian is [28]

$$ \hat{H} = \frac{U_s}{2N} \hat{S}^2 - q \hat{N}_0, \quad (1) $$

where $U_s > 0$ is the spin interaction energy per atom$^5$, $q > 0$ is the QZ energy, $\hat{S}$ is the total spin operator and $\hat{N}_\alpha$ is the number operator in Zeeman state $\alpha$. 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 \sim 2-5$ nK, while $q$ can be varied from zero to values much larger than $U_s$ by changing the magnetic field [25, 26].

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, 11]. The energy eigenstates are thus simply the total spin eigenstates $|N, S, M\rangle$, 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 [2, 5, 6] (see also appendix A).

When $q \neq 0$, since $[\hat{S}_z, \hat{N}_0] = 0$, the magnetic quantum number $M$ (eigenvalue of $\hat{S}_z$) remains a good quantum number. One can diagonalize $\hat{H}$ by block in each $M$ sector. For each $M$, the energy eigenstates can be expressed in the angular momentum basis

$$ |\phi_M\rangle = \sum_{S=|M|}^{N} c_{S,M} |N, S, M\rangle. \quad (2) $$

To express the Hamiltonian in (1) in the $|N, S, M\rangle$ basis, we need to compute the action of $\hat{N}_0$. The non-vanishing matrix elements of $\hat{N}_0$ are $\langle N, S, M | \hat{N}_0 | N, S, M\rangle$, $\langle N, S \pm 2, M | \hat{N}_0 | N, S, M\rangle$ (see appendix A). The Schrödinger equation then takes the form of a tridiagonal matrix equation

$$ h_{S,S+2}^{M} c_{S+2,M} + h_{S,S-2}^{M} c_{S-2,M} + h_{S,S}^{M} c_{S,M} = E c_{S,M} \quad (3) $$

with $E$ the energy eigenvalue and where the coefficients $h_{S,S'}^M$ are easily obtained from the expressions given in appendix A.

## 3. 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 in this section is to propose an analytical approximation to better understand the structure of the spectrum and the eigenstates. The discussion allows one to describe how the

<footnote>$^5$ The spin interaction energy $U_s$ can be calculated from $U_s = g_s \int dx |\psi(x)|^4$, where $\psi(x)$ is the spatial orbital of the condensate.</footnote>
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 \neq 0$ provided its value is not too large ($|M| \ll N$).

### 3.1. 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\rangle$ basis obey $1 \ll S \ll 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 \gg U_s/N$. In this limit, the matrix elements $h_{S,S}$, $h_{S,S\pm2}$ can be simplified. We obtain to lowest order in $1/S$, $S/N$ (see appendix B)

$$ -J(x+\epsilon)c(x+\epsilon) - J(x-\epsilon)c(x-\epsilon) + \frac{NU_s}{2}x^2c(x) = \left(E + \frac{Nq}{2}\right)c(x), \quad (4) $$

where we have set $x=S/N$, $\epsilon=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) = Nq(1-x^2/2)/4$ and a harmonic potential strength $NU_s$. Boundary conditions confine the particle to $0 \le x \le 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

$$ -\frac{q}{N}c''(x) + \frac{N}{4}(q+2U_s)x^2c(x) = (E+Nq)c(x). \quad (5) $$

The boundary condition $c(0) = 0$ selects the eigenstates of the standard harmonic oscillator with odd parity. The mass $m$ and oscillation frequency $\omega$ of the fictitious oscillator are found from $\hbar^2/2m \equiv q/N$ and $m\omega^2 \equiv N(q+2U_s)/2$. The oscillator frequency is thus

$$ \hbar\omega = \sqrt{q(q+2U_s)}. \quad (6) $$

This collective spectrum was also obtained by the Bogoliubov approach of [10, 12].

### 3.2. 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 [12] obtained by using a different method. The ground state of the truncated fictitious harmonic oscillator (with boundary condition $c_0(0) = 0$) is given by

$$ c_0(x) = \frac{1}{\pi^{1/4}\sigma^{1/2}} \frac{x}{\sigma} \exp\left(-\frac{x^2}{2\sigma^2}\right) \quad (7) $$

with the quantum harmonic oscillator size

$$ \sigma = \sqrt{\frac{\hbar}{m\omega}} = \sqrt{\frac{2}{N}} \left( \frac{q}{q+2U_s} \right)^{1/4}. \quad (8) $$
The continuum approximation is valid only if $c(x)$ varies smoothly on the scale of the discretization step $\epsilon$, or equivalently when $\sigma \gg 1/N$. This gives the validity criterion for this approximation

$$q \gg \frac{U_s}{N^2} \qquad (9)$$

For $q < U_s/N^2$, the ground state is very close to the singlet state, with a width $\sigma \ll 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 \gg 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 \gg U_s$, the true ground state is a superposition of $\sim N\sigma \approx \sqrt{2N}$ total spin eigenstates. In this regime, we can compute the moments of $N_0$ by expressing the depletion operator $N - \hat{N}_0$ in terms of the ladder operators $\hat{b}$ and $\hat{b}^\dagger$ associated with the fictitious harmonic oscillator. We find

$$N - \langle N_0 \rangle = \frac{U_s + q}{2\sqrt{q(q + 2U_s)}}, \qquad (10)$$

$$\Delta N_0^2 = \frac{U_s^2}{2q(q + 2U_s)}. \qquad (11)$$

For $U_s/N^2 \ll q \ll U_s$, the depletion $N - \langle N_0 \rangle$ and variance $\Delta N_0^2$ are larger than unity but small compared to $N, N^2$, respectively, while for $q \gg U_s$, they become less than one particle: in the latter case, the ground state approaches the Fock state $(\hat{a}_0^\dagger)^N |\text{vac}\rangle$ expected from the mean field theory. We indicate both regimes as ‘BEC $m = 0$’ in figure 1, without marking the distinction.

### 3.3. Excited states for $M = 0$

We now turn to the description of the excited states, still limiting ourselves to the case $M = 0$ for simplicity. The tight-binding model (4) is characterized by a tunneling parameter $J = Nq(1 - x^2/2)/4$ and a harmonic potential strength $\kappa = N U_s$. Let us examine two limiting cases. For $q = 0$ (no hopping), the energy eigenstates coincide with the ‘position’ eigenstates with energy $E(S) \approx (U_s/2N)S^2$ for $S \gg 1$. Conversely, when $U_s = 0$ the energy eigenstates are delocalized states, which form an allowed energy band of width $\sim 4J \sim Nq$. 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, \kappa \neq 0$, the eigenstates can be divided into two groups [29, 30]:

* 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 \ll 4J = Nq$ (the requirement $q \gg U_s/N^2$ found before still holds); and

* high-energy states with $E \gg 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) \approx \sqrt{2E/NU_s}$ with a width $\sim 1/N$. As a result, they are very similar to the angular momentum eigenstates $|N, S, 0\rangle$ for the corresponding value of $S$. For these states, the continuum approximation does not hold.
<figure><img src="image_4.png" /></figure>

**Figure 2.** Probability densities (amplitude shown as gray scale) of the eigenstates of the spin-1 Hamiltonian (1) as a function of ‘position’ $x = S/N$ and energy $E/U_s$. The plot corresponds to $N = 2000$ and $Nq/U_s = 10$. At low energy, the eigenstates explore the whole available region, from the turning point down to $x = 0$. Conversely, at high energies, the eigenstates are more and more localized around the diagonal, as expected for potential energy eigenstates. We show as insets the probability densities for the 10th (b) and 500th (c) excited states for illustration.

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 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, $\epsilon_n \approx \hbar\omega(2n + 3/2)$ with $n$ integer. For high energies, the energy eigenstates are localized around $x_n = n/N$, with a spectrum given by $\epsilon_n \approx U_s n^2/2N$ with $n$ integer. Both expressions agree well with the numerical result in their respective domains of validity.

**4. 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 the 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 [25, 26], 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 \gg U_s/N$. We will compute the first two moments of $N_0$ at finite temperatures, $\langle N_0 \rangle_T$ and $(\Delta N_0^2)_T = \langle N_0^2 \rangle_T - \langle N_0 \rangle_T^2$, and use these quantities to study the fragmented to single condensate crossover.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Energy spectrum for $N = 2000$ and $Nq/U_s = 10$. The black solid line is the spectrum calculated by numerical diagonalization of the Hamiltonian (1), shifted up by $qN$. The red dashed line corresponds to $E(S) = (U_s/2N)S^2 + qN/2$, the blue dotted line to the harmonic oscillator approximation.</figcaption></figure>

### 4.1. 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 \ll N$. This is best seen by considering values of $S$ such that $1 \ll S \ll N$. In this limit, we find

$$ \langle \hat{N}_0 \rangle_{SM} \approx (N^2 - S^2)(S^2 - M^2)/8, \qquad (12) $$

$$ \langle \hat{N}_0^2 \rangle_{SM} \approx (3N^2 - S^2)(S^2 - M^2)^2/8, \qquad (13) $$

$$ (\Delta N_0^2)_{SM} \approx (N^2 - S^2)(S^2 - M^2)^2/8, \qquad (14) $$

where $\langle \hat{N}_0^p \rangle_{SM} = \langle N, S, M | N_0^p | N, S, M \rangle$. Hence, we find super-Poissonian fluctuations for $M \ll S \ll N$, which eventually vanish as $S$ (resp. $M$) increases to its maximum value $N$ (resp. $S$).

We calculate now the thermally averaged $\langle n_0 \rangle_T$ and $(\Delta n_0^2)_T$ in the canonical ensemble. The average population in $m = 0$ is given by

$$ \langle N_0 \rangle_T = \frac{1}{Z} \sum_{S,M} e^{-\beta'S(S+1)} \langle N_0 \rangle_{SM}. \qquad (15) $$

The second moment $\langle N_0^2 \rangle_T$ and the variance $(\Delta N_0^2)_T$ are given by similar expressions. Here, $Z$ is the partition function and $\beta' = U_s/2Nk_B T$. Assuming that the temperature is large compared to the level spacing $(k_B T \gg U_s/N)$, the thermodynamic sums over energy levels is dominated by states with large $S \gg 1$. There are two regimes to consider.

At intermediate temperatures, states with $1 \ll S \ll 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$, $\langle S \rangle \sim (Nk_B T/U_s)^{1/2}$, shows that the condition $1 \ll S \ll N$ corresponds to the boundaries

$$ \frac{U_s}{N} \ll k_B T \ll N U_s. \qquad (16) $$
In this regime, we find

$$ \langle N_0 \rangle_T \approx \frac{N}{3}, \qquad (17) $$

$$ (\Delta N_0^2)_T = \langle N_0^2 \rangle_T - \langle N_0 \rangle_T^2 \approx \frac{4N^2}{45}, \quad k_B T \ll N U_s. \qquad (18) $$

We note that to the 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/NU_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 \le 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

$$ (\Delta N_0^2)_{T \gg N U_s} \approx N^2/18. \qquad (19) $$

To summarize (see figure 1), for $q=0$ we always find large depletion and super-Poissonian fluctuations ($\Delta N_0^2 \sim \langle N_0 \rangle^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 $\Delta N_0/N \approx 2/3\sqrt{5} \approx 0.298$ for $k_B T \ll N U_s$, and changes to $1/3\sqrt{2} \approx 0.236$ for very large temperatures $k_B T > N U_s$ where all the states are occupied with equal probability.

## 4.2. Bogoliubov approximation for $q \neq 0$

For large $q > 0$ (and $\langle S_z \rangle$ 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 [12]), which extends to any $M$ the harmonic oscillator approximation made earlier for the $M=0$ sector. One sets $\hat{a}_0 \approx \sqrt{N} \hat{a}_0$, and expresses the fluctuations $\hat{a}_{\pm 1}$ in terms of new operators $\hat{a}^{\pm}$,

$$ \hat{a}_{\pm} = u \hat{a}_{\pm 1} - v \hat{a}_{\pm 1}^{\dagger}. \qquad (20) $$

Here, the Bogoliubov amplitudes $u, v$ defined by

$$ u \pm v = \left( \frac{q}{2U_s + q} \right)^{\pm 1/4} \qquad (21) $$

are chosen to put the Hamiltonian in diagonal form

$$ H_{\text{Bogo}} = \sum_{\mu=\pm} \hbar\omega \left( \hat{a}_{\mu}^{\dagger} \hat{a}_{\mu} + \frac{1}{2} \right) - (g+q). \qquad (22) $$

The energy $\hbar\omega$ of the Bogoliubov mode is identical to the one previously found in the harmonic oscillator approximation for $M=0$ (equation (6)). Note that we have now two such modes (instead of only one in the case $M=0$).$^6$

<footnote>$^6$ 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 the 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 canceled—corresponding to magnetization fluctuations which are explicitly forbidden, thus leaving only one excitation mode.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Depletion (left) and standard deviation (right) of $N_0$. The solid line shows the exact numerical result for $N = 1000$ and $T = 10U_s$, the blue dashed line is the result calculated for $q = 0$ by using equations (17), (18) and the red dashed line shows the Bogoliubov approximation. Deviations are observed for $q/T \sim N$ (not visible at the scale of the figure), which is expected from our approximation: this regime corresponds to a depletion of one atom or less, and corrections $\propto 1/N$ that we neglect become important.</figcaption></figure>

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 \gg \hbar\omega$:

$$1 - \langle n_0 \rangle = \frac{2(U_s + q) k_B T}{q + 2U_s Nq}, \qquad (23)$$

$$\Delta n_0^2 = \frac{2[(U_s + q)^2 + U_s^2]}{(q + 2U_s)^2} \left(\frac{k_B T}{Nq}\right)^2. \qquad (24)$$

The prefactors take values of order unity, and both the depletion $1 - \langle n_0 \rangle$ and standard deviation $\Delta n_0$ scale as $k_B T/Nq$. The above expressions are valid provided they describe small corrections to a regular polar condensate where almost all the atoms accumulate in $m = 0$ ($\langle n_0 \rangle = 1$), or in other words for temperatures

$$k_B T \ll Nq. \qquad (25)$$

### 4.3. Comparison between the different approximations

We compare in figure 4 the predictions for the moments of $N_0$ obtained from the various approximations discussed in 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 by using the exact spectrum and eigenstates.

When $Nq/k_B T \ll 1$, the localized states of section 3.3, which are dominated by their potential energy, will be populated. Since 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 continuous blue line in figure 4). On the other hand,
for $Nq/k_B T \gg 1$, thermal states mostly populate states with $E \sim qN$, i.e. 'delocalized' states within the low-energy 'Bloch band' of width $\sim Nq$. Those states correspond to small depletion and fluctuations, and they are well described by 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 \ll k_B T/N$ or a thermal state of Bogoliubov-like excitations for $q \gg k_B T/N$.

We note to conclude this section that in the regime $U_s/N \ll k_B T \ll N U_s$, the tight-binding model defined in equation (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 $Nq/U_s$. We found that the physical quantities $\langle N_0 \rangle$, $(\Delta N_0)$ depend only on their ratio $Nq/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.

## 5. 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 spin-nematic states, defined as

$$|N : \mathbf{\Omega}\rangle = \frac{1}{\sqrt{N!}} (\mathbf{\Omega} \cdot \hat{\mathbf{a}})^N |\text{vac}\rangle, \qquad (26)$$

where the vector $\mathbf{\Omega}$ reads in the standard basis

$$\mathbf{\Omega} = e^{i\gamma} \begin{pmatrix} \frac{1}{\sqrt{2}} \sin(\theta) e^{i\phi} \\ \cos(\theta) \\ -\frac{1}{\sqrt{2}} \sin(\theta) e^{-i\phi} \end{pmatrix}. \qquad (27)$$

For a single particle, the states $|\mathbf{\Omega}\rangle = \sum_{i=0,\pm 1} \mathbf{\Omega}_i |m=i\rangle$ form a continuous family of spin 1 wave functions with vanishing average spin. In fact, $|\mathbf{\Omega}\rangle$ is the eigenvector with zero eigenvalue of the operator $\mathbf{\Omega} \cdot \hat{\mathbf{s}}$, with $\hat{\mathbf{s}}$ the spin 1 operator. The states $|N : \mathbf{\Omega}\rangle$ correspond to a many-body wave function where all the particles occupy the single-particle state $|\mathbf{\Omega}\rangle$. As a result, one has $\langle N : \mathbf{\Omega}|\hat{\mathbf{s}}|N : \mathbf{\Omega}\rangle = 0$.

### 5.1. Zero temperature

It is interesting to connect the spin nematic states to the angular momentum eigenstates. The spin nematic states form an overcomplete basis of the bosonic Hilbert space. On writing the states $|N, S, M\rangle$ in this basis, one finds [2, 3, 12]

$$|N, S, M\rangle \propto \int d\mathbf{\Omega} Y_{S,M}(\mathbf{\Omega}) |N : \mathbf{\Omega}\rangle, \qquad (28)$$

where $Y_{S,M}$ denotes the usual spherical harmonics and where $d\mathbf{\Omega} = \sin(\theta)d\theta d\phi$. In particular, the singlet ground state $|N, 0, 0\rangle$ appears to be a coherent superposition with equal weights
of the nematic states. Consider now the average value in the singlet state $\langle \hat{O}_k \rangle_{\text{singlet}} = \langle N, 0, 0 | \hat{O}_k | N, 0, 0 \rangle$ of a k-body operator $\hat{O}_k$,

$$ \langle \hat{O}_k \rangle_{\text{singlet}} = \frac{1}{4\pi} \int d\Omega d\Omega' \langle N : \Omega' | \hat{O}_k | N : \Omega \rangle. \qquad (29) $$

As shown in [2], for a few-body operators with $k \ll N$ this expectation value can be approximated to order $1/N$ by the much simpler expression

$$ \langle \hat{O}_k \rangle_{\text{singlet}} \approx \frac{1}{4\pi} \int d\Omega \langle N : \Omega | \hat{O}_k | N : \Omega \rangle. \qquad (30) $$

This approximation shows that the system can be equally well described by a statistical mixture of spin-nematic states described by the density matrix [2, 3]

$$ \hat{\rho}_{\text{BS}} = \frac{1}{4\pi} \int d\Omega |N : \Omega\rangle\langle N : \Omega|. \qquad (31) $$

At zero temperature and zero field, there is no preferred direction for the vector $\Omega$ 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 $\langle N : \Omega' | N : \Omega \rangle$ between two spin-nematic states vanishes very quickly with the distance $|\Omega - \Omega'|$. This allows one to use the approximation $\langle N : \Omega' | N : \Omega \rangle \approx 4\pi\delta(\Omega - \Omega') + \mathcal{O}(1/N)$, which leads to

$$ \langle \hat{O}_k \rangle_{\text{BS}} = \text{Tr}[\hat{\rho}_{\text{BS}} \hat{O}_k] \approx \langle \hat{O}_k \rangle_{\text{singlet}} + \mathcal{O}(1/N). \qquad (32) $$

This result can be written as a general statement concerning the average values of a few-body observables with $k \ll N$ [2]: to the 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$.

It is worth noting the difference between individual states and the ensemble. The moments of $N_0$ in the state $|N : \Omega\rangle$ are given by

$$ \langle N : \Omega | \hat{N}_0 | N : \Omega \rangle = Nu^2, $$

$$ \langle N : \Omega | \hat{N}_0^2 | N : \Omega \rangle = N(N-1)u^4 + Nu^2, $$

where $u = \cos(\theta)$. The variance of $N_0$ for a system prepared in a single spin-nematic state, $Nu^2(1 - u^2)$, is thus Poissonian, as expected for a regular condensate. On the other hand, computing the ensemble averages over $\hat{\rho}_{\text{BS}}$ gives

$$ \langle N_0 \rangle = N \int_0^1 du \, u^2 = \frac{N}{3}, $$

$$ \langle N_0^2 \rangle = \int_0^1 du (N(N-1)u^4 + Nu^2) = \frac{3N^2 + 2N}{15}, $$

$$ \Delta N_0^2 = \frac{4N^2 + 6N}{45}. $$
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 $\propto N$. This is in agreement with the general statement made above.

## 5.2. 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 : \Omega\rangle$. To the leading order in $1/N$, these states have zero interaction energy$^7$ and a mean QZ energy given by $-Nq \cos^2(\theta)$. In the spirit of the mean-field approximation, we replace the Boltzmann factor by its mean value and write the density matrix as

$$ \hat{\rho}_{BS} \approx \frac{1}{Z} \int d\Omega |N : \Omega\rangle\langle N : \Omega| e^{N\beta q \Omega_z^2} \quad (33) $$

with $\beta = 1/k_B T$. The partition function can then be expressed as

$$ Z = \int_{0}^{2\pi} d\phi \int_{0}^{\pi} \sin(\theta) d\theta \, e^{N\beta q \cos^2(\theta)} = 2\pi F_{-1/2}(N\beta q). \quad (34) $$

Here, we introduced the family of functions

$$ F_{\alpha}(y) = \int_{0}^{1} x^{\alpha} e^{y \cdot x} dx, \quad (35) $$

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 the leading order in $N$ as

$$ \langle n_0^m \rangle = \frac{F_{m-1/2}(N\beta q)}{F_{-1/2}(N\beta q)}. \quad (36) $$

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 the leading order in $N$, the moments of $N_0$ obey a universal curve depending only on $Nq/k_B T$ and not on $q/U_s$ or $T/U_s$ separately.

From the properties of the functions $F_\alpha$, we recover the results established in the previous section. When $x \to 0$, one finds $F_\alpha(x) \sim 1/(\alpha+1)$ and $\langle n_0^m \rangle \sim 1/(2m+1)$. By using this result we recover for $q=0$ the previous results, i.e. $\langle n_0 \rangle = 1/3$ and $\Delta n_0^2 = 4/45$. When $x \to \infty$, $F_\alpha(x) \sim e^x/x \times [1-\alpha/x+\alpha(\alpha-1)/x^2]$. This leads to the asymptotic behavior $\langle n_0^m \rangle \sim 1-m/(N\beta q)+m(m-3/2)/(N\beta q)^2+\dots$ when $N\beta q \gg 1$, which reproduces the Bogoliubov results (23), (24) for $q \ll 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

<footnote>$^7$ Explicitly, one has $\langle N : \Omega | \hat{S}^2 | N : \Omega \rangle = N(1+\cos^2(\theta))$, so that the interaction energy of the state $|N : \Omega\rangle$ is given by $U_s \cos^2(\theta) \sim \mathcal{O}(1)$ compared to the QZ energy $\sim \mathcal{O}(N)$. The same argument applies to the off-diagonal matrix elements $\langle N : \Omega' | \hat{S}^2 | N : \Omega \rangle$.</footnote>

<footnote>$^8$ For $q \gg U_s$, the depletion and standard deviation of $\hat{n}_0$ are of the order $1/N$, of the same order as the error made by using the broken symmetry approach.</footnote>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Exact diagonalization (red solid line) versus broken symmetry approach (black dots) for $N = 1000$, $k_B T/U_s = 10$.</figcaption></figure>

low temperatures: although these states are not true eigenstates of the Hamiltonian (1), the action of $\hat{H}$ yields off-diagonal matrix elements scaling as $1/N$ [31], and thus vanishing in the thermodynamic limit. At high temperatures ($k_B T \sim N U_s$), where all the high-energy states are populated the broken-symmetry ansatz is no longer adequate.

## 6. 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 the 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 the 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 paper, we focused on the 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$ [32]. This problem, which can be linked to the more general question of thermalization of closed quantum systems [33], provides an interesting direction for future work.

## Acknowledgments

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).
## Appendix A. Total spin eigenstates

The general expression of the states $|N, S, M\rangle$ in the Fock basis is

$$|N, S, M\rangle = \frac{1}{\sqrt{N}} (\hat{S}^{(-)})^P (\hat{A}^{\dagger})^Q (\hat{a}_{+1}^{\dagger})^S |\text{vac}\rangle. \quad (A.1)$$

Here, $P = S - M$, $2Q = N - S$, $\hat{S}_-$ is the spin lowering operator and $\hat{A}^\dagger = \hat{a}_0^\dagger - 2\hat{a}_{-1}^\dagger \hat{a}_{+1}^\dagger$ is the singlet creation operator. The two operators commute. The normalization constant reads

$$\mathcal{N} = \frac{S!(N-S)!!(N+S+1)!!(S-M)!(2S)!}{(2S+1)!!(S+M)!}, \quad (A.2)$$

where !! indicates a double factorial.

The action of $\hat{a}_0$ on the angular momentum eigenstates is

$$\hat{a}_0|N, S, M\rangle = \sqrt{A_-(N, S, M)}|N-1, S-1, M\rangle + \sqrt{A_+(N, S, M)}|N-1, S+1, M\rangle, \quad (A.3)$$

where $\hat{a}_0$ is the annihilation operator of a boson in the Zeeman state $m=0$, and where the coefficients $A_\pm$ are given by

$$A_-(N, S, M) = \frac{(S^2 - M^2)(N + S + 1)}{(2S - 1)(2S + 1)}, \quad (A.4)$$

$$A_+(N, S, M) = \frac{((S+1)^2 - M^2)(N-S)}{(2S+1)(2S+3)}. \quad (A.5)$$

The non-zero matrix elements of $\hat{N}_0$ are

$$\langle S | \hat{N}_0 | S \rangle = (A_+(N, S, M) + A_-(N, S, M)), \quad (A.6)$$

$$\langle S + 2 | \hat{N}_0 | S \rangle = \sqrt{A_-(N, S+2, M) A_+(N, S, M)}, \quad (A.7)$$

$$\langle S - 2 | \hat{N}_0 | S \rangle = \sqrt{A_+(N, S - 2, M) A_-(N, S, M)}, \quad (A.8)$$

where we abbreviated the notation for the state $|N, S, M\rangle$ as $|S\rangle$ to simplify the notation. We then obtain the matrix elements of $\hat{H}_0$ in the $|N, S, M\rangle$ basis as

$$h_{S,S}^{M} = \frac{U_S}{2N} S(S+1) - q \langle S | \hat{N}_0 | S \rangle, \quad (A.9)$$

$$h_{S,S+2}^{M} = -q \langle S+2 | \hat{N}_0 | S \rangle, \quad (A.10)$$

$$h_{S,S-2}^{M} = -q \langle S-2 | \hat{N}_0 | S \rangle. \quad (A.11)$$
## Appendix B. Continuum approximation

We expand the matrix elements $h_{S,S}$, $h_{S,S \pm 2}$ to the first order in $1/S$, $S/N$, $M^2/S^2$, and obtain

$$h_{S,S \pm 2}^{M} \approx \frac{N}{4} \left[ 1 - \frac{M^2}{(S \pm 1)^2} \right] - \frac{1}{8N} \left[ (S \pm 1)^2 - M^2 \right], \qquad (B.1)$$

$$h_{S,S}^{M \neq 0} \approx \frac{N}{2} \left(1 - \frac{M^2}{S^2}\right), \quad h_{S,S}^{M=0} \approx \frac{N}{2}. \qquad (B.2)$$

For $M=0$, we obtain

$$-\frac{Nq}{4} \left[ \left(1 - \frac{(x+\epsilon)^2}{2}\right) c_{S+2} + \left(1 - \frac{(x-\epsilon)^2}{2}\right) c_{S-2} \right] + \frac{NU_s}{2} x^2 c_s = \left(E + \frac{Nq}{2}\right) c_s, \qquad (B.3)$$

where we have set $x = S/N$ and $\epsilon = 2/N$. We now take the continuum limit, where $\epsilon \ll 1$ is taken as a discretization step and $c_s$ becomes a continuous function $c(S)$. We write

$$\frac{N^2}{4} (c_{S+2} + c_{S-2}) \approx \Delta c(s) + \frac{N^2}{2} c(s). \qquad (B.4)$$

By substituting in (B.4) and neglecting a term $\propto (qx^2/N)\Delta 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 $\sim 1/\sqrt{N}$ for $q \gg U_s$. For the thermal states, the width is $\sim \sqrt{k_B T / [N(2U_s+q)]}$, which gives the condition $k_B T \ll N(U_s+q)$. Finally, the cross-term $\propto (qx^2/N)\Delta c$ is of the 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.

## References

[1] Pitaevskii L and Stringari S 2003 *Bose Einstein Condensation* (Oxford: Oxford University Press)

[2] Castin Y and Herzog C 2001 *C. R. Acad. Sci. Paris* **2** 419–43

[3] Mueller E J, Ho T L, Ueda M and Baym G 2006 *Phys. Rev. A* **74** 033612

[4] Nozières P and Saint James D 1982 *J. Phys. France* **42** 1133–48

[5] Law C K, Pu H and Bigelow N P 1998 *Phys. Rev. Lett.* **81** 5257

[6] Ho T L and Yip S K 2000 *Phys. Rev. Lett.* **84** 4031

[7] Kuklov A B and Svistunov B V 2002 *Phys. Rev. Lett.* **89** 170403

[8] Ashhab S and Leggett A J 2003 *Phys. Rev. A* **68** 063612

[9] Zhou F 2001 *Phys. Rev. Lett.* **87** 080401

[10] Cui X, Wang Y and Zhou F 2008 *Phys. Rev. A* **78** 050701

[11] Barnett R, Sau J D and Das Sarma S 2010 *Phys. Rev. A* **82** 031602

[12] Barnett R, Hui H Y, Lin C H, Sau J D and Das Sarma S 2011 *Phys. Rev. A* **83** 023613

[13] Lamacraft A 2011 *Phys. Rev. A* **83** 033605

[14] Tasaki H 2013 *Phys. Rev. Lett.* **110** 230402

[15] Foot C 2005 *Atomic Physics* (Oxford: Oxford University Press)

[16] Ho T L 1998 *Phys. Rev. Lett.* **81** 742
[17] Ohmi T and Machida T 1998 *J. Phys. Soc. Japan* **67** 1822

[18] Leib E and Matthis D 1962 *J. Math. Phys.* **3** 749

[19] Kaiser C and Peschel I 1989 *J. Phys. A: Math. Gen.* **22** 4257

[20] Kaplan T A, von der Linden W and Horsch P 1990 *Phys. Rev. B* **42** 4663–9

[21] Bernu B, Lhuillier C and Pierre L 1992 *Phys. Rev. Lett.* **69** 2590–3

[22] Azaria P, Delamotte B and Mouhanna D 1993 *Phys. Rev. Lett.* **70** 2483–6

[23] van Wezel J, van den Brink J and Zaanen J 2005 *Phys. Rev. Lett.* **94** 230401

[24] van Wezel J, Zaanen J and van den Brink J 2006 *Phys. Rev. B* **74** 094430

[25] Black A T, Gomez E, Turner L D, Jung S and Lett P D 2007 *Phys. Rev. Lett.* **99** 070403

[26] Jacob D, Shao L, Corre V, Zibold T, De Sarlo L, Mimoun E, Dalibard J and Gerbier F 2012 *Phys. Rev. A* **86** 061601

[27] Yi S, Müstecaplıoğlu Ö E, Sun C P and You L 2002 *Phys. Rev. A* **66** 011601

[28] Stamper-Kurn D M and Ueda M 2013 *Rev. Mod. Phys.* **85** 1191–244

[29] Hooley C and Quintanilla J 2004 *Phys. Rev. Lett.* **93** 080404

[30] Ott H, de Mirandes E, Ferlaino F, Roati G, Türck V, Modugno G and Inguscio M 2004 *Phys. Rev. Lett.* **93** 120407

[31] Anderson P W 1952 *Phys. Rev.* **86** 694–701

[32] Pu H, Law C K, Raghavan S, Eberly J H and Bigelow N P 1999 *Phys. Rev. A* **60** 1463–70

[33] Polkovnikov A, Sengupta K, Silva A and Vengalattore M 2011 *Rev. Mod. Phys.* **83** 863–83