Integrable model of a $p$-wave bosonic superfluid

TL;DR

This paper introduces an exactly solvable $p$-wave pairing model for two bosonic species, revealing a quantum phase transition between a gapless singlet pair condensate and a pair Bose superfluid, with unique critical state properties.

Contribution

It presents a novel exactly solvable bosonic $p$-wave pairing model that shares features with fermionic models but exhibits distinct phase transition characteristics.

Findings

Identifies a third order quantum phase transition in the model.

Describes the critical state as a pair condensate similar to Moore-Read state.

Shows the transition separates a gapless condensate from a superfluid.

Abstract

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.

Equations

$$K_{\mathbb{k},\mathbb{Q}}^{+}$$

$$K_{\mathbb{k},\mathbb{Q}}^z$$

$$R_{\mathbb{k}}$$

$$r_{\mathbb{k}}=d_{\mathbb{k}}\left[1+2\lambda \sum _{{\mathbb{k}}^{\prime }(\ne \mathbb{k})>0}{d}_{{\mathbb{k}}^{\prime }}\frac{{\eta }_{\mathbb{k}}^2+{\eta }_{{\mathbb{k}}^{\prime }}^2}{{\eta }_{\mathbb{k}}^2-{\eta }_{{\mathbb{k}}^{\prime }}^2}-2\lambda \sum _{\alpha =1}^M\frac{e_{\alpha }+{\eta }_{\mathbb{k}}^2}{e_{\alpha }-{\eta }_{\mathbb{k}}^2}\right],$$

$$\sum _{\mathbb{k}>0}\frac{d_{\mathbb{k}}}{{\eta }_{\mathbb{k}}^2-e_{\alpha }}+\sum _{\beta =1(\beta \ne \alpha )}^M\frac{1}{e_{\beta }-e_{\alpha }}+\frac{\tilde{Q}}{e_{\alpha }}=0,$$

$$\tilde{Q}=-\frac{1}{4\lambda }+\frac{M-1+{\sum }_{\mathbb{k}>0}{d}_{\mathbb{k}}}{2}.$$

$$|{\Phi }_{M,\nu }⟩=\prod _{\alpha =1}^M\left(\sum _{\mathbb{k}>0}\frac{{\eta }_{\mathbb{k}}}{{\eta }_{\mathbb{k}}^2-e_{\alpha }}{K}_{\mathbb{k}}^{+}\right)|\nu ⟩,$$

$$⟨{\Phi }_{M,\nu }|K_{\mathbf{k}}^z|{\Phi }_{M,\nu }⟩=d_{\mathbf{k}}\left(1-2{\lambda }^2\sum _{\alpha =1}^M\frac{2{\eta }_{\mathbf{k}}^2}{\left({\eta }_{\mathbf{k}}^2-e_{\alpha }\right)}\frac{\partial e_{\alpha }}{\partial \lambda }\right),$$

$$H$$

$$H=\sum _{k>0}{\eta }_k^2{\hat{N}}_k-G\sum _{k,k^{\prime }>0}{\eta }_k{\eta }_{k^{\prime }}{K}_k^{+}{K}_{k^{\prime }}^{-},$$

$$H=2x\sum _{k>0}{\eta }_k^2{R}_k-2Z-G\sum _{k>0}{d}_k(1-d_k){\eta }_k^2,$$

$$E=\sum _{k>0}{\eta }_k^2{\nu }_k+2\sum _{\alpha =1}^M{e}_{\alpha }.$$

$$G_c=\frac{2}{2L+M-1},$$

$${|{\Phi }_M⟩}_{\mathsf{BMR}}={\left(\sum _{k>0}\frac{1}{{\eta }_k}{K}_k^{+}\right)}^M|0⟩,$$

$$|\Phi ⟩=a_0^{†M}{b}_0^{†M}|0⟩$$

$$|\Phi ⟩=(K_{k_{\mathsf{min}}}^{+}{)}^M|0⟩={\left(b_{k_{\mathsf{min}}}^{†}{a}_{-k_{\mathsf{min}}}^{†}-a_{k_{\mathsf{min}}}^{†}{b}_{-k_{\mathsf{min}}}^{†}\right)}^M|0⟩,$$

$$G_{\infty }=\frac{2}{M-1}$$

$$|{\Phi }_M{⟩}_{G_{\infty }}={\left(\sum _k{\eta }_k{K}_k^{+}\right)}^M|0⟩.$$

$$\frac{\pi }{g}={\int }_0^{\pi }\frac{{\eta }_k^2}{E_k}dk,\rho =\frac{2}{\pi }{\int }_0^{\pi }{v}_k^{2}dk,$$

$$E_k=\sqrt{{\left({\eta }_k^2-\mu \right)}^2-4{\eta }_k^2{\Delta }^2},v_k^2=\frac{1}{2}\left(\frac{{\eta }_k^2-\mu }{E_k}-1\right),$$

$$\mathcal{E}=-\frac{4{\Delta }^2}{g}-1+\frac{2}{\pi }{\int }_0^{\pi }\frac{{\eta }_k^2({\eta }_k^2-\mu )}{E_k}dk.$$

$$\mathcal{E}\approx -\frac{{\pi }^2}{3}{\left(\frac{\tilde{g}}{g}\right)}^2\left(\frac{\tilde{g}}{g}-12\frac{g-1}{g{e}^{2(2-g_c)}}{e}^{\frac{2(g_c-1)}{\tilde{g}}}\right),$$

$$\pi (\rho +1)$$

$$\frac{\pi }{g}$$

$${\int }_0^{1}dx\frac{1}{\sqrt{x-x^2}\sqrt{x^2-2ax+{\mu }^2}}\stackrel{\delta \to 0^{+}}{=}{I}_d$$

$$I_d=\sqrt{\frac{1}{a-a^2}}\left(\log 16+2\log (a-a^2)-\log \delta \right).$$

$${\int }_0^{1}dx\frac{x-a}{\sqrt{x-x^2}\sqrt{x^2-2ax+{\mu }^2}}$$

$$\pi \rho$$

$$\frac{\pi }{g}$$

$$\rho +1$$

$$\mu$$

$$\Gamma$$

$$\mathcal{E}$$

Full text

I Introduction

Integrable Richardson-Gaudin (RG) models Amico2001 ; Dukelsky2001 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 Duk00 , where the original Richardson’s exact solution of the BCS model Richardson1963 was rediscovered, to their generalization to a broad range of phenomena in interacting quantum many-body systems Duk04 ; Ortiz2005 . The rational or XXX family of integrable RG models has been extensively studied, and includes the constant pairing Hamiltonian (BCS model) Rich1966 ; Ortiz2005-2 ; Duk2006 , the central spin model Bortz2010 , generalized Tavis-Cummings models Duke2004 , and more recently, open quantum systems Row2018 . The hyperbolic or XXZ family is much less investigated. The notable $p+ip$ model of $p$-wave fermionic pairing Sierra2009 ; Rombouts2010 ; VanNeck2014 is an exception, having the Moore-Read (MR) Pfaffian, proposed for the non-Abelian quantum Hall fluid with filling fraction 5/2 Moo91 ; Rea00 , 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 non-trivial superfluids phases Ortiz2014 . 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 Ortiz2013 .

Contrary to the fermionic case, $su$(1,1) bosonic RG models are unexplored territory. Richardson introduced the bosonic constant pairing Hamiltonian Rich1967 , later generalized to study condensate fragmentation for repulsive pairing interactions Schuck2001 , and the transition from spherical to $\gamma$-unstable nuclei in the nuclear interacting boson model Pan ; Pittel . The hyperbolic $su$(1,1) RG, proposed in Dukelsky2001 , has only been employed to demonstrate the integrability of the celebrated Lipkin-Meshkov-Glick model in the Schwinger boson representation Ortiz2005 ; Pan2 ; Lerma2013 . 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 85Rb and 87Rb atomic mixtures Papp2008 ; Dong2016 that could lead to stable thermodynamics phases dominated by $p$-wave attractive interactions. Mean-field studies based on a two-channel model predict three phases Rad2009 : a) an atomic 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 atomic-molecular phase for certain densities giving rise to the formation of polar superfluid droplets Li2019 . Our exactly-solvable 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 ${\bf k}$ space is based on the interspecies pair operators

[TABLE]

where $n^{b}_{\bf k}\!=b^{\dagger}_{\mathbb{k}}b^{\;}_{\mathbb{k}}$ and $n^{a}_{\bf k}\!=a^{\dagger}_{\mathbb{k}}a^{\;}_{\mathbb{k}}$. In order to satisfy the $su$(1,1) algebra $[K_{\mathbb{k,Q}}^{-},K_{\mathbb{k^{\prime},Q}}^{+}]=2\delta_{\mathbb{k},\mathbb{k^{\prime}}}K_{\mathbb{k,Q}}^{z}$ and $[K_{\mathbb{k,Q}}^{z},K_{\mathbb{k^{\prime},Q}}^{\pm}]=\pm\delta_{\mathbb{k},\mathbb{k^{\prime}}}K_{\mathbb{k,Q}}^{\pm}$, and to avoid double counting, we restrict momenta $\bf k$ and $\bf 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 $\bf k$ values in the Brillouin zone. The operator $K_{\mathbb{k,Q}}^{+}$, that creates a two-species pair with center-of-mass momentum $\bf 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_{\mathbb{k,Q}}^{+}$ creates a singlet state. The pseudo-spin $1/2$ bosons define an independent and commuting $su$(2) spin algebra generated by $S_{\mathbb{k,Q}}^{z}=(n^{b}_{\bf k+Q}+n^{b}_{\bf-k}-n^{a}_{\bf k+Q}-n^{a}_{\bf-k})/2$, $S_{\mathbb{k,Q}}^{+}=b_{\mathbb{k+Q}}^{\dagger}a^{\;}_{\mathbb{k+Q}}+b_{-\mathbb{k}}^{\dagger}a^{\;}_{-\mathbb{k}}$, $S_{\mathbb{k,Q}}^{-}=(S_{\mathbb{k,Q}}^{+})^{\dagger}$. Although we will focus on the ${\bf Q}={\bf 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 Duk2006 for fermionic systems.

In terms of the $su$(1,1) generators (1), the hyperbolic integrals of motion for ${\bf Q}={\bf 0}$ are Dukelsky2001 ; Ortiz2005

[TABLE]

where $\eta_{\mathbb{k}}$ are arbitrary odd functions of ${\bf k}$. The sum $\mathbb{k}^{\prime}>0$ means that the component $k^{\prime}_{x}$ should be positive.

For a fixed number of bosons $N=2M+\nu$, where $M$ is the number of singlet boson pairs and $\nu$ the total number of unpaired bosons, the eigenvalues of the integrals of motion are

[TABLE]

where $d_{\mathbb{k}}=\nu_{\mathbb{k}}/2+1$, $\nu_{\mathbb{k}}$ is the seniority quantum number (number of unpaired bosons) of level $\mathbb{k}$, and $\nu=\sum_{\mathbb{k}>0}\nu_{\mathbb{k}}$. The spectral parameters $e_{\alpha}$, so-called pairons, are roots of the Richardson equations $(\alpha=1,\dots,M)$

[TABLE]

with

[TABLE]

Each independent solution of the Richardson equations (3) defines a common eigenstate of the integrals of motion (2):

[TABLE]

where the state $|\nu\rangle$, with $\nu$ unpaired bosons, satisfies $\hat{K}_{\mathbb{k}}^{-}|\nu\rangle=0$ for all $\mathbb{k}$, and $\hat{K}_{\mathbb{k}}^{z}|\nu\rangle=d_{\mathbb{k}}|\nu\rangle$.

By combining the integrals of motion $R_{\bf k}$ with the Hellmann-Feynman theorem Rombouts2010 , the occupation probabilities can be obtained from the expectation value

[TABLE]

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 one-dimensional version of the $p$-wave model. It is straightforward to extend our model to higher dimensions as has been done in the fermionic case Sierra2009 ; Rombouts2010 .

III The $p$-wave Bose Hamiltonian

The $p$-wave pairing Bose Hamiltonian we want to study is given by

[TABLE]

where $\eta_{k}=\sin(k/2)$ and $\eta_{k}^{2}=(1-\cos k)/2$. Assuming antiperiodic boundary conditions the allowed $k$ values are $k=\pm\pi/2L,\pm 3\pi/2L,\dots,\pm(2\pi L-\pi)/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

[TABLE]

can be derived from the hyperbolic $su$(1,1) RG integrals of motion (2), by using the linear combination

[TABLE]

where $x=[1+2\lambda(M+L+(\nu/2)-1)]^{-1}$, $Z=\sum_{k>0}\eta_{k}^{2}$, and $G=4\lambda/(1+2\lambda(M-1+\sum_{k>0}d_{k}))$.

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}=\sum_{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}\neq 0$) is contained in our exact solution. For instance, an excess of $a$ atoms manifests through the seniorities $\nu_{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

[TABLE]

Let us analyze next the way the ground state evolves as a function of coupling strength $G\geq 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 $\nu=0$) for weak coupling $G$ has the pairons distributed in the real interval between zero and the minimum $\eta_{\pi/2L}^{2}=\sin^{2}(\pi/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

[TABLE]

where the exact (non-normalized) ground state becomes a condensate of singlet pairs

[TABLE]

which is algebraically analogous to the MR state of the $p+ip$ fermionic model Sierra2009 ; Rombouts2010 , 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\leq G\leq G_{c}$, corresponds to a zero-momentum condensate for each species

[TABLE]

since, as we will see, the quasi-particle gap $\Delta$ 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 Kuklow ; Leggett , which in momentum space becomes

[TABLE]

with $k_{\sf min}=\pi/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,\eta^{2}_{\pi/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_{\sf 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 Leggett .

For $G>G_{c}$ the pairons distribute along an arc that expands in the complex plane as $G$ increases (Fig. 2). At

[TABLE]

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\sum_{k>0}\eta^{2}_{k}$. Infinite pairon energies have been observed previously in fermionic hyperbolic models Ortiz2014 and they were related to a duality associated to the particle-hole symmetry Links . At this point the exact ground state can be expressed as a different pair condensate

[TABLE]

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 $\rho=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 $\eta^{2}_{\pi/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_{\infty}$, 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\rightarrow\infty$ with constant density $\rho=N/(2L)$ and rescaled interaction strength $g=GL$. In this limit, the Richardson equations (3) transform into the boson gap and number equations Rombouts2010 ; Ortiz2005-2

[TABLE]

with quasi-boson energies $E_{k}$ and occupation probabilities $v_{k}^{2}$

[TABLE]

where $\mu$ is the chemical potential and $\Delta$ 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 $\mu\approx-\gamma_{1}g$ and $\Delta\approx\gamma_{2}g$, with $\gamma_{1,2}$ positive constants satisfying $4\gamma_{2}^{2}<\gamma_{1}^{2}$. This latter condition guarantees that the quasi-boson energies, given by $E_{k}\approx g\gamma_{1}\sqrt{1-(4\gamma_{2}^{2}/\gamma_{1}^{2})\eta_{k}^{2}}$ are always real, even in the limit $g\rightarrow\infty$. The ground state energy density ${\cal E}\equiv E/L$ for a given density $\rho$ in the thermodynamic limit is given by

[TABLE]

The critical coupling of the exact solution in the finite-size case, becomes $g_{c}=\lim_{L,N\rightarrow\infty}[G_{c}L]=2/(2+\rho)$ in the thermodynamic limit. The gap $\Delta$ is zero at weak pairing up to the critical value $g_{c}$. The inset of Fig. 1 shows the behavior of the gap for $\rho=0.2$. It increases monotonically for $g>g_{c}$. In the same thermodynamic limit, the coupling where all pairons diverge becomes $g_{\infty}=\lim_{L,N\rightarrow\infty}[G_{\infty}L]=2/\rho$ (Fig. 1).

We are interested in establishing the nature of the non-analyticities of ${\cal E}$ at the critical point. It turns out that ${\cal 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 Rombouts2010 ; Lerma2011 . Close to $g-g_{c}\approx 0^{+}$, it behaves as

[TABLE]

where $\tilde{g}=(g-g_{c})/g_{c}$ (See Appendix A). Interestingly, the behavior of $\cal E$ close to $g_{c}$ depends on $\rho$ 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 $\partial_{g}^{3}{\cal E}\rfloor_{g-g_{c}\rightarrow 0^{+}}=-2\pi^{2}/g_{c}^{6}$, signaling a discontinuity of third order. This is illustrated in Fig. 3 for $\rho=0.2$ where, moreover, ${\cal 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 $\mu+2\Delta^{2}\leq 0$, $E_{k}$ is a monotonous increasing function with minimum at $k_{\sf min}=0$ and energy $E_{k_{\sf min}}=|\mu|$ ($g=1.8$). The previous condition is fulfilled in the superfluid phase only for small densities $\rho<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 Sierra2009 defined by a superfluid with the minimum quasi-boson energy at $k=0$, is given by $\mu+2\Delta^{2}=0$. For $0<\mu+2\Delta^{2}<1$, $E_{k}$ has a minimum at $k_{\sf min}=2\arcsin(\sqrt{\mu+2\Delta^{2}})$, satisfying $0<k_{\sf min}<\pi$ ($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 $\mu+2\Delta^{2}\geq 1$ (area with horizontal lines in Fig. 1), the quasi-boson dispersion is a monotonous decreasing function with minimum at $k_{\sf min}=\pi$ ($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_{\sf min}$ resulting in a delta distribution in the thermodynamic limit. At $g_{c}$, in that limit, the macroscopic occupation at $k_{\sf min}\rightarrow 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 Rombouts2010 ; Ortiz2014 and the $s$-$d$ RG model of Ref. Claeys2018 . In the fermionic case, this fact has been linked to a topological phase transition Rombouts2010 ; Ortiz2014 . 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=\pi$.

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 ${\bf Q}$ has been studied in depth. Imbalanced binary mixtures ($\nu\neq 0$, $m_{a}\neq m_{b}$) and finite ${\bf 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 Sierra2009 ; Rombouts2010 ; VanNeck2014 ; Links 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 Ortiz2014 . 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 Petrov , 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 Kadau and Bose (potassium) mixtures Semeghini ; Cabrera . Can one obtain a quantum liquid phase in $p$-wave Bose systems? This question has been recently addressed in Li2019 , 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.

Acknowledgments— S.L.-H. acknowledges financial support from the Mexican CONACyT project CB2015-01/255702. J.D. is supported by the Spanish Ministerio de Ciencia, Innovación y Universidades, and the European regional development fund (FEDER) under Projects No. FIS2015-63770-P and PGC2018-094180-B-I00, S.L.-H. and J.D. acknowledges financial support .from the Spanish collaboration Grant I-COOP2017 Ref:COOPB20289. G.O acknowledges support from the US Department of Energy grant DE-SC0020343.

Appendix A Non-analytic behavior at the quantum critical point

One can write the boson gap and number equations (13), in the thermodynamic limit, as

[TABLE]

where the following change of variables has been performed: $x=\eta_{k}^{2}$, and $a=\mu+2\Delta^{2}$.

We are interested in characterizing the behavior of physical quantities, such as the chemical potential $\mu$, superfluid gap $\Delta$, and ground state energy density $\cal E$, near the phase transition $g\approx g_{c}$ where a non-analyticity develops. Close to the transition, and for couplings $g>g_{c}$, $\mu<0$ and $a>0$, such that $a>\delta=\mu^{2}-a^{2}>0$. We need to determine the behavior of the above integrals in the limit $\delta\rightarrow 0^{+}$. A few algebraic steps lead to:

[TABLE]

Similar manipulations result in

[TABLE]

Therefore, the resulting gap and number equations close to the critical point become

[TABLE]

or equivalently

[TABLE]

and whose consistency can be checked by taking the limit $a\rightarrow 0$, $g\rightarrow g_{c}$. This gives $\rho+1=2/g_{c}$, as expected from the exact solution. On the other hand, we would like to determine the behavior of the gap $\Delta$ and chemical potential $\mu$ as a function of $\rho$ and $g$ close to the transition. It turns out to be more convenient to write $a-\Delta^{2}=\mu+\Delta^{2}=\mu\,\Gamma$, and find solutions for $\mu$ and $\Gamma$

[TABLE]

What is the behavior of the ground state energy density $\cal E$, Eq. (15),

[TABLE]

close to the phase transition? Following the same strategy, close to the transition point,

[TABLE]

and to first order in powers of $\Gamma$, it results

[TABLE]

where $\tilde{g}=(g-g_{c})/g_{c}$, displaying a discontinuity of third order as indicated in the main text.