Integrable model of a p-wave bosonic superfluid

We present an exactly-solvable p-wave pairing model for two bosonic species. The model is solvable in any spatial dimension and shares some commonalities with the p+ip Richardson-Gaudin fermionic model, such as a third order quantum phase transition. However, contrary to the fermionic case, in the bosonic model the transition separates a gapless fragmented singlet pair condensate from a pair Bose superfluid, and the exact eigenstate at the quantum critical point is a pair condensate analogous to the fermionic Moore-Read state.


I. INTRODUCTION
Integrable Richardson-Gaudin (RG) models [1,2] based on the su(2) fermion pair algebra have attracted a lot of attention in recent years. Starting with studies of the metal to superconductor transition in ultrasmall grains [3], where the original Richardson's exact solution of the BCS model [4] was rediscovered, to their generalization to a broad range of phenomena in interacting quantum many-body systems [5,6]. The rational or XXX family of integrable RG models has been extensively studied, and includes the constant pairing Hamiltonian (BCS model) [7][8][9], the central spin model [10], generalized Tavis-Cummings models [11], and more recently, open quantum systems [12]. The hyperbolic or XXZ family is much less investigated. The notable p + ip model of p-wave fermionic pairing [13,14,17] is an exception, having the Moore-Read (MR) Pfaffian, proposed for the non-Abelian quantum Hall fluid with filling fraction 5/2 [18,19], as ground state at a given coupling strength.
Another recent finding is a number conserving version of the Kitaev wire which hosts topologically trivial and nontrivial superfluids phases [15]. Interestingly, its repulsive version in the strong coupling limit has been shown to be related to the quantum Hall Hamiltonian projected onto the lowest Landau level subspace [20].
Contrary to the fermionic case, su(1,1) bosonic RG models are unexplored territory. Richardson introduced the bosonic constant pairing Hamiltonian [21], later generalized to study condensate fragmentation for repulsive pairing interactions [22], and the transition from spherical to γ-unstable nuclei in the nuclear interacting boson model [23,24]. The hyperbolic su(1,1) RG, proposed in [2], has only been employed to demonstrate the integrability of the celebrated Lipkin-Meshkov-Glick model in the Schwinger boson representation [6,25,26]. In this Letter we derive an integrable two bosonic species p-wave pairing Hamiltonian, and study its quantum phase diagram. We are motivated by the recent experimental observation of broad p-wave resonances in ultracold 85 Rb and 87 Rb atomic mixtures [27,28] that could lead to stable thermodynamics phases dominated by p-wave attractive interactions. Mean-field studies based on a twochannel model predict three phases [29]: a) an atomic  6) , in the space g = GL and ρ = N/(2L), where N is the number of bosons and 2L the size of the lattice. The gray area indicates the singlet pair fragmented BEC phase, separated from the pair Bose superfluid phase (PBS) by a critical line gc. In the superfluid phase the Volovik line, depicted by the lower dashed line, delineates the phase region where the minimum of the quasi-particle energy is at momentum k min = 0, while to the right of the upper dashed line it is at k min = π. Between these two lines 0 < k min < π. All pairons diverge to infinity at g∞. Symbols indicate the couplings used below at density ρ = 0.2. Horizontal dashed line is ρ = 0.299433. The inset shows the behavior of the pairing gap ∆ versus coupling strength g at density ρ = 0.2.
Bose-Einstein condensate (BEC) for large negative detuning, b) a molecular BEC for large positive detuning, and c) an atomic-molecular BEC for intermediate detuning. Quantum fluctuations may stabilize the atomicmolecular phase for certain densities giving rise to the formation of polar superfluid droplets [30]. Our exactlysolvable attractive one-channel p-wave model displays two phases (see Fig. 1), a gapless fragmented BEC of singlet pairs, where each of the species condenses into the lowest finite momentum (grey area), and a gapped pair Bose superfluid (white area).

II. HYPERBOLIC su(1,1) RG INTEGRALS OF MOTION
The hyperbolic su(1,1) model for two bosonic species a and b in momentum k space is based on the interspecies pair operators In order to satisfy the su(1,1) algebra [K − k,Q , K + k ,Q ] = 2δ k,k K z k,Q and [K z k,Q , K ± k ,Q ] = ±δ k,k K ± k,Q , and to avoid double counting, we restrict momenta k and Q to have the component along one of the dimensions, for instance k x , larger than zero, k x > 0 and Q x > 0. As we will see below this does not restrict the k values in the Brillouin zone. The operator K + k,Q , that creates a two-species pair with center-of-mass momentum Q, is antisymmetric under the exchange of species. If we interpret both species as the two components of a pseudo-spin 1/2, the pair operator K + k,Q creates a singlet state. The pseudo-spin 1/2 bosons define an independent and commuting su(2) spin algebra generated by S z k, Although we will focus on the Q = 0 case, these commuting algebras can be exploited to describe Larkin-Ovchinnikov-Fulde-Ferrell-type phases and/or mass imbalance two-component cold atom gases as described in [9] for fermionic systems.
In terms of the su(1,1) generators (1), the hyperbolic integrals of motion for Q = 0 are [2, 6] where η k are arbitrary odd functions of k. The sum k > 0 means that the component k x should be positive. For a fixed number of bosons N = 2M + ν, where M is the number of singlet boson pairs and ν the total number of unpaired bosons, the eigenvalues of the integrals of motion are where d k = ν k /2 + 1, ν k is the seniority quantum number (number of unpaired bosons) of level k, and ν = k>0 ν k . The spectral parameters e α , so-called pairons, are roots of the Richardson equations (α = 1, . . . , M ) Each independent solution of the Richardson equations (3) defines a common eigenstate of the integrals of motion (2): where the state |ν , with ν unpaired bosons, satisfieŝ K − k |ν = 0 for all k, andK z k |ν = d k |ν . By combining the integrals of motion R k with the Hellmann-Feynman theorem [14], the occupation probabilities can be obtained from the expectation value where the pairon derivatives can be obtained from the derivatives of Eq. (3) leading to a linear set of equations. For ease of presentation we consider next a onedimensional version of the p-wave model. It is straightforward to extend our model to higher dimensions as has been done in the fermionic case [13,14].

III. THE p-WAVE BOSE HAMILTONIAN
The p-wave pairing Bose Hamiltonian we want to study is given by where η k = sin(k/2) and η 2 k = (1 − cos k)/2. Assuming antiperiodic boundary conditions the allowed k values are k = ±π/2L, ±3π/2L, . . . , ±(2πL − π)/2L, with 2L the size of the chain and L the number of su(1,1) copies. We have chosen antiperiodic boundary conditions to explicitly exclude the k = 0 state. This state cannot support singlet pairs and, therefore, it will be excluded from the dynamics of p-wave pair scattering. This model Hamiltonian, which written in terms of the su(1, 1) generators is  Fig. 1. Coupling strengths G are indicated by numbers close to the respective symbols. Horizontal lines depict η 2 k levels. See [31] for an animation of the pairons evolution as a function of G for 0 < G < G∞.
can be derived from the hyperbolic su(1,1) RG integrals of motion (2), by using the linear combination Our p-wave Hamiltonian (7) has an explicit U(1) symmetry, i.e., conservation of the total number of bosons, and a pseudospin invariance that basically preserves the polarization S z = k S z k , that is, the difference between the number of bosonic species. Here we will focus on an unpolarized mixture of atoms characterized by S z = 0, although the polarized case (S z = 0) is contained in our exact solution. For instance, an excess of a atoms manifests through the seniorities ν k specifying the k states occupied by the unpaired a atoms.
Eigenvalues of (7) can be determined from the integrals of motion, using the same linear combination, which, after using Eq. (3), gives Let us analyze next the way the ground state evolves as a function of coupling strength G ≥ 0 (Fig. 2). Each independent solution of the Richardson equations (3) provides a set of M pairons that define both the energy eigenvalue (8) and the corresponding eigenstate (4). The ground state (with ν = 0) for weak coupling G has the pairons distributed in the real interval between zero and the minimum η 2 π/2L = sin 2 (π/4L). At G n = 2/(2L + 2M − n − 1), n pairons collapse to zero. In between collapses, the n pairons expand as complex conjugate pairs forming an arc in the complex plane around zero. The whole set of M pairons collapses to zero at the critical point where the exact (non-normalized) ground state becomes a condensate of singlet pairs which is algebraically analogous to the MR state of the p + ip fermionic model [13,14], and, therefore, we will call it Bose Moore-Read (BMR) state. Naively, in an extended system, one would expect that the ground state of the BEC phase, 0 ≤ G ≤ G c , corresponds to a zero-momentum condensate for each species since, as we will see, the quasi-particle gap ∆ vanishes. This state has maximum spin S = M . For mesoscopic systems, it has been shown that the correct ground state at weak coupling is a fragmented singlet pair BEC [32,33], which in momentum space becomes with k min = π/2L for the antiperiodic chain. Note that in this phase, the exact ground state has a mixture of complex pairons close to zero and real pairons in the interval [0, η 2 π/2L ]. For large L the pairons will cluster around zero and the exact ground state (4) will tend to the BMR state (9) which is representative of the whole phase. The BMR state is controlled by k min , and therefore it converges to the singlet pair condensate in the large L limit. Interestingly, in the thermodynamic limit the states (11) and (10), as well as condensates with other spin quantum numbers S, become degenerate. A weak repulsive interaction may destabilize those degenerate spin states against the singlet pair condensate [33].
For G > G c the pairons distribute along an arc that expands in the complex plane as G increases (Fig. 2). At the absolute value of all pairons diverges to infinity. This divergence does not affect the energy since imaginary parts cancel out pairwise in (8) and the real parts combine to give E = 2GM k>0 η 2 k . Infinite pairon energies have been observed previously in fermionic hyperbolic models [15] and they were related to a duality associated to the particle-hole symmetry [16]. At this point the exact ground state can be expressed as a different pair condensate In turn, in the fermionic case we find that this state appears as the highest energy eigenstate in the repulsive pairing region. In Fig. 1, at density ρ = 0.2, we show five distinct symbols covering all distinct areas of the phase diagram, at couplings g = 0.5, 1.2, 1.8, 5.0, 6.8, with g = GL. Figure 2 displays pairons of a finite-sized system with M = 10 and L = 50, for these same five values. As discussed above, the first point with G < G c has 10 pairons distributed in the real positive axis below η 2 π/2L (see the right inset). After the pairons collapse to zero at G c , they form an arc in the complex plane that expands for increasing values of G. This is the case for the remaining four couplings that lay in between G c and G ∞ , two of them can be seen in the left inset while the other two in the central figure.

IV. QUANTUM PHASE DIAGRAM
The thermodynamic limit is obtained in the limit of N, L → ∞ with constant density ρ = N/(2L) and rescaled interaction strength g = GL. In this limit, the Richardson equations (3) transform into the boson gap and number equations [8,14] with quasi-boson energies E k and occupation probabili- where µ is the chemical potential and ∆ the gap. Though E k in (14) may, in principle, be complex, we have numerically verified that in the large attractive g limit, Eqs. (13) have solutions µ ≈ −γ 1 g and ∆ ≈ γ 2 g, with γ 1,2 positive constants satisfying 4γ 2 2 < γ 2 1 . This latter condition guarantees that the quasi-boson energies, given by E k ≈ gγ 1 1 − (4γ 2 2 /γ 2 1 )η 2 k are always real, even in the limit g → ∞. The ground state energy density E ≡ E/L for a given density ρ in the thermodynamic limit is given by The critical coupling of the exact solution in the finitesize case, becomes g c = lim L,N →∞ [G c L] = 2/(2 + ρ) in the thermodynamic limit. The gap ∆ is zero at weak pairing up to the critical value g c . The inset of Fig. 1 shows the behavior of the gap for ρ = 0.2. It increases monotonically for g > g c . In the same thermodynamic limit, the coupling where all pairons diverge becomes g ∞ = lim L,N →∞ [G ∞ L] = 2/ρ (Fig. 1).
We are interested in establishing the nature of the nonanalyticities of E at the critical point. It turns out that E = 0, for 0 < g < g c and is non-analytic at g = g c with a third order phase transition to a pair superfluid phase [14,34]. Close to g − g c ≈ 0 + , it behaves as whereg = (g − g c )/g c (See Appendix A). Interestingly, the behavior of E close to g c depends on ρ only through its critical value g c . The first and second-order derivatives at the critical point are zero, while the third-order derivative is ∂ 3 g E g−gc→0 + = −2π 2 /g 6 c , signaling a discontinuity of third order. This is illustrated in Fig. 3 for ρ = 0.2 where, moreover, E is compared with the exact energy density for M = 10 and L = 50.

V. NATURE OF EXCITATIONS
In Fig. 4 we show the quasi-boson energies for the five values of g indicated in Fig. 1. The quasi-boson energies change from E k = sin 2 (k/2) in the gapless pair condensate phase (g = 0.5), to a complex dispersion in the pair Bose superfluid phase. For µ + 2∆ 2 ≤ 0, E k is a monotonous increasing function with minimum at k min = 0 and energy E k min = |µ| (g = 1.8). The previous condition is fulfilled in the superfluid phase only for small densities ρ < 0.299433 in a finite coupling interval. The region is indicated by the area with diagonal lines in Fig. 1. The boundary of this region, the so called Volovik line [13] defined by a superfluid with the minimum quasi-boson energy at k = 0, is given by µ + 2∆ 2 = 0. For 0 < µ + 2∆ 2 < 1, E k has a minimum at k min = 2 arcsin( µ + 2∆ 2 ), satisfying 0 < k min < π (g = 1.1, 5.0). The region of the phase diagram where E k has this dispersion is indicated by the white area in Fig.  1. The previous condition is fulfilled for any density, and gives the form of the quasi-boson dispersion immediately after the quantum phase transition. For µ + 2∆ 2 ≥ 1 (area with horizontal lines in Fig. 1), the quasi-boson dispersion is a monotonous decreasing function with minimum at k min = π (g = 6.8).
The occupation probabilities in momentum space are displayed in the lower inset of Fig. 4 for the five values of g indicated in Fig. 1. Continuous lines are the thermodynamic limit solution and symbols correspond to the exact solution for the finite-size case calculated using Eq. (5). For g = 0.5 the system is condensed in k min resulting in a delta distribution in the thermodynamic limit. At g c , in that limit, the macroscopic occupation at k min → 0 jumps to zero and then the maximum of the distribution moves to finite k values. This jump in the k = 0 momentum state resembles the one observed in the p + ip and RG Kitaev models [14,15] and the s-d RG model of Ref. [35]. In the fermionic case, this fact has been linked to a topological phase transition [14,15]. For g = 1.1 and 1.8 the profiles broaden and maxima get displaced to larger values of k. Finally, for g = 5.0 and 6.8 the profiles are inverted with a maximum occupation at k = π.

VI. OUTLOOK
We introduced an exactly-solvable two species p-wave bosonic model and established its quantum phase diagram in the attractive sector. Only the case of a balanced mixture with equal masses (m a = m b ) and zero center-of-mass momentum Q has been studied in depth. Imbalanced binary mixtures (ν = 0, m a = m b ) and finite Q pairs are contained within our exactly solvable model. The exact, finite and thermodynamic limit, treatments of the p-wave pairing Bose Hamiltonian (7), although seemingly similar, have profound physical differences to its fermionic counterpart [13,14,16,17] despite the fact that both cases share a third-order quantum phase transition. In the fermionic case the latter separates two gapped superfluid phases and has a topological character [15]. In the bosonic case one of the phases is gapless and displays a fragmented BEC condensate with macroscopic occupations of both species in the lowest finite momentum pair states (−k, k), while the other is a gapped pair Bose superfluid (PBS). Moreover, while for fermions the critical coupling takes place at the Read-Green point, with one pairon at zero energy and the other M − 1 pairons with real and negative energies, for bosons the phase transition takes place at the equivalent of the fermionic Moore-Read point with all pairons collapsing to zero energy. It is at this critical point that the exact bosonic ground state is a pair condensate with amplitudes fixed by the single particle energies.
Motivated by a theoretical prediction [36], recent experiments unveiled a new type of ultradilute quantum liquid in ultracold bosonic systems. Apparently, there is no unique mechanism leading to such a liquid state since it has been observed in single-species dipolar systems [37] and Bose (potassium) mixtures [38,39]. Can one obtain a quantum liquid phase in p-wave Bose systems? This question has been recently addressed in [30], and answered in the affirmative for a particular model. Our PBS represent a (fixed-point) number-conserving candidate for such quantum liquid phase. The pairing interaction in (7) may thus provide a new effective mechanism for its emergence. Although the superfluid gap protects that state from expansion in finite geometries, further studies in trapped potentials are required to identify a possible self-bound quantum liquid droplet. On the experimental side, it is crucial to have a precise understanding of the spectrum of excitations to compare to our theoretical predictions.