Spin fragmentation of Bose–Einstein condensates with antiferromagnetic interactions

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_BT ≫ 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)

$$\begin{equation} \fl -J(x+\epsilon) c(x+\epsilon)-J(x-\epsilon)c(x-\epsilon) +\frac{NU_{\rm s}}{2} x^2 c(x)=\left(E +\frac{Nq}{2}\right) c(x),\qquad \end{equation} \tag{ 4 }$$

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) = Nq (1 − x²/2)/4 and a harmonic potential strength NU_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

$$\begin{equation} -\frac{q}{N} c''(x)+\frac{N}{4}\left(q+2 U_{\rm s}\right) x^2 c(x)=\left(E+Nq\right)c(x). \end{equation} \tag{ 5 }$$

The boundary condition c(0) = 0 selects the eigenstates of the standard harmonic oscillator with odd parity. The mass m and oscillation frequency ω of the fictitious oscillator are found from ℏ²/2m ≡ q/N and mω² ≡ N(q + 2U_s)/2. The oscillator frequency is thus

$$\begin{equation} \hbar \omega = \sqrt{q(q+2 U_{\rm s})}. \end{equation} \tag{ 6 }$$

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

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) is given by

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

with the quantum harmonic oscillator size

$$\begin{equation} \sigma =\sqrt{\frac{\hbar}{m \omega}} = \sqrt{\frac{2}{N}} \left( \frac{q}{q+2U_{\rm s}} \right)^{1/4}. \end{equation} \tag{ 8 }$$

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

$$\begin{eqnarray} &&q \gg\frac{U_{\rm s}}{N^2}. \end{eqnarray} \tag{ 9 }$$

For q < U_s/N², 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², 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 ${\sim} N\sigma \approx \sqrt {2N}$ total spin eigenstates. In this regime, we can compute the moments of N₀ 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

$$\begin{equation} N-\langle N_0\rangle = \frac{U_{\rm s}+q}{2\sqrt{q(q+2U_{\rm s})}}, \end{equation} \tag{ 10 }$$

$$\begin{equation} \Delta N_0^2 = \frac{U_{\rm s}^2}{2q(q+2U_{\rm s})}. \end{equation} \tag{ 11 }$$

For U_s/N² ≪ q ≪ U_s, the depletion N − 〈N₀〉 and variance ΔN²₀ are larger than unity but small compared to N,N², respectively, while for q ≫ U_s, they become less than one particle: in the latter case, the ground state approaches the Fock state $( \hat {a}_0^\dagger )^N \vert {\mathrm { vac}} \rangle$ expected from the mean field theory. We indicate both regimes as 'BEC m = 0' in figure 1, without marking the distinction.

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)/4 and a harmonic potential strength κ = NU_s. Let us examine two limiting cases. For q = 0 (no hopping), the energy eigenstates coincide with the 'position' eigenstates with energy E(S) ≈ (U_s/2N)S² for S ≫ 1. Conversely, when U_s = 0 the energy eigenstates are delocalized states, which form an allowed energy band of width  ∼ 4J ∼ 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, κ ≠ 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 ≪ 4J = Nq (the requirement q ≫ U_s/N² found before still holds); and
• 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)\approx \sqrt {2 E/NU_{\mathrm { 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)|² 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, _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_sn²/2N with n integer. Both expressions agree well with the numerical result in their respective domains of validity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

$$\begin{equation} \langle \skew3\hat{N}_0 \rangle_{SM} \approx (N^2-S^2)(S^2-M^2)/8, \end{equation} \tag{ 12 }$$

$$\begin{equation} \langle \skew3\hat{N}_0^2 \rangle_{SM} \approx (3N^2-S^2)(S^2-M^2)^2/8, \end{equation} \tag{ 13 }$$

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

where $\langle \hat {N}_0^p \rangle _{SM} = \langle N,S,M \mid N_{0}^{\mathrm {p}} \mid N,S,M \rangle$. 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₀〉_T and $(\Delta n_{0}^2)_{T}$ in the canonical ensemble. The average population in m = 0 is given by

$$\begin{equation} \langle N_{0}\rangle_{T} = \frac{1}{Z} \sum_{S,M}\mathrm{e}^{-\beta' S(S+1)} \langle N_{0}\rangle_{SM}. \end{equation} \tag{ 15 }$$

The second moment 〈N²₀〉_T and the variance $(\Delta N_{0}^2)_T$ are given by similar expressions. Here, Z is the partition function and β' = U_s/2Nk_BT. Assuming that the temperature is large compared to the level spacing (k_BT ≫ 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〉 ∼ (Nk_BT/U_s)^1/2, shows that the condition 1 ≪ S ≪ N corresponds to the boundaries

$$\begin{equation} \frac{U_{\rm s}}{N} \ll k_{\rm B}T \ll N U_{\rm s}. \end{equation} \tag{ 16 }$$

In this regime, we find

$$\begin{equation} \langle N_{0}\rangle_{T} \approx \frac{N}{3}, \end{equation} \tag{ 17 }$$

$$\begin{equation} ( \Delta N_{0}^2 )_{T } =\langle N_{0}^2\rangle_{T}-\langle N_{0}\rangle_{T}^2 \approx \frac{4 N^2}{45},\quad k_{\rm B} T \ll NU_{\rm s}. \end{equation} \tag{ 18 }$$

We note that to the leading order in 1/N, the moments of N₀ are identical for those found in the singlet state.

The second regime arises when the temperature becomes very large (k_BT/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 ⩽ 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

$$\begin{equation} (\Delta N_{0}^2)_{T\gg NU_{\rm s}} \approx N^2/18. \end{equation} \tag{ 19 }$$

To summarize (see figure 1), for q = 0 we always find large depletion and super-Poissonian fluctuations (ΔN²₀ ∼ 〈N₀〉²). 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_BT ≪ NU_s, and changes to $1/3\sqrt {2}\approx 0.236$ for very large temperatures k_BT > NU_s where all the states are occupied with equal probability.

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 [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}_{0}$, and expresses the fluctuations ${\hat a}_{\pm 1}$ in terms of new operators $\hat {\alpha }^\pm$,

$$\begin{equation} \skew3\hat{\alpha}_{\pm} = u {\hat a}_{\pm 1}-v {\hat a}_{\mp 1}^\dagger. \end{equation} \tag{ 20 }$$

Here, the Bogoliubov amplitudes u,v defined by

$$\begin{equation} u\pm v = \left( \frac{q}{2U_{\rm s}+q}\right)^{\pm 1/4} \end{equation} \tag{ 21 }$$

are chosen to put the Hamiltonian in diagonal form

$$\begin{equation} H_{\rm Bogo} = \sum_{\mu=\pm} \hbar \omega \left( \skew3\hat{\alpha}_{\mu}^\dagger \skew3\hat{\alpha}_{\mu} +\frac{ 1}{2}\right)-(g+q). \end{equation} \tag{ 22 }$$

The energy ℏω 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).⁶

In the Bogoliubov approximation, the moments of N₀ can be obtained analytically. The quantum (T = 0) depletion of N₀ is smaller than one atom. The thermal part of the depletion and variance of n₀ = N₀/N read for k_BT ≫ ℏω:

$$\begin{equation} 1-\langle n_0\rangle = \frac{2(U_{\rm s}+q)}{q+2U_{\rm s}} \frac{k_{\rm B} T}{Nq}, \end{equation} \tag{ 23 }$$

$$\begin{equation} \Delta n_0^2 = \frac{2\left[(U_{\rm s}+q)^2+U_{\rm s}^2\right]}{(q+2U_{\rm s})^2} \left(\frac{k_{\rm B} T}{Nq}\right)^2. \end{equation} \tag{ 24 }$$

The prefactors take values of order unity, and both the depletion 1 − 〈n₀〉 and standard deviation Δn₀ scale as k_BT/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 (〈n₀〉 = 1), or in other words for temperatures

$$\begin{equation} k_{\rm B} T \ll Nq. \end{equation} \tag{ 25 }$$

We compare in figure 4 the predictions for the moments of N₀ 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.

Zoom In Zoom Out Reset image size

When Nq/k_BT ≪ 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_BT ≫ 1, thermal states mostly populate states with E ∼ qN, i.e. 'delocalized' states within the low-energy 'Bloch band' of width ∼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 ≪ k_BT/N or a thermal state of Bogoliubov-like excitations for q ≫ k_BT/N.

We note to conclude this section that in the regime U_s/N ≪ k_BT ≪ NU_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_BT/U_s and Nq/U_s. We found that the physical quantities 〈N₀〉, $\left (\Delta N_0 \right )$ depend only on their ratio Nq/k_BT, 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.

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〉 in this basis, one finds [2, 3, 12]

$$\begin{equation} \vert N, S, M\rangle \propto \int \mathrm{d}{\boldsymbol{\Omega}} ~Y_{S,M} ({\boldsymbol{\Omega}}) \vert N: {\boldsymbol{\Omega}} \rangle, \end{equation} \tag{ 28 }$$

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 $\langle \hat {O}_k \rangle _{\mathrm { singlet}} =\langle N,0,0 \vert \hat {O}_k \vert N,0,0 \rangle$ of a k-body operator $\hat {O}_k$,

$$\begin{equation} \langle \skew3\hat{O}_k \rangle_{\rm singlet} = \frac{1}{4\pi} \int \mathrm{d}{\boldsymbol{\Omega}} \,\mathrm{d}{\boldsymbol{\Omega}'} \langle N: {\boldsymbol{\Omega}'} \vert \skew3\hat{O}_k \vert N: {\boldsymbol{\Omega}}\rangle. \end{equation} \tag{ 29 }$$

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

$$\begin{equation} \langle \skew3\hat{O}_k \rangle_{\rm singlet} \approx \frac{1}{4\pi} \int \mathrm{d}{\boldsymbol{\Omega}} \langle N: {\boldsymbol{\Omega}} \vert \skew3\hat{O}_k \vert N: {\boldsymbol{\Omega}}\rangle. \end{equation} \tag{ 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]

$$\begin{equation} \skew3\hat{\rho}_{\rm BS} = \frac{1}{4\pi} \int \mathrm{d}{\boldsymbol{\Omega}} \vert N: {\boldsymbol{\Omega}} \rangle \langle N: {\boldsymbol{\Omega}} \vert. \end{equation} \tag{ 31 }$$

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 $\langle N: {\boldsymbol {\Omega }'} \vert N: {\boldsymbol {\Omega }}\rangle \approx 4\pi\delta ({\boldsymbol {\Omega }-\boldsymbol {\Omega }'})+ {\cal O}(1/N)$, which leads to

$$\begin{equation} \langle \skew3\hat{O}_k \rangle_{\rm BS} = {\rm Tr}[ \skew3\hat{\rho}_{\rm BS} \skew3\hat{O}_k] \approx \langle \skew3\hat{O}_k \rangle_{\rm singlet} + {\cal O}(1/N). \end{equation} \tag{ 32 }$$

This result can be written as a general statement concerning the average values of a few-body observables with k ≪ 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₀ in the state |N:Ω〉 are given by

$$\begin{eqnarray*} \langle N:{\boldsymbol{\Omega}}\vert \skew3\hat{N}_0 \vert N:{\boldsymbol{\Omega}} \rangle &= &N u^2,\\ \langle N:{\boldsymbol{\Omega}}\vert \skew3\hat{N}_0^2 \vert N:{\boldsymbol{\Omega}} \rangle &=& N(N-1)u^4+N u^2, \end{eqnarray*}$$

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

$$\begin{eqnarray*} \langle N_0 \rangle &=&N \int_0^1 \mathrm{d}u \,u^2 = \frac{N}{3},\\ \langle N_0^2 \rangle &=& \int_0^1 \mathrm{d}u \left( N(N-1) u^4 +N u^2 \right)=\frac{3N^2+2N}{15},\nonumber \\ \Delta N_0^2&=&\frac{4N^2+6N}{45}.\end{eqnarray*}$$

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.

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 the leading order in 1/N, these states have zero interaction energy⁷ and a mean QZ energy given by −Nq cos²(θ). In the spirit of the mean-field approximation, we replace the Boltzmann factor by its mean value and write the density matrix as

$$\begin{equation} \skew3\hat{\rho}_{\rm BS} \approx \frac{1}{{\cal Z} } \int \mathrm{d}{\boldsymbol{\Omega}} \vert N: {\boldsymbol{\Omega}} \rangle \langle N: {\boldsymbol{\Omega}} \vert\mathrm{ e}^{N \beta q \Omega_z^2} \end{equation} \tag{ 33 }$$

with β = 1/k_BT. The partition function can then be expressed as

$$\begin{equation} {\cal Z} = \int_{0}^{2\pi}\mathrm{d}\phi\int_{0}^{\pi} \sin(\theta)\mathrm{d}\theta~ \mathrm{e}^{N \beta q \cos^2(\theta)}=2\pi F_{-1/2}\left(N\beta q\right). \end{equation} \tag{ 34 }$$

Here, we introduced the family of functions

$$\begin{equation} F_\alpha(y)=\int_0^1 x^{\alpha} \,\mathrm{e}^{y x}\, \mathrm{d}x, \end{equation} \tag{ 35 }$$

which are related to the lower incomplete gamma functions. In a similar way, we can compute the moments of n₀ = N₀/N to the leading order in N as

$$\begin{equation} \langle n_0^m \rangle = \frac{F_{m-1/2}\left( N\beta q\right)}{F_{-1/2} \left( N\beta q\right)}. \end{equation} \tag{ 36 }$$

From this result, one can easily deduce the average and variance of n₀. This calculation provides an explicit proof of the numerical evidence that, to the leading order in N, the moments of N₀ obey a universal curve depending only on Nq/k_BT 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, one finds F_α(x) ∼ 1/(α + 1) and 〈n^m₀〉 ∼ 1/(2m + 1). By using this result we recover for q = 0 the previous results, i.e. 〈n₀〉 = 1/3 and Δn²₀ = 4/45. When x → ∞, $F_\alpha (x) \sim \mathrm {e}^x/x\times [1-\alpha /x+\alpha (\alpha -1)/x^2]$. This leads to the asymptotic behavior 〈n^m₀〉 ∼ 1 − m/(Nβq) + m(m − 3/2)/(Nβq)² + ···  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 $\hat {H}$ yields off-diagonal matrix elements scaling as 1/N [31], and thus vanishing in the thermodynamic limit. At high temperatures (k_BT ∼ NU_s), where all the high-energy states are populated the broken-symmetry ansatz is no longer adequate.

Zoom In Zoom Out Reset image size