Spontaneous formation of polar superfluid droplets in a p-wave interacting Bose gas

TL;DR

This paper investigates how quantum fluctuations stabilize a novel polar superfluid phase in a p-wave interacting Bose gas, leading to spontaneous symmetry breaking and the formation of anisotropic superfluid droplets.

Contribution

It reveals that quantum corrections in p-wave Bose gases can stabilize unstable mean-field phases and predicts the emergence of anisotropic polar superfluid droplets.

Findings

Quantum fluctuations stabilize the phase above a critical density.

The superfluid order parameter has opposite finite momenta for atomic species.

Predicted formation of anisotropic polar superfluid droplets.

Abstract

We study the quantum fluctuations in the condensates of a mixture of bosonic atoms and molecules with interspecies p-wave interaction. Our analysis shows that the quantum phase of coexisting atomic and molecular condensates is unstable at the mean-field level. Unlike the mixture of s-wave interaction, the Lee-Huang-Yang correction of p-wave interaction is unexpectedly found here to exhibit an opposite sign with respect to its mean-field term above a critical particle density. This quantum correction to the mean-field energy provides a remarkable mechanism to self-stabilize the phase. The order parameter of this superfluid phase carries opposite finite momenta for the two atomic species while the molecular component is a polar condensate. Such a correlated order spontaneously breaks a rich set of global U(1) gauge, atomic spin, spatial rotation and translation, and time-reversal symmetries. For potential experimental observation, the phenomenon of anisotropic polar superfluid droplets is predicted to occur, when the particle number is kept finite.

Tables

Table 1. Table 1: Symmetry transformation. U N ​ ( 1 ) subscript 𝑈 𝑁 1 U_{N}(1) : θ ∈ [ 0 , 2 ​ π ) 𝜃 0 2 𝜋 \theta\in[0,2\pi) is an arbitrary angle. This symmetry correponds to the total number conservation. [ S ​ U ​ ( 2 ) / U y ​ ( 1 ) ] delimited-[] 𝑆 𝑈 2 subscript 𝑈 𝑦 1 [SU(2)/U_{y}(1)] with spin rotation symmetry U y ​ ( 1 ) subscript 𝑈 𝑦 1 U_{y}(1) generated by σ y subscript 𝜎 𝑦 \sigma_{y} : θ x subscript 𝜃 𝑥 \theta_{x} and θ z subscript 𝜃 𝑧 \theta_{z} are arbitrary angles. Here σ x , y , z subscript 𝜎 𝑥 𝑦 𝑧 \sigma_{x,y,z} are the Pauli matrices. S ​ O ​ ( 3 ) 𝑆 𝑂 3 SO(3) : λ x , y , z subscript 𝜆 𝑥 𝑦 𝑧 \lambda_{x,y,z} are defined in Eq. ( 2 ) and θ ¯ x , y , z subscript ¯ 𝜃 𝑥 𝑦 𝑧 \overline{\theta}_{x,y,z} are arbitrary rotation angles. T ​ r 𝑇 𝑟 Tr : r ′ superscript r ′ \textbf{r}^{\prime} is an arbitrary displacement vector in 3D spatial coordinate. 𝒯 𝒯 \mathcal{T} time-reversal: We use momentum representation to expand ψ ^ 1 subscript ^ 𝜓 1 \hat{\psi}_{1} and ψ ^ 2 subscript ^ 𝜓 2 \hat{\psi}_{2} fields. Due to momentum conservation, the momentum of molecule fields is restricted to p 1 + p 2 subscript p 1 subscript p 2 \textbf{p}_{1}+\textbf{p}_{2} .

symmetry	${\hat{\psi }}_1\left(\text{r}\right)$	${\hat{\psi }}_2\left(\text{r}\right)$	${\hat{\varphi }}_{x,y,z}\left(\text{r}\right)$	$\nabla$
$U_N\left(1\right)$	$e^{i\theta }{\hat{\psi }}_1$	$e^{i\theta }{\hat{\psi }}_2$	$e^{2i\theta }{\hat{\varphi }}_{x,y,z}$	$-$
$\left[SU\left(2\right)/U_y\left(1\right)\right]$	$e^{i{\theta }_x{\sigma }_x+i{\theta }_z{\sigma }_z}{({\hat{\psi }}_1,{\hat{\psi }}_2)}^T$		$-$	$-$
$SO\left(3\right)$	$-$	$-$	$e^{i{\sum }_{i=x,y,z}{\overline{\theta }}_i{\lambda }_i}\hat{\varphi }$	$e^{i{\sum }_{i=x,y,z}{\overline{\theta }}_i{\lambda }_i}\nabla$
$Tr$	${\hat{\psi }}_1\left(\text{r}+{\text{r}}^{\prime }\right)$	${\hat{\psi }}_2\left(\text{r}+{\text{r}}^{\prime }\right)$	${\hat{\varphi }}_{x,y,z}\left(\text{r}+{\text{r}}^{\prime }\right)$	$-$
$𝒯$	${\sum }_{p_1}{e}^{-i{\text{p}}_1\cdot \text{r}}{\hat{a}}_{1,-{\text{p}}_1}$	${\sum }_{p_2}{e}^{-i{\text{p}}_2\cdot \text{r}}{\hat{a}}_{2,-{\text{p}}_2}$	$e^{-i\left({\text{p}}_1+{\text{p}}_2\right)\cdot \text{r}}{\hat{b}}_{i,-{\text{p}}_1-{\text{p}}_2}$	$-$

Table 2. Table 2: Table of ground state phases. Here we have three phases by setting different detuning. ASF, AMSF and MSF are the atomic, atomic-molecular and molecular condensate phases, respectively.

Phase	${ϵ}_0$	$n_M$	$n_1=n_2$	$Q$	$\mu$	$Z$	$E_0/V$
ASF		$0$	$\frac{1}{2}n$	$0$	$0$	$0$	$0$
AMSF		$\frac{1}{4}n+\frac{{ϵ}_0}{2{\overline{g}}^{2}m}$	$\frac{1}{4}n-\frac{{ϵ}_0}{2{\overline{g}}^{2}m}$	$-\overline{g}m\sqrt{\frac{1}{4}n+\frac{{ϵ}_0}{2{\overline{g}}^{2}m}}$	$-\frac{1}{8}{\overline{g}}^{2}mn-\frac{1}{4}{ϵ}_0$	$0$	$-\frac{1}{16{\overline{g}}^{2}m}{\left({\overline{g}}^{2}mn+2{ϵ}_0\right)}^2$
MSF	${ϵ}_0>\frac{1}{2}{\overline{g}}^{2}mn$	$\frac{1}{2}n$	$0$	$-\frac{1}{\sqrt{2}}\overline{g}m\sqrt{n}$	$-\frac{1}{4}{\overline{g}}^{2}mn$	$0$	$-\frac{1}{2}{ϵ}_{0}n$

Equations

$$\begin{array}{rl}\mathcal{H}=& \sum _{\sigma =1,2}{\hat{\psi }}_{\sigma }^{†}(-\frac{{\nabla }^2}{2m}){\hat{\psi }}_{\sigma }+\sum _{i=x,y,z}{\hat{\varphi }}_i^{†}(-\frac{{\nabla }^2}{4m}-{ϵ}_0){\hat{\varphi }}_i\\ & +\sum _{i=x,y,z}[\frac{\overline{g}}{2}{\hat{\varphi }}_i^{†}({\hat{\psi }}_1,{\hat{\psi }}_2){\sigma }_y{\partial }_i({\hat{\psi }}_1,{\hat{\psi }}_2{)}^T+h.c.],\end{array}$$

$$\begin{array}{rl}& {\lambda }_x=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& -i\\ 0& i& 0\end{array}\right),{\lambda }_y=\left(\begin{array}{ccc}0& 0& i\\ 0& 0& 0\\ -i& 0& 0\end{array}\right),{\lambda }_z=\left(\begin{array}{ccc}0& -i& 0\\ i& 0& 0\\ 0& 0& 0\end{array}\right).\end{array}$$

$$N_1+N_2+2N_M=N,N_1-N_2=\delta N,$$

$${\Psi }_1={\Psi }_{1,Q}{e}^{-i\text{Q}\cdot \text{r}},{\Psi }_2={\Psi }_{2,-Q}{e}^{i\text{Q}\cdot \text{r}},$$

$$\begin{array}{rl}& \Psi =\sqrt{n_A}{e}^{i\theta }{e}^{i({\theta }_x{\sigma }_x+{\theta }_z{\sigma }_z)}\left(\begin{array}{c}\cos {\chi }_A{e}^{-i\text{Q}\cdot \text{r}}\\ \sin {\chi }_A{e}^{i\text{Q}\cdot \text{r}}\end{array}\right),\\ & \text{Q}=e^{i{\sum }_{i=x,y,z}{\overline{\theta }}_i{\lambda }_i}{\text{Q}}_0,\\ & \Phi =\sqrt{n_M}{e}^{i(2\theta +{\theta }_M)}{e}^{i{\sum }_{i=x,y,z}{\overline{\theta }}_i{\lambda }_i}\left(\begin{array}{c}\cos {\chi }_M\\ i\sin {\chi }_M\\ 0\end{array}\right),\end{array}$$

$$\begin{array}{rl}F/V=& \sum _{\sigma =1,2}\frac{Q^2}{2m}{n}_{\sigma }-{ϵ}_0{n}_M-\mu (n_1+n_2+2n_M-n)\\ & +\frac{\overline{g}}{2}{n}_A\sqrt{n_M}\sin 2{\chi }_A[e^{-i{\theta }_M}(\cos {\chi }_M,-i\sin {\chi }_M,0)\\ & \cdot {\text{Q}}_0+h.c]-Z(n_1-n_2),\end{array}$$

$${\hat{\psi }}_{\sigma }={\Psi }_{\sigma }+\delta {\hat{\psi }}_{\sigma },{\hat{\varphi }}_i={\Phi }_i+\delta {\hat{\varphi }}_i,$$

$$\delta {\hat{\psi }}_{\sigma }=\frac{1}{\sqrt{V}}\sum _{\text{k}}\delta {\hat{a}}_{\sigma ,\text{k}}{e}^{-i\text{k}\cdot \text{r}},\delta {\hat{\varphi }}_i=\frac{1}{\sqrt{V}}\sum _{\text{k}}\delta {\hat{b}}_{i,\text{k}}{e}^{-i\text{k}\cdot \text{r}},$$

$$\begin{array}{rl}{H}_f=& \frac{1}{2}\sum _{\text{k}}\{\sum _{\sigma =1,2}{\epsilon }_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}\delta {\hat{a}}_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}^{†}\delta {\hat{a}}_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}\\ & +\sum _i{\omega }_{i,k}\delta {\hat{b}}_{i,k}^{†}\delta {\hat{b}}_{i,k}-2\sum _{\sigma ,i}{\alpha }_{i,\overline{\sigma },\text{k}}\delta {\hat{b}}_{i,k}^{†}\delta {\hat{a}}_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}\\ & +h.c.\},\end{array}$$

$$\begin{array}{rl}& {\epsilon }_{\sigma ,\text{k}}={ϵ}_k=\frac{k^2}{2m},{\omega }_{i,\text{k}}=\frac{1}{2}{ϵ}_k+\frac{1}{2}{\overline{g}}^{2}mn,\\ & {\alpha }_{i,\sigma ,\text{k}}=\mp \frac{1}{2\sqrt{2}}\overline{g}\sqrt{n}{k}_i.\end{array}$$

$$E_i^{ASF}=\frac{1}{4}{\overline{g}}^{2}mn{f}_i^{ASF}(\tilde{x},\tilde{y},\tilde{z}),i=1,\cdots ,5$$

$$\begin{array}{rl}{H}_f=& \sum _{\text{k}}\{\sum _{\sigma =1,2}\frac{1}{2}{\epsilon }_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}\delta {\hat{a}}_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}^{†}\delta {\hat{a}}_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}\\ & +\sum _i\frac{1}{2}{\omega }_{i,k}\delta {\hat{b}}_{i,k}^{†}\delta {\hat{b}}_{i,k}+t_{\text{k}+\text{Q}}\delta {\hat{a}}_{1,\text{k}+\text{Q}}\delta {\hat{a}}_{2,-\text{k}-\text{Q}}\\ & -\sum _{\sigma ,i}{\alpha }_{i,\overline{\sigma },\text{k}}\delta {\hat{b}}_{i,k}^{†}\delta {\hat{a}}_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}+h.c.\},\end{array}$$

$$\begin{array}{rl}& {\epsilon }_{\sigma ,\text{k}}={ϵ}_k+\frac{1}{8}{\overline{g}}^{2}mn+\frac{1}{4}{ϵ}_0,{\omega }_{i,\text{k}}=\frac{1}{2}{ϵ}_k+\frac{1}{4}{\overline{g}}^{2}mn-\frac{1}{2}{ϵ}_0,\\ & t_{\text{k}}=-\overline{g}\sum _i{\Phi }_i^{\ast }{k}_i,{\alpha }_{i,\overline{\sigma },\text{k}}=\pm \overline{g}\sqrt{n_{\sigma }}(Q_{\sigma ,i}-k_i/2),\end{array}$$

$$E_i^{AMSF}=\frac{1}{4}{\overline{g}}^{2}mn{f}_i^{AMSF}(\tilde{x},\tilde{y},\tilde{z},\Delta ),i=1,\cdots ,5$$

$$\begin{array}{rl}{H}_f=& \sum _{\text{k}}\{\sum _{\sigma =1,2}\frac{1}{2}{\epsilon }_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}\delta {\hat{a}}_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}^{†}\delta {\hat{a}}_{\sigma ,\text{k}+{\text{Q}}_{\sigma }}\\ & +t_{\text{k}+\text{Q}}\delta {\hat{a}}_{1,\text{k}+\text{Q}}\delta {\hat{a}}_{2,-\text{k}-\text{Q}}\\ & +\sum _i\frac{1}{2}{\omega }_{i,k}\delta {\hat{a}}_{i,k}^{†}\delta {\hat{a}}_{i,k}+h.c.\},\end{array}$$

$$\begin{array}{rl}& E_1^{MSF}=\frac{1}{4}{\overline{g}}^{2}mn\sqrt{{\tilde{r}}^2(2+{\tilde{r}}^2-2\sqrt{2}\tilde{r}\cos \gamma )},\\ & E_2^{MSF}=\frac{1}{4}{\overline{g}}^{2}mn\sqrt{{\tilde{r}}^2(2+{\tilde{r}}^2+2\sqrt{2}\tilde{r}\cos \gamma )},\\ & E_{3,4,5}^{MSF}=\frac{1}{8}{\overline{g}}^{2}mn{\tilde{r}}^2.\end{array}$$

$$E_i^{MSF}=\frac{1}{4}{\overline{g}}^{2}mn{f}_i^{MSF}(\tilde{x},\tilde{y},\tilde{z},\Delta ),i=1,\cdots ,5$$

$$\begin{array}{rl}& -i{T}_{\text{k},{\text{k}}^{\prime }}^{(l_z)}(k)=D^{(0)}(k)(-i\overline{g}{)}^2{k}^2{Y}_{1,l_z}(\hat{\text{k}})Y_{1,l_z}^{\ast }({\hat{\text{k}}}^{\prime })\\ & +D^{(0)2}(k)(-i\overline{g}{)}^{4}2{\Pi }_{l_z}(k)k^2{Y}_{1,l_z}(\hat{\text{k}})Y_{1,l_z}^{\ast }({\hat{\text{k}}}^{\prime })+...\\ & =D(k)(-i\overline{g}{)}^2{k}^2{Y}_{1,l_z}(\hat{\text{k}})Y_{1,l_z}^{\ast }({\hat{\text{k}}}^{\prime }).\end{array}$$

$$D^{(0)}(k)=\frac{i}{k^2/m+{ϵ}_0+i0+},$$

$${\Pi }_{l_z}(k)=\int \frac{d^{3}p}{(2\pi {)}^3}\frac{i{p}^2|Y_{1,l_z}(\hat{\text{p}}){|}^2}{k^2/m-p^2/m+i0+}.$$

$$D^{-1}(k)=[D^{(0)}(k){]}^{-1}-(-i\overline{g}{)}^2{\Pi }_{l_z}(k).$$

$$\frac{{\widetilde{ϵ}}_0}{{\tilde{\overline{g}}}^2}=\frac{{ϵ}_0}{{\overline{g}}^2}+\int \frac{d^{3}p}{(2\pi {)}^3}m|Y_{1,l_z}(\hat{\text{p}}){|}^2,$$

$$\frac{1}{{\tilde{\overline{g}}}^2}=\frac{1}{{\overline{g}}^2}-\int \frac{d^{3}p}{(2\pi {)}^3}{m}^2\frac{|Y_{1,l_z}(\hat{\text{p}}){|}^2}{p^2},$$

$$E_{0,r}^{ASF}/V=0,$$

$$E_{0,r}^{AMSF}/V=\int \frac{d^{3}k}{(2\pi {)}^3}\frac{1}{12}{\overline{g}}^{2}mn+\frac{1}{6}{ϵ}_0+\frac{{\overline{g}}^4{m}^2{n}^2-4{ϵ}_0^2}{72k^2},$$

$$E_{0,r}^{MSF}/V=\int \frac{d^{3}k}{(2\pi {)}^3}\frac{1}{6}{\overline{g}}^{2}mn+\frac{{\overline{g}}^{2}mn-2{ϵ}_0}{12k^2}.$$

$$\begin{array}{rl}{E}_{LHY}^{ASF}/V=& -{\overline{g}}^5{m}^4{n}^{2.5}\int \frac{d^3\tilde{r}}{(2\pi {)}^3}\sum _{i=1}^5\frac{1}{4}{f}_i^{ASF}\\ & -\frac{7}{8}{\tilde{r}}^2-\frac{3}{8},\end{array}$$

$$\begin{array}{rl}{E}_{LHY}^{AMSF}/V=& -{\overline{g}}^5{m}^4{n}^{2.5}\int \frac{d^3\tilde{r}}{(2\pi {)}^3}\sum _{i=1}^5\frac{1}{4}{f}_i^{AMSF}\\ & -\frac{7}{8}{\tilde{r}}^2-\frac{13}{24}+\frac{5}{6}\Delta +\frac{1-16{\Delta }^2}{72{\tilde{r}}^2},\end{array}$$

$$\begin{array}{rl}{E}_{LHY}^{MSF}/V=& -{\overline{g}}^5{m}^4{n}^{2.5}\int \frac{d^3\tilde{r}}{(2\pi {)}^3}\sum _{i=1}^5\frac{1}{4}{f}_i^{MSF}\\ & -\frac{7}{8}{\tilde{r}}^2-\frac{1}{3}+\frac{1-4\Delta }{12{\tilde{r}}^2}.\end{array}$$

$$E_g^{ASF}/V=0,$$

Full text

I Introduction

Quantum fluctuation is one of the most intrinsic properties of quantum mechanics, which is responsible for many fascinating physical phenomena, such as Casimir effect and abundant quantum phase transitions. Recently, Petrov showed that quantum fluctuation reflected by Lee-Huang-Yang (LHY) correction can prevent a mean-field-unstable Bose gas from collapsing Petrov (2015). The competition between the mean-field attraction and LHY repulsion stabilizes the Bose gas into a self-bound liquidlike droplet state. Subsequently, several experimental groups reported this novel quantum state with the prediction of Petrov Semeghini et al. (2018); Cheiney et al. (2018); Cabrera et al. (2018). In order to protrude the action of LHY correction, which is typically small in the dilute limit, Petrov suggested to subtly balance the inter- and intra-species interactions at the mean-field level. Owing to its unique formation mechanism, the self-bound state shows many interesting features, such as the quantum droplet is self-trapped and evaporated without external potential Petrov (2015).

The properties of quantum droplet are linked to the properties of interaction between particles. It is natural to ask if quantum droplet can be stabilized with other types of interaction and what their properties might be. It was also found that quantum droplets can be stabilized in a dipolar Bose gas benefiting from the competition between the dipolar interaction and s-wave contact interaction Ferrier-Barbut et al. (2016); Kadau et al. (2016); Chomaz et al. (2016); Schmitt et al. (2016). The quantum droplets in a dipolar Bose gas are anisotropic and form a regular array, as a consequence of the dipolar interaction is anisotropic and long-ranged. Morover, it is also predicted quantum droplets can be stabilized with the assistance of three-body interaction Bulgac (2002); Sekino and Nishida (2018) and spin-orbit coupling Cui (2018).

Here we study the beyond-mean-field ground state of a p-wave interacting Bose gas, and predict the existence of finite-momentum anisotropic self-stabilized quantum droplet. At the mean-field level, this p-wave interacting Bose gas typically has three ground-state phases: atomic superfluid (ASF) phase with only the atomic condensate, atomic-molecular superfluid (AMSF) phase with both atomic and molecular condensates, and molecular superfluid (MSF) phase with only the molecular condensate. We find AMSF phase is unstable and tends to collapse. Unlike pure s-wave interaction T. D. Lee and Yang (1957), we find the sign of the LHY correction of p-wave interaction may be different from that of the mean-field term when varying particle densities. A balance between the mean-field part and LHY correction exists for certain particle density, which gives rise to a self-stabilized (-bound) state without external potential. It is shown the self-stabilized state even survives in the dilute limit estimated with scattering volume. In addition to the U(1) global phase symmetry, the rotation, translation and time-reversal symmetries are found to be spontaneously broken by the presence of finite momentum of the order parameters. The result ground state is predicted to be an anisotropic quantum droplet with finite momentum for a system with finite particle number.

II Model

Inspired by the experimental observations of p-wave Feshbach resonance in the mixture of 85Rb and 87Rb atoms Papp et al. (2008); Dong et al. (2016), we consider a mixture of two distinguishable species of bosonic atoms respectively created by $\hat{\psi}_{1}^{\dagger}(\boldsymbol{r})$ and $\hat{\psi}_{2}^{\dagger}(\boldsymbol{r})$ with interspecies p-wave interaction. The p-wave interaction arises from a p-wave Feshbach resonance by coupling with three closed molecular channels denoted by $l_{z}=-1,0,1$. Here $l_{z}\hbar$ are the magnetic angular momentum carried by the molecules on the closed channels, which are created by $\hat{\phi}_{l_{z}=-1,0,1}^{\dagger}(\boldsymbol{r})$ respectively. It will be convenient to discuss the physics with bases $\hat{\phi}_{i=x,y,z}^{\dagger}$, which are related with $\hat{\phi}_{l_{z}=-1,0,1}^{\dagger}$ through $\phi_{\pm 1}^{\dagger}=(\phi_{x}^{\dagger}\pm i\phi_{y}^{\dagger})/\sqrt{2}$, and $\phi_{0}^{\dagger}=\phi_{z}^{\dagger}$. To focus on the physics arising from p-wave interaction, we will restrict our attention on the case where the closed channels are degenerate and background (non-resonant) interactions are neglectable. The system we consider is characterized by Hamiltonian density

[TABLE]

where the atomic masses have been assumed to be the same, i.e. $m_{1}=m_{2}=m$, $\epsilon_{0}$ is the detuning of molecule channels, $\bar{g}$ represents the strength of p-wave interaction, and $\sigma_{y}$ is the Pauli matrix. Here the reduced Plank constant $\hbar$ has been set as $1$.

Our model possesses $U_{N}(1)\times[SU(2)/U_{y}(1)]\times SO(3)\times Tr\times\mathcal{T}$ symmetries, where $U_{N}(1)$ is the global gauge symmetry, $[SU(2)/U_{y}(1)]$ the spin rotation symmetry around $x$ and $z$ directions, $SO(3)$ the 3-dimensional spatial rotation symmetry, $Tr$ the translation symmetry in the absence of an external field, and $\mathcal{T}$ the time reversal symmetry. The symmetry transformations are listed in Tab. 1. It is worth noting that spin-rotation symmetry $[SU(2)/U_{y}(1)]$ is reduced to a spin-rotation symmetry $U_{z}(1)$ generated by $\sigma_{z}$ in presence of intraspecies s-wave interaction Radzihovsky and Choi (2009); Choi and Radzihovsky (2011). In $SO(3)$ rotation symmetry, the atom fields are scalar fields, so they remain constant under $SO(3)$ transformation. However, molecular field $\hat{\phi}$ and gradient operator $\boldsymbol{\nabla}$ are all vector fields, and they are transformed by a 3D spatial rotation. In Tab. 1, the generators of rotation symmetry $\lambda_{x,y,z}$ are given by,

[TABLE]

Time-reversal symmetry $\mathcal{T}$ is given by reversing the momentum of atomic and molecular field operators, i.e. transforming $\hat{a}_{1,\textbf{p}_{1}},\hat{a}_{2,\textbf{p}_{2}}$, and $\hat{b}_{i,\textbf{p}_{1}+\textbf{p}_{2}}$ as $\hat{a}_{1,-\textbf{p}_{1}},\hat{a}_{2,-\textbf{p}_{2}}$, and $\hat{b}_{i,-\textbf{p}_{1}-\textbf{p}_{2}}$, respectively.

The total particle number $N$ and atomic number difference $\delta N$ are defined as below,

[TABLE]

where we use $N_{1,2}=\int d^{3}r\langle\hat{\psi}_{1,2}^{\dagger}\hat{\psi}_{1,2}\rangle$ and $N_{M}=\sum_{i=x,y,z}\int d^{3}r\langle\hat{\phi}^{\dagger}_{i}\hat{\phi}_{i}\rangle$ to denote the numbers of atoms and molecules, respectively. Here $\langle\cdots\rangle$ represents the average over the ground state. Obviously $N$ and $\delta N$ are conserved in our model, which correspond to the $U_{N}(1)$ and $[SU(2)/U_{y}(1)]$ symmetries.

III Mean-Field Ground State

As the foundation of beyond-mean-field study, we need to characterize the ground state at the mean-field level at first. We use the mean fields $\Psi_{1}=\langle\hat{\psi}_{1}\rangle$, $\Psi_{2}=\langle\hat{\psi}_{2}\rangle$ and $\Phi_{i}=\langle\hat{\phi}_{i}\rangle$ to describe the atomic and molecular condensates. The mean-field ground state of a p-wave resonant Bose gas including considerable large intraspecies s-wave interaction has been systematically discussed before Radzihovsky and Choi (2009); Choi and Radzihovsky (2011). Three mean-field stable phases for the ground states: atomic (ASF), atomic-molecular (AMSF) and molecular (MSF) superfluid, are found. Typically, the atomic condensates carry finite momentum due to the p-wave interaction in AMSF phase. Actually, the ground-sate phase diagram of our model is similar to the case there. While due to the lack of intraspecies s-wave interaction (or due to weak intraspecies s-wave interaction), it is shown the ground states of our model may be unstable in the mean field level.

As the typical feature of p-wave interaction, the atomic condensates generally carry finite momentum due to the shift of energy minimum in momentum space by the interaction terms Radzihovsky and Choi (2009); Choi and Radzihovsky (2011). Although a general description of atomic order parameters should be written as $\Psi_{1}=\sum_{\textbf{Q}_{n}}\Psi_{1,\textbf{Q}_{n}}e^{-i\textbf{Q}_{n}\cdot\textbf{r}}$ and $\Psi_{2}=\sum_{\textbf{Q}_{n}}\Psi_{2,-\textbf{Q}_{n}}e^{i\textbf{Q}_{n}\cdot\textbf{r}}$, it is shown that the assumption $\textbf{Q}_{n}=\textbf{Q}$ is sufficient to capture the ground state in presence of intraspecies s-wave interaction Radzihovsky and Choi (2009); Choi and Radzihovsky (2011), i.e.

[TABLE]

due to the lack of $\sigma_{x}$ spin-rotation symmetry for atomic components. Correspondingly the molecular components are space-independent, since the molecular fields only feel a homogeneous potential by atoms. Considering the symmetries of our model, we have the following ground-state ansatz

[TABLE]

where $\theta,\theta_{M}\in[0,2\pi)$ are $U(1)$ phases, $\theta_{x,z}\in[0,2\pi)$ are $[SU(2)/U_{y}(1)]$ spin rotation angles, $\overline{\theta}_{x,y,z}$ are $SO(3)$ rotation angles, $\chi_{A},\chi_{M}\in[0,2\pi)$, $\textbf{Q}_{0}=(Q_{0,x},Q_{0,y},Q_{0,z})^{T}$ is an arbitrary real three dimentional vector, $n_{A}=(N_{1}+N_{2})/V,n_{M}=N_{M}/V$ with system volume $V$ are the total atomic density and molecular density respectively.

Furthermore, we derive the free energy density by substituting the above ansatz (5) to the Hamiltonian density (1)

[TABLE]

where $n_{1,2}=N_{1,2}/V$, $\mu$ and $Z$ are the Lagrange multipliers set for the conservations of the total particle number and atom-number difference. For simplicity, we only consider a nonpolarized situation in this paper, i.e. $n_{1}=n_{2}=n_{A}/2$, and fix the total particle number. The free energy density does not depend on $\theta,\theta_{x,z},\overline{\theta}_{x,y,z}$. To minimize the free energy, we obtain the optimal values for the parameters: $\theta_{M}=0,\chi_{A}=\pi/4,\chi_{M}=0,Q_{0,x}=|\textbf{Q}|,Q_{0,y,z}=0$, from which we can see that $\Phi$ is real and parallel to Q by setting $\theta=0$. To be more convenient, we set $\overline{\theta}_{y}=\pi/2,\overline{\theta}_{x,z}=0$ so that Q and $\Phi$ are aligned to $z$ direction. Without loss of generosity, we choose $\overline{g}$ to be negative (if $\overline{g}>0$, Q will be opposite to $\Phi$, however, it gives us the same phases and LHY corrections as we obtain below). Gross-Pitaevskii (GP) equations can be derived from the free energy density formula, and we obtain the optimized solutions to minimize the free energy.

Similar to previous literatures Radzihovsky and Choi (2009); Choi and Radzihovsky (2011), the ground state phase diagram of our model is also divided into three phases for different detuning $\epsilon_{0}$, where the ground state phases are listed in the Tab. 2. Here ASF refers to the atomic superfluid phase, where only atomic condensates exist. Note that there is no superfluidity here due to the absence of background atom-atom interaction, where the name of phase is only taken to be consistent with previous convention Radzihovsky and Choi (2009); Choi and Radzihovsky (2011). AMSF refers to the atomic-molecular superfluid phase, where atom and molecular condensates are present in the same phase. MSF with only molecular condensate is the molecular superfluid phase.

In ASF phase, the condensate in both atomic species stays stationary due to vanishing Q and the two condensates do not interact. The atomic chemical potential remains zero. In AMSF phase, the rotation and time-reversal symmetries are all broken due to the finite-momentum condensates. The $SO(3)$ rotation symmetry is spontaneously broken into $SO(2)$ symmetry. In MSF phase, although the atomic condensates density is zero, we still have non-zero Q. This results in a MSF excitation spectrum translated in momentum space by Q as we will see in section IV.

From Tab. 2, we can also find the total energy $E_{0}$ is proportional to particle number $N=nV$ in phases ASF and MSF, which is due to the lack of background atom-atom and molecule-molecule interactions in these phases, respectively. It means the total energy $E_{0}$ is constant, such that the ground state is stable, for a system with fixed total particle number. However, we can find it is energetically favorable to increase density $n$ to reach lower total energy $E_{0}$ in AMSF phase. It implies that in this phase the mean-field ground state is unstable and tends to collapse into a state with smaller volume but large particle density when the total particle number is fixed. The instability of the ground state in AMSF phase also manifests itself in the fact that the excitation mode becomes complex in the long-wavelength limit $k\to 0$ Aybar and Oktel (2019). It will be shown the ground state collapses into a small droplet after considering LHY correction T. D. Lee and Yang (1957). In order to calculate this correction, we need to analyze the Bogoliubov excitation spectrum at first.

IV Bogoliubov excitation spectrum

We will study the Bogoliubov excitation spectrum in this section. Following Bogoliubov’s theory Bogolyubov (1947); Pethick and Smith (2008), we expand the atomic and molecular fields around the ground-state mean fields,

[TABLE]

with the fluctuation fields $\delta\hat{\psi}_{\sigma}$ and $\delta\hat{\phi}_{i}$. For convenience, we furthermore expand $\delta\hat{\psi}_{\sigma}$ and $\delta\hat{\phi}_{i}$ with the Fourier transformation

[TABLE]

where $\delta\hat{a}_{\sigma,\textbf{k}}$ and $\delta\hat{b}_{i,\textbf{k}}$ are the corresponding quantum fluctuation fields in momentum space. Substituting Eqs. (7) and (8) into Eq. (1) and keeping only the second-order terms (the first-order terms are vanished due to the saddle-point solution and higher-order terms will be neglected), we can derive the Bogoliubov Hamiltonian. The Bogoliubov excitation spectrum can be extracted by diagonalizing Bogoliubov Hamiltonian.

IV.1 ASF phase

This phase has only atomic condensates, i.e. $n_{1}=n_{2}=n/2$, $n_{M}=0$ and the zero atomic condensates momentum $\textbf{Q}=0$. The Bogoliubov Hamiltonian can be written as

[TABLE]

where $\overline{\sigma}=3-\sigma,\sigma=1,2$, and the parameters are given as

[TABLE]

The corresponding Bogoliubov excitation spectrum is given by

[TABLE]

where $f_{i}^{ASF}$ is a dimensionless function, and $\tilde{x}=\frac{k_{x}}{\overline{g}m\sqrt{n}},\tilde{y}=\frac{k_{y}}{\overline{g}m\sqrt{n}},\tilde{z}=\frac{k_{z}}{\overline{g}m\sqrt{n}}$.

We show $f^{ASF}_{i}$ along the radial direction in Fig. 1(a) and (d). The spectrum is symmetric in all directions and has two gapless atomic modes. The quadratic dispersions of gapless mode are due to the absence of atom-atom interaction.

IV.2 AMSF phase

In AMSF phase, particles are condensed into both the atomic and molecular channels, and the atomic condensates carry opposite finite momentums. The directions of atomic momentum Q and molecular condensates order parameter $\Phi=(\Phi_{x},\Phi_{y},\Phi_{z})$ are parallel in mean-field ground state, where the direction of $\Phi$ is defined by the it three spatial components. For convenience, we build the coordinate so that this direction is aligned along $z$ axis. The Bogoliubov Hamiltonian is written as

[TABLE]

where the parameters are given by

[TABLE]

with $\epsilon_{k}=\frac{k^{2}}{2m}$, $\sigma=1,2$ (correspondingly $\bar{\sigma}=2,1$), $\textbf{Q}_{1}=\textbf{Q}$ and $\textbf{Q}_{2}=-\textbf{Q}$. The Bogoliubov excitation spectrum can be written as

[TABLE]

where $f_{i}^{AMSF}$ is a dimensionless function, $\tilde{x}=\frac{k_{x}}{\overline{g}m\sqrt{n}},\tilde{y}=\frac{k_{y}}{\overline{g}m\sqrt{n}},\tilde{z}=\frac{k_{z}}{\overline{g}m\sqrt{n}}$ and $\Delta=\frac{\epsilon_{0}}{2\overline{g}^{2}mn}$.

The schematic plots of $f_{i}^{AMSF}$ are shown in Fig. 1(b) and (e) along $z$ and $x$ directions respectively. As we can see from the two figures, the blue-dashed curve shows imaginary mode consistent with the instability of the mean-field ground state Aybar and Oktel (2019), which is absent when the ground state is stable Radzihovsky and Choi (2009); Choi and Radzihovsky (2011). Actually, the true ground state is lost due to the homogeneous assumption (the system with finite particle number will collapse into a droplet shape that breaks the spatial translation symmetry) and the absence of LHY correction. On the other hand, the inverse of the largest momentum carried by imaginary modes is expected to be comparable with the size of droplet Aybar and Oktel (2019). The minima on the blue-solid curve in Fig. 1(b) corresponds to the nonvanishing momentum $2\textbf{Q}$ in AMSF phase, where the atomic condensates locate. That the spectrum softens to zero at $k_{z}=2Q$ implies our ansatz correctly captures the feature of the ground state.

IV.3 MSF phase

In this phase, we have $n_{M}=n/2$ and $n_{1}=n_{2}=0$. The Bogoliubov Hamiltonian is given as

[TABLE]

where $\varepsilon_{\sigma,\textbf{k}}=\epsilon_{k}+\frac{1}{4}\overline{g}^{2}mn,\quad\omega_{i,\textbf{k}}=\frac{1}{2}\epsilon_{k}$ and $\quad t_{\textbf{k}}=-\overline{g}\sum_{i}\Phi_{i}^{\ast}k_{i}$.

Fortunately, we can derive analytical formulas for the excitation modes in this phase, i.e.

[TABLE]

where $\tilde{r}^{2}=(k_{x}^{2}+k_{y}^{2}+k_{z}^{2})/\overline{g}^{2}m^{2}n$ and $\gamma$ is the angle between $z$ axis and unit vector $\hat{\textbf{k}}$ as we have aligned $\hat{\textbf{Q}}$ along $z$. Similar to what we defined in ASF and AMSF phases, we rewrite the excitation modes in this formula

[TABLE]

Fig. 1(c) and (f) are the corresponding $f_{i}^{MSF}$ along $z$ and $x$ directions. In Fig. 1(b), the red dotted and green dashed curves are the two atom modes respectively, and the minima on green dashed curve corresponds to the nonvanishing momentum $2\textbf{Q}$. The blue curve denotes the triply-degenerated molecule modes. In Fig. 1(f), the red curve denotes the doubly-degenerated atom modes and blue curve denotes the triply-degenerated molecule modes.

V LHY correction

The LHY correction is the leading-order correction of quantum fluctuation. It is composed of Bogoliubov excitation energies, commutation energies which appear due to the commutation relations of Nambu basis, and energy correction due to the interaction renormalization. Here the interaction renormalization is employed to remove the energy divergence arising in collecting the energy of quantum fluctuation T. D. Lee and Yang (1957). Let us review the renormalization procedure before going ahead.

To remove the divergence appears in the calculation of LHY correction, we need to renormalize the interaction parameter $\overline{g}$ and detuning $\epsilon_{0}$ T. D. Lee and Yang (1957); Aybar and Oktel (2019). The two body T matrix for p wave interaction is given by Gurarie and Radzihovsky (2007)

[TABLE]

where the index $l_{z}$ denotes different interacting channels $l_{z}=-1,0,1$. $Y_{1,l_{z}}(\hat{\textbf{k}})$ is the $l_{z}$-th channel of the first order spherical harmonics. $D^{(0)}(k)$ is the p-wave scattering propagator, and $\Pi_{l_{z}}(k)$ is the polarization bubble for channel $l_{z}$, which are given by

[TABLE]

and

[TABLE]

Using Eq. (18), we yield

[TABLE]

Comparing $k^{0}$ term and $k^{2}$ term on both sides of Eq. (21), we obtain the renormalization relations Qin (2018),

[TABLE]

and

[TABLE]

where $\tilde{\epsilon}_{0}$ and $\tilde{\overline{g}}$ are the renormalized detuning and p-wave interacting strength respectively.

Applying these renormalization relations into the ground state energy in different phases, one obtains the renormalized mean-field ground state energies

[TABLE]

[TABLE]

and

[TABLE]

From the analysis of Bogoliubov spectrum and interaction renormalization, we obtain the LHY correction densities in different phases,

[TABLE]

[TABLE]

and

[TABLE]

Combining the mean-field ground state energy densities and LHY corrected energy densities yields the total ground state energy density $E_{g}/V=E_{0}/V+E_{LHY}/V$ for different phases as follows,

[TABLE]

[TABLE]

and

[TABLE]

where $F(\Delta)$ is depicted in Fig. 2 numerically.

We plot the total energy density versus particle density for different detuning in Fig. 3. As we can see, for $\epsilon_{0}>0$, the minimum energy density is well defined, and lies in the AMSF phase. We expect there exists a self-stabilized state at around the minimum. If the particle number is finite, it forms a quantum droplet Petrov (2015). We also depict the dependence between the particle density of the self-stabilized state $n_{s}$ and the detuning $\epsilon_{0}$ in Fig. 4. It is shown the stabilized density is almost linearly proportional to $\epsilon_{0}$. However, if $\epsilon_{0}<0$, the energy density is degenerated inside ASF phase, but it can be broken by introducing an atom-atom s-wave interaction. Typically, the atom-atom s-wave interaction is repulsive and the corresponding LHY correction is also positive T. D. Lee and Yang (1957). Therefore, the lowest energy density lies at $n=0$ inside ASF phase. For this reason, we do not expect a self-stabilized state when the detuning $\epsilon_{0}<0$.

The diluteness of p-wave interacting gas can be characterized by the product between the particle density and the scattering volume $v_{p}$ Qin (2018); Luciuk et al. (2016), i.e. $nv_{p}=\overline{g}^{2}mn/(16\pi^{2}\epsilon_{0})$. Therefore, we can rewrite the ground-state energy given by mean-field theory (MFT) and the LHY correction in terms of the diluteness as

[TABLE]

[TABLE]

and

[TABLE]

The diluteness of the self-stabilized state with respect to detuning is shown in Fig. 5. As detuning approaches zero, the diluteness tends to diverge, which may indicate that higher order corrections beside MFT and LHY are needed. But for a large detuning regime, the mixture is dilute, so that it is reasonable to characterize our model with only first order beyond-mean-field calculation.

VI Quantum Droplets

According to the above analysis, we find the mean-field collapsing state becomes self-stabilized after considering beyond-mean-field correction. This self-stabilized state forms a quantum droplet when particle number is finite Petrov (2015). To figure out the density distribution of the quantum droplet, we will derive an effective theory to characterize the density profile. Here we employ function $\xi(\textbf{r})$ to characterize the droplet density profile. If the system size is infinite, we have solution $\xi(\textbf{r})=1$ as it should be a uniform gas. However, if the system size is finite, the density profile will be inhomogeneous.

As a qualitative analysis, we will take the local density approximation (LDA). With this approximation, the order parameters can be rewritten as Aybar and Oktel (2019),

[TABLE]

and $\Phi=\sqrt{n_{s,M}}\xi(\textbf{r})\hat{\textbf{z}}$, where

[TABLE]

Here we have chosen $\hat{\textbf{z}}$ direction due to the spontaneous breaking of $SO(3)$ rotation symmetry by $\textbf{Q}=-\overline{g}m\sqrt{n_{s,M}}\xi(\textbf{r})\hat{\textbf{z}}$, where $n_{s,1,2,M}$ are the stabilized densities for the two atomic species and molecule respectively.

To access the analytical form of effective Hamiltonian, an approximative form of $F(\Delta)$ at around the stable point is considered. For $\epsilon_{0}>0$, we find a linearized formula for $F(\Delta)$, which captures its behavior at around the minimum energy density inside the AMSF phase ($0<\Delta<\approx 0.08$) [see Fig. 2]. It is written as

[TABLE]

According to Eq. (31), the approximated total ground-state energy in AMSF phase is given by

[TABLE]

Furthermore, by substituting Eqs. (36) and (37) to Eq. (1), we derive the effective Hamiltonian

[TABLE]

The chemical potential $\tilde{\mu}$ is fixed by the normalization condition $\int d^{3}r|\xi|^{2}=N/n_{s}$, where $N$ is the total number of particle and $n_{s}$ is the stabilized total density. The profile function $\xi(\textbf{r})$ is determined by the GP equation

[TABLE]

which is derived by minimizing the effective Hamiltonian.

The above GP equation is solved numerically by using imaginary time evolution method. The solutions for different detuning and particle numbers are shown in Fig. 6. We can find the quantum droplet is typically suppressed in the $z$ direction. The degrees of suppression decreases for a larger $\epsilon_{0}$. Hence the droplet looks like a pancake when $N/n_{s}$ is large enough but $\epsilon_{0}$ is small (see the upper right subfigure of Fig. 6). We also show the section of the solution where $\overline{g}=-1,m=1,\epsilon_{0}=0.5,N/n_{s}=10^{6}$ in Fig. 7. The density is found to suddenly fall to zero in the horizontal directions ($x$ or $y$ directions), while gently decreasing to zero in the $z$ direction. Except for the boundary regime, the profile varies smoothly everywhere, which implies that LDA could qualitatively catch the features of quantum droplet here. In fact, the anisotropy of quantum droplet arises from the spontaneous breaking of SO(3) rotation symmetry by finite-momentum atomic condensates. It is intrinsically different from the anisotropic quantum droplets in the presence of dipolar interaction Ferrier-Barbut et al. (2016); Kadau et al. (2016); Chomaz et al. (2016); Schmitt et al. (2016) or spin-orbit coupling Cui (2018), where the anisotropy arises from external fields. As we can see, the value on the plateau remains almost constant and close to 1, which will be exactly 1 when the system size goes to infinity. Another special feature of quantum droplet here is that the atomic components carry finite momentums due to the breaking of time-reversal symmetry.

VII Conclusion

In this paper, we study the quantum fluctuation correction to the ground states of a p-wave interacting Bose gas. Beginning with the mean-field analysis of the ground states, it is found that the ground states can be divided into three typical phases for different detunings of molecule channel, i.e. the ASF, AMSF and MSF phases, where particles are condensed into only the atomic, both the molecular and atomic, and only the molecular channels, respectively. Particularly, we find the ground state is unstable in phase AMSF. The unstability of the ground state in the phase AMSF also manifests itself in the emergence of imaginary long-wavelength Bogoliubov excitation modes. Furthermore, we calculate the LHY correction with the Bogoliubov excitations. We find the LHY correction can stabilize the ground state in the mean-field-unstable regime. That means that the p-wave interacting Bose gas is self-stabilized at certain density. Finally, we construct an effective Hamiltonian to characterize the ground state of a finite system. By solving the corresponding GP equation, we find self-stabilized quantum-droplet solutions. Unlike the s-wave case, the quantum droplet is anisotropic and carries finite momentums because the spatial rotation and the time-reversal symmetries are spontaneously broken. Although only the interspecies p-wave interaction is considered here, our results could be extended into the case with weak background s-wave interactions and may be observed in systems like ${}^{85}Rb-^{87}Rb$ Bose mixture Papp et al. (2008); Dong et al. (2016).

Acknowledgements

The authors are indebted to Wei Yi, Chao Gao, Jing-Bo Wang and Fang Qin for helpful discussion. This work is supported by the AFOSR Grant No. FA9550-16-1-0006, the MURI-ARO Grant No. W911NF17-1-0323, the ARO Grant No. W911NF-11-1-0230 (Z. L. and W.V. L.), the National Postdoctoral Program for Innovative Talents of China (Grant No. BX201700156), the National Natural Science Foundation of China (Grant No. 11804221) (J.-S. P.), the Science and Technology Commission of Shanghai Municipality (Grants No.16DZ2260200) and National Natural Science Foundation of China (Grants No.11655002) (J.-S. P. and W.V. L.), and the Overseas Scholar Collaborative Program of NSF of China No. 11429402 sponsored by Peking University (W.V. L.). J.-S. P. also acknowledges the support by Xiong-Jun Liu’s group during his visit to Peking University.