The theoretical description of Bose-Einstein condensates (BECs) is of fundamental importance for capturing their quantum nature [1, 2]. So far, the coherent-state-based Gross-Pitaevskii equation (GPE) [3, 4], along with the perturbative Bogoliubov excitations around ground states [5, 6], has achieved great success when describing the BECs of dilute atomic gases [2, 7]. The bosonic version of Hartree-Fock-Bogoliubov (HFB) formalism [5, 8-11], which provides a self-consistent treatment of the fluctuations, is a natural generalization of the GPE plus Bogoliubov theory. Although the HFB formalism can provide satisfactory results in both zero- and finite-temperature problems [7, 12-14], it fails to yield the gapless excitation spectrum [11, 15], which violates the Hugenholtz-Pines theorem [9]. By adopting a Gaussian-state ansatz (i.e., squeezed-coherent state for a bosonic gas), the Gaussian-state theory (GST) [16], which is equivalent to the HFB theory for the ground-state solutions, takes a completely different route to treat the quantum fluctuations. Specifically, when treating the fluctuation around a ground state, it also incorporates the fluctuations of the correlation functions, which successfully produces the gapless excitation spectrum [17-19]. Physically, the fluctuations of the correlation functions represent two-particle excitations.

As an application of GST to BECs, we considered the simplest scenario by studying the ground states and excitations of a trapped atomic gas with contact interaction [18]. Interestingly, it was found that the ground state of a condensate with weakly attractive interaction becomes a single-mode squeezed vacuum state (SVS), as opposed to the coherent state (CS) under repulsive interaction. Other systems that fit GST studies well are the self-bound quantum droplets in both dipolar and two-component condensates [20-26]. These systems have been widely investigated using the extended GPE [27-30], which is a generalized GPE that includes the contribution of the fluctuation through the perturbative Lee-Huang-Yang correction [31]. Along this direction, we explored the quantum phases of both dipolar and binary droplets in the presence of three-body repulsion [32, 33]. Generally, in the low-density regime, where the interaction is purely attractive, the condensate belongs to the SVS phase, whereas in the high-density regime, where the three-body repulsion dominates, the condensate is in the CS phase. Moreover, for the parameter regime between the SVS and CS phases, the system falls into a single-mode squeezed-coherent state (SCS). Finally, we also studied the collapse dynamics of a trapped single-component condensate using the Gaussian state theory (GST) for mixed states [33], which revealed a distinct type of collapse assisted by fluctuations through amplifying the attractive interaction between atoms.

Technically, to obtain the ground state, one needs to simultaneously evolve, in imaginary time, the normal and anomalous correlation functions of the condensate wave function. This poses a serious challenge for numerical calculations, as it requires tremendous computing resources and even becomes computationally intractable for three-dimensional systems without sufficient symmetries. To reduce the demand for computing resources, we note that a Gaussian state can always be factorized into a multimode coherent-squeezed state [32, 33]. Since, in practice, only a few squeezed modes are highly occupied, we may approximate the Gaussian state by one that involves the coherent mode and these highly occupied squeezed modes. In the present work, we develop the effective Gaussian-state theory (EGST) by deriving the imaginary-time evolution equations for the squeezed modes and the occupation numbers. We then calibrate these effective equations in a Bose gas of self-bound dipolar droplets. We show that the EGST results converge to those obtained via the full GST when the number of squeezed modes included in the calculation is sufficiently large. Particularly, for the SVS and SCS phases of dipolar droplets in which the number of notably occupied squeezed modes is small, the EGST approach has a substantial advantage over the full GST.

The rest of this paper is organized as follows. In section 2, we briefly recall GST. In section 3, we present an alternative approach to factorize the Gaussian state into the form of a multimode squeezed-coherent state. We derive, in section 4, the imaginary-time evolution equations for the highly populated squeezed modes and the corresponding occupation numbers. These equations constitute the effective theory for GST, and their solution yields an approximate ground state very close to that obtained via the full GST. In section 5, we validate the effective Gaussian-state approach by revisiting the system of dipolar droplets. We show that the effective theory improves the efficiency of numerical calculations, in particular when the number of squeezed modes involved in the ground state is small. Finally, we conclude in section 6.

We begin by briefly reviewing GST for bosons [34]. For a many-body quantum system, a pure Gaussian state takes the form where ${\widehat{{\mathscr{U}}}}_{{\rm{GS}}}$ is a unitary operator, θ is a global phase, $\widehat{{\rm{\Psi }}}({\boldsymbol{r}})=\left(\begin{array}{c}\widehat{\psi }({\boldsymbol{r}})\\ {\widehat{\psi }}^{\dagger }({\boldsymbol{r}})\end{array}\right)$ is the bosonic field operator satisfying the commutation relation $[\widehat{{\rm{\Psi }}}({\boldsymbol{r}}),{\widehat{{\rm{\Psi }}}}^{\dagger }({{\boldsymbol{r}}}^{{\prime} })]={{\rm{\Sigma }}}_{z}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })$ [$\equiv {\sigma }_{z}\otimes \delta ({\boldsymbol{r}}-{{\boldsymbol{r}}}^{{\prime} })$ with σ_z being the Pauli matrix], and ${\rm{\Phi }}({\boldsymbol{r}})=\left(\begin{array}{c}\phi ({\boldsymbol{r}})\\ {\phi }^{* }({\boldsymbol{r}})\end{array}\right)$ and ${\rm{\Xi }}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })=\left(\begin{array}{cc}{\xi }_{11}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} }) & {\xi }_{12}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })\\ {\xi }_{12}^{\dagger }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} }) & {\xi }_{11}^{{\rm{T}}}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })\end{array}\right)$ are variational parameters. Particularly, the unitarity of ${\widehat{{\mathscr{U}}}}_{{\rm{GS}}}$ requires that ξ₁₁ and ξ₁₂ are, respectively, Hermitian and symmetric matrices with respect to the indices r and ${{\boldsymbol{r}}}^{{\prime} }$, i.e., ${\xi }_{11}^{{\rm{T}}}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })={\xi }_{11}^{* }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })$ and ${\xi }_{12}^{\dagger }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })={\xi }_{12}^{* }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })$. In equation (1), we have introduced shorthand notations, in which the products in the exponents denote matrix products in the Nambu space and integration in the coordinate space, i.e.,

For convenience, we introduce two unitary operators, i.e., the displacement operator $\widehat{{\mathscr{D}}}({\rm{\Phi }})\equiv {{\rm{e}}}^{{\widehat{{\rm{\Psi }}}}^{\dagger }{{\rm{\Sigma }}}_{z}{\rm{\Phi }}}$ and the squeezing operator $\widehat{{\mathscr{S}}}(\xi )\equiv {{\rm{e}}}^{\frac{{\rm{i}}}{2}{\widehat{{\rm{\Psi }}}}^{\dagger }{\rm{\Xi }}\widehat{{\rm{\Psi }}}}$. One can easily verify that where ${\rm{\Phi }}=\langle \widehat{{\rm{\Psi }}}\rangle $ and ${\mathbb{S}}(\xi )={{\rm{e}}}^{{\rm{i}}{{\rm{\Sigma }}}_{z}{\rm{\Xi }}}$. It is convenient to rewrite ${\mathbb{S}}(\xi )$ as a matrix of the form where ${{\mathbb{S}}}_{21}={{\mathbb{S}}}_{12}^{* }$ and ${{\mathbb{S}}}_{22}={{\mathbb{S}}}_{11}^{* }$. Moreover, from the bosonic commutation relation of $\widehat{{\rm{\Psi }}}$, we have ${\widehat{{\mathscr{S}}}}^{\dagger }(\xi )[\widehat{{\rm{\Psi }}},{\widehat{{\rm{\Psi }}}}^{\dagger }]\widehat{{\mathscr{S}}}(\xi )={{\rm{\Sigma }}}_{z}$, and thereforewhich indicates that ${\mathbb{S}}$ is symplectic and represents a canonical transformation. Alternatively, using the Taylor expansion of ${\mathbb{S}}$, it can be shown that which also gives rise to equation (5).

It is important to note that there is a gauge redundancy in the variational wave function, as different ξ's may lead to the same Gaussian state. In fact, any unitary operator $\widehat{{\mathscr{V}}}$ satisfying $\widehat{{\mathscr{V}}}| 0\rangle =| 0\rangle $ leaves the Gaussian state unchanged, since ${\widehat{{\mathscr{U}}}}_{{\rm{GS}}}\widehat{{\mathscr{V}}}| 0\rangle ={\widehat{{\mathscr{U}}}}_{{\rm{GS}}}| 0\rangle $. Thus, ${\widehat{{\mathscr{U}}}}_{{\rm{GS}}}\widehat{{\mathscr{V}}}$ defines a new Gaussian-state operator. As a concrete example, we consider an operator $\widehat{{\mathscr{V}}}={{\rm{e}}}^{{\rm{i}}{\widehat{\psi }}^{\dagger }{\chi }_{11}\widehat{\psi }}$ with ${\chi }_{11}^{\dagger }={\chi }_{11}$, which, in the Nambu space, can be rewritten as $\widehat{{\mathscr{V}}}={{\rm{e}}}^{\frac{{\rm{i}}}{2}{\widehat{{\rm{\Psi }}}}^{\dagger }{{\rm{\Xi }}}_{0}\widehat{{\rm{\Psi }}}}$, where ${{\rm{\Xi }}}_{0}=\left(\begin{array}{cc}{\chi }_{11} & 0\\ 0 & {\chi }_{11}^{* }\end{array}\right)$. Since there exists ${{\rm{\Xi }}}^{{\prime} }$ such that ${{\rm{e}}}^{\frac{{\rm{i}}}{2}{\widehat{{\rm{\Psi }}}}^{\dagger }{{\rm{\Xi }}}^{{\prime} }\widehat{{\rm{\Psi }}}}={{\rm{e}}}^{\frac{{\rm{i}}}{2}{\widehat{{\rm{\Psi }}}}^{\dagger }{\rm{\Xi }}\widehat{{\rm{\Psi }}}}{{\rm{e}}}^{\frac{{\rm{i}}}{2}{\widehat{{\rm{\Psi }}}}^{\dagger }{{\rm{\Xi }}}_{0}\widehat{{\rm{\Psi }}}}$, the redundancy in Ξ is then confirmed.

To remove this redundancy, we replace the variational parameter Ξ by the covariance matrix ${\rm{\Gamma }}\equiv \langle \{\delta \widehat{{\rm{\Psi }}},\delta {\widehat{{\rm{\Psi }}}}^{\dagger }\}\rangle $, where $\delta \widehat{{\rm{\Psi }}}\equiv \widehat{{\rm{\Psi }}}-{\rm{\Phi }}$ is the fluctuation field and {, } represents the anti-commutator. Γ and Ξ are related as follows: The matrix ${\mathbb{S}}$ has a unitary freedom, as Γ remains unchanged under the transformation ${\mathbb{S}}\to {\mathbb{S}}\left(\begin{array}{cc}{ \mathcal W } & 0\\ 0 & {{ \mathcal W }}^{* }\end{array}\right)$, with ${ \mathcal W }$ being an arbitrary unitary matrix. This unitary freedom is associated with the redundancy of Ξ, which, as shown below, can be removed by requiring that Ξ only contains off-diagonal blocks in Nambu space. It is also worth noting that Γ satisfies

For the ground state of a system described by Hamiltonian $\widehat{H}$, the Gaussian state wave function can be determined by minimizing the energy ${E}_{0}=\langle {\rm{GS}}| \widehat{H}| {\rm{GS}}\rangle $ with respect to Φ and Γ. This is equivalent to evolve a trial Gaussian state ∣GS⟩ in imaginary time τ until it converges. As shown in [34], the imaginary-time evolution of a wave function is governed by the equation Here ${\widehat{H}}_{{\rm{MF}}}$ is the mean-field Hamiltonian, which is obtained by replacing $\widehat{\psi }$ in the Hamiltonian $\widehat{H}$ with $\phi +\delta \widehat{\psi }$, applying Wick's theorem with respect to the Gaussian state, and retaining terms up to the quadratic order in $\delta \widehat{\psi }$. Eventually, ${\widehat{H}}_{{\rm{MF}}}$ can be generally expressed as where ϒ(r) = $\left(\begin{array}{c}\eta ({\boldsymbol{r}})\\ {\eta }^{* }({\boldsymbol{r}})\end{array}\right)$, ${\mathbb{H}}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })=\left(\begin{array}{cc}{ \mathcal E }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} }) & {\rm{\Delta }}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })\\ {{\rm{\Delta }}}^{\dagger }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} }) & {{ \mathcal E }}^{{\rm{T}}}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })\end{array}\right)$, and $:\widehat{O}{:}_{{\rm{GS}}}$ denotes the normal ordered operator with respect to the Gaussian state. Clearly, to find the explicit expressions for η, ${ \mathcal E }$, and Δ, one has to specify the Hamiltonian $\widehat{H}$ as demonstrated in section 5.

To proceed further, we derive, from equation (9), the imaginary-time equations of motion (ITEOM) for the variational parameters Φ and Γ. For this purpose, we first evaluate the left-hand side of equation (9) and find where $:\widehat{O}{:}_{0}$ denotes the normal ordered operator with respect to the vacuum state ∣0⟩. After some straightforward calculations, the right-hand side of equation (9) reduces to For terms on the right-hand side of equations (11) and (12), only those proportional to ${\widehat{\psi }}^{\dagger }$ and ${\widehat{\psi }}^{\dagger }{\widehat{\psi }}^{\dagger }$ give rise to nonzero results after acting on the vacuum state. Therefore, we can find the ITEOM for Φ and Γ by separately matching the coefficients of the terms ${\widehat{\psi }}^{\dagger }$ and ${\widehat{\psi }}^{\dagger }{\widehat{\psi }}^{\dagger }$.

For the ${\widehat{\psi }}^{\dagger }$ term, we have where the subscript '1' ('2') denotes the upper (lower) vector in Nambu space. To establish the relation between the full-sized vectors, we consider the lower halves of the vectors in equation (13). Making use of equation (4), we obtain which, by combining with equation (13), immediately lead to Equations (13) and (15), combining with (6), then give rise to

For EOM of Γ, we match the ${\widehat{\psi }}^{\dagger }{\widehat{\psi }}^{\dagger }$ terms in equations (11) and (12), which gives rise to where the subscript 'ij' denotes the (i, j) block in the Nambu space. Using the result of equation (5), ${{\mathbb{S}}}^{\dagger }{{\rm{\Sigma }}}_{z}{\partial }_{\tau }{\mathbb{S}}$ is anti-Hermitian, which further yields Then for the full matrix in Nambu space, we have where, since the first term on the right-hand side only contains the off-diagonal blocks of ${{\mathbb{S}}}^{\dagger }{{\rm{\Sigma }}}_{z}{\partial }_{\tau }{\mathbb{S}}$, ${\mathbb{Q}}$ is the block diagonal matrix accounting for the diagonal blocks. As a result, ${\mathbb{Q}}$ is also an anti-Hermitian matrix. Now, multiplying equation (19) by Σ_z from the left and by ${{\mathbb{S}}}^{\dagger }$ from the right, we obtain which, by making use of equation (7), eventually leads to Equations (16) and (21) constitute the ITEOM for the Gaussian states [34].

Before moving to the next section, we would like to discuss the physical significance of the steady state, say ${{\mathbb{S}}}_{0}$, obtained from the ITEOM. When ${{\mathbb{S}}}_{0}$ is reached, we have ${\left({{\mathbb{S}}}_{0}^{\dagger }{{\rm{\Sigma }}}_{z}{\partial }_{\tau }{{\mathbb{S}}}_{0}\right)}_{12}=0$, which, according to equation (17), implies that the off-diagonal blocks of ${{\mathbb{S}}}_{0}^{\dagger }{\mathbb{H}}{{\mathbb{S}}}_{0}$ vanish. Furthermore, utilizing equation (4), it can be verified that both diagonal blocks of ${{\mathbb{S}}}_{0}^{\dagger }{\mathbb{H}}{{\mathbb{S}}}_{0}$ are Hermitian and satisfy ${\left({{\mathbb{S}}}_{0}^{\dagger }{\mathbb{H}}{{\mathbb{S}}}_{0}\right)}_{11}={\left({{\mathbb{S}}}_{0}^{\dagger }{\mathbb{H}}{{\mathbb{S}}}_{0}\right)}_{22}^{{\rm{T}}}$. Therefore, ${{\mathbb{S}}}_{0}^{\dagger }{\mathbb{H}}{{\mathbb{S}}}_{0}$ can always be diagonalized by a unitary matrix, say ${\overline{{\mathbb{V}}}}^{{\prime} }\equiv \left(\begin{array}{cc}{\overline{V}}^{{\prime} } & 0\\ 0 & {\overline{V}}^{{\prime} * }\end{array}\right)$, as where σ₀ is the 2 × 2 unit matrix and ${\mathsf{\Omega }}$ is a diagonal matrix formed by the eigenenergies. Due to the unitary freedom of ${{\mathbb{S}}}_{0}$, we can always absorb ${\overline{{\mathbb{V}}}}^{{\prime} }$ into ${{\mathbb{S}}}_{0}$ to form ${{\mathbb{S}}}_{0}^{{\prime} }$ = ${{\mathbb{S}}}_{0}{\overline{{\mathbb{V}}}}^{{\prime} }$. Clearly, ${{\mathbb{S}}}_{0}^{{\prime} }$ gives rise to the same physics as ${{\mathbb{S}}}_{0}$ since ${{\mathbb{S}}}_{0}^{{\prime} }{{\mathbb{S}}}_{0}^{{\prime} \dagger }={{\mathbb{S}}}_{0}{{\mathbb{S}}}_{0}^{\dagger }={{\rm{\Gamma }}}_{0}$. The quadratic term of the mean-field Hamiltonian can now be diagonalized as from which we identify ${{\mathbb{S}}}_{0}^{{\prime} -1}\delta \widehat{{\rm{\Psi }}}$ as the Bogoliubov quasiparticle. Indeed, it can be verified that the Gaussian state ∣GS⟩ is the vacuum state of ${{\mathbb{S}}}_{0}^{{\prime} -1}\delta \widehat{{\rm{\Psi }}}$. Moreover, making use of equation (5), we obtain the Bogoliubov equation from equation (22), i.e., which, in the conventional notation, takes the form Here $\overline{U}=({u}_{1}({\boldsymbol{r}}),{u}_{2}({\boldsymbol{r}}),\cdots \,)$ and $\overline{V}=({v}_{1}({\boldsymbol{r}}),{v}_{2}({\boldsymbol{r}}),\cdots \,)$ are Bogoliubov modes, and ${\mathsf{\Omega }}$ are the energies of the Bogoliubov quasiparticle. Conversely, if ${\mathbb{H}}$ is diagonalized by ${{\mathbb{S}}}_{0}^{{\prime} }$ as equation (22), we then have which implies that ∂_τΓ₀ = 0 and the Γ₀ corresponds to a steady (ground) state.

Previously, it was shown in [32] that a Gaussian state can be factorized into a multimode squeezed-coherent state. Here we present an alternative derivation for this factorization. For this purpose, we rewrite the covariance matrix as where ${ \mathcal G }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })=\langle \delta {\widehat{\psi }}^{\dagger }({{\boldsymbol{r}}}^{{\prime} })\delta \widehat{\psi }({\boldsymbol{r}})\rangle $ and ${ \mathcal F }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })=\langle \delta \widehat{\psi }({\boldsymbol{r}})\delta \widehat{\psi }({{\boldsymbol{r}}}^{{\prime} })\rangle $ are the normal and the anomalous correlation functions, respectively, and ${ \mathcal I }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })=\delta ({\boldsymbol{r}}-{{\boldsymbol{r}}}^{{\prime} })$ is the Dirac function. Apparently, ${ \mathcal G }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })$ is a positive-semidefinite matrix with respect to the indices r and ${{\boldsymbol{r}}}^{{\prime} }$ and ${ \mathcal F }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })$ is a complex symmetric one. As a result, we may diagonalize Γ₁₁ using a unitary matrix $\overline{A}$, i.e., ${\overline{A}}^{\dagger }{{\rm{\Gamma }}}_{11}\overline{A}={{\mathsf{\Gamma }}}_{11}^{(d)}$ with ${{\mathsf{\Gamma }}}_{11}^{(d)}$ being a diagonal matrix with all elements being non-negative. We then rewrite Γ as where ${{\mathsf{\Gamma }}}_{12}^{{\prime} }={\overline{A}}^{\dagger }{{\rm{\Gamma }}}_{12}{\overline{A}}^{* }$ is also a complex symmetric matrix. After substituting equation (28) into the relation equation (8), we find which indicates that the matrix ${{\mathsf{\Gamma }}}_{12}^{{\prime} }$ can be diagonalized simultaneously in the non-degenerate subspace of ${{\mathsf{\Gamma }}}_{11}^{(d)}$. Even in the degenerate subspace of ${{\mathsf{\Gamma }}}_{11}^{(d)}$, the matrix ${{\mathsf{\Gamma }}}_{12}^{{\prime} }$ can be diagonalized by a unitary matrix ${\mathsf{B}}$ as ${{\mathsf{B}}}^{\dagger }{{\mathsf{\Gamma }}}_{12}^{{\prime} }{{\mathsf{B}}}^{* }={{\mathsf{\Gamma }}}_{12}^{(d)}$ according to the Autonne-Takagi diagonalization [35, 36], where ${{\mathsf{\Gamma }}}_{12}^{(d)}$ is a diagonal matrix. We note that the elements of ${{\mathsf{\Gamma }}}_{12}^{(d)}$ can always be chosen to be real and non-negative as the phase of a complex eigenvalue, say ${({{\mathsf{\Gamma }}}_{12}^{(d)})}_{j}=| {({{\mathsf{\Gamma }}}_{12}^{(d)})}_{j}| {{\rm{e}}}^{{\rm{i}}{\varphi }_{j}}$, can be absorbed into the corresponding mode function. With the redefined unitary matrix $W=\overline{A}{\mathsf{B}}$, Γ can be factorized as Namely, Γ₁₁ and Γ₁₂ are now diagonalized simultaneously.

To proceed further, we substitute equation (30) into equation (8) and find which allows us to parameterize them as ${{\mathsf{\Gamma }}}_{11}^{(d)}=\cosh 2{\mathsf{D}}$ and ${{\mathsf{\Gamma }}}_{12}^{(d)}=\sinh 2{\mathsf{D}}$ with ${\mathsf{D}}$ being a diagonal matrix. Here ${\mathsf{I}}$ is the identity of the mode space, i.e., ${{\mathsf{I}}}_{ij}={\delta }_{ij}$. Then, under a suitable gauge, ${\mathbb{S}}$ takes the form where ${\rm{\Xi }}=\left(\begin{array}{cc}0 & -{\rm{i}}\overline{W}{\mathsf{D}}{\overline{W}}^{{\rm{T}}}\\ {\rm{i}}{\overline{W}}^{* }{\mathsf{D}}{\overline{W}}^{\dagger } & 0\end{array}\right)$. Here the redundancy in Ξ has been removed by requiring that Ξ only contains off-diagonal blocks. Furthermore, the covariance matrix can now be parameterized as from which we identify that where $\overline{W}=({\bar{\varphi }}_{1}({\boldsymbol{r}}),{\bar{\varphi }}_{2}({\boldsymbol{r}}),\ldots )$ with column vectors ${\bar{\varphi }}_{j}({\boldsymbol{r}})$ representing the normalized functions for the squeezed modes and ${\mathsf{D}}={\rm{diag}}({d}_{1},{d}_{2},\ldots )$ are squeezing factors. Then, from equation (34), it can be identified that ${N}_{s,j}\equiv {\sinh }^{2}{d}_{j}$ is the particle occupation number in the jth squeezed mode and N_s = ∑_jN_s,j is the total number of squeezed atoms. For convenience, we always assume that d_j and, equivalently, N_s,j are sorted in descending order. We remark that, for the matrix $\overline{W}$, r and j should be regarded as row and column indices, respectively. Furthermore, the unitarity of $\overline{W}$ implies that which, in the form of mode functions, can be expressed, respectively, as and

After removing the redundancy in Ξ, the Gaussian state wave function reduces to where N_c = ∫dr∣φ(r)∣² is the occupation number in the coherent mode, ${\widehat{a}}^{\dagger }={N}_{c}^{-1/2}\int {\rm{d}}{\boldsymbol{r}}\phi ({\boldsymbol{r}}){\widehat{\psi }}^{\dagger }({\boldsymbol{r}})$ and ${\widehat{b}}_{j}^{\dagger }\,=\int {\rm{d}}{\boldsymbol{r}}{\bar{\varphi }}_{j}({\boldsymbol{r}}){\widehat{\psi }}^{\dagger }({\boldsymbol{r}})$ are the creation operators of the coherent and the j-th squeezed modes, respectively.

This factorized form, equation (38), is particularly convenient for characterizing the quantum states of condensates [32, 33]. To see this, let us first introduce the fraction of the coherent atoms, i.e., f_c = N_c/N with N = N_c + N_s being the total number of atoms. The quantum state of a condensate can be characterized as follows. A conventional condensate described by a CS satisfies N_c/N ≈ 1 and N_s/N ≪ 1, for which the squeezed atoms are the quantum depletion. In this case, none of the squeezed modes is significantly occupied. In the opposite limit with N_s,1/N ≈ 1, N_c/N ≈ 0, and N_s,j>1/N ≈ 0, the condensate is in a macroscopic single-mode squeezed vacuum state, i.e., SVS, which, as shown previously [18, 32, 33], possesses completely different statistical properties compared to the coherent state. Finally, in the intermediate case, where both φ and ${\bar{\varphi }}_{s,1}$ are macroscopically occupied, the condensate is in SCS.

In this section, we derive the ITEOM for the squeezed modes ${\bar{\varphi }}_{j}({\boldsymbol{r}})$ and the occupation number N_s,j. Although these equations are equivalent to equation (21) and should lead to the same ground state, their practical advantage appears when the number of occupied squeezed modes is small. In practice, there always exists a truncation S on the number of squeezed modes such that N_s,j are negligibly small for j  >  S. As a result, we only need to find the ITEOM governing the effective modes $\widetilde{W}\equiv ({\bar{\varphi }}_{1},{\bar{\varphi }}_{2},\ldots ,{\bar{\varphi }}_{S})$ and the squeezing factors $\widetilde{{\mathsf{D}}}=\,\rm{diag}\,({d}_{1},{d}_{2},\ldots ,{d}_{S})$. The reduction of the squeezed modes can significantly lower the numerical challenges.

For this purpose, we consider equation (17), which, after making use of equation (32), reduces to First, we extract the ITEOM for the squeezing factors ${\mathsf{D}}$, which are contained in the real part of the diagonal elements of the matrix equation (39). Taking the real part of the diagonal elements with j≤S, we obtain where $\widetilde{{\mathsf{E}}}\equiv {\widetilde{W}}^{\dagger }{ \mathcal E }\widetilde{W}$ and $\widetilde{{\mathsf{\Delta }}}\equiv {\widetilde{W}}^{\dagger }{\rm{\Delta }}{\widetilde{W}}^{* }$ are two S × S matrices, explicitly they are ${\widetilde{{\mathsf{E}}}}_{jj}=\int {\rm{d}}{\boldsymbol{r}}{\rm{d}}{\boldsymbol{r}}^{\prime} {\bar{\phi }}_{j}^{* }({\boldsymbol{r}}){ \mathcal E }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} }){\bar{\phi }}_{j}({{\boldsymbol{r}}}^{{\prime} })$ and ${\widetilde{{\mathsf{\Delta }}}}_{jj}=\int {\rm{d}}{\boldsymbol{r}}{\rm{d}}{{\boldsymbol{r}}}^{{\prime} }{\bar{\phi }}_{j}^{* }({\boldsymbol{r}}){\rm{\Delta }}({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} }){\bar{\phi }}_{j}^{* }({{\boldsymbol{r}}}^{{\prime} })$. We note that, to derive equation (40), we have used the fact that both ${\overline{W}}^{\dagger }{\partial }_{\tau }\overline{W}$ and ${\overline{W}}^{{\rm{T}}}{\partial }_{\tau }{\overline{W}}^{* }$ are anti-Hermitian and their diagonal elements are purely imaginary.

To further derive the ITEOM of the effective squeezed modes, we introduce the projection operator $\widetilde{{\mathsf{P}}}=\left(\begin{array}{cc}\tilde{{\mathsf{I}}} & 0\\ 0 & 0\end{array}\right)$, where $\tilde{{\mathsf{I}}}$ is the S × S unit matrix. Clearly, $\widetilde{{\mathsf{P}}}$ projects a vector onto the subspace formed by the effective modes. The time derivatives of the effective squeezed modes can be expressed as ${\partial }_{\tau }\widetilde{W}={\partial }_{\tau }\overline{W}\widetilde{{\mathsf{P}}}$. In order to utilize equation (39), we multiply ${\partial }_{\tau }\overline{W}\widetilde{{\mathsf{P}}}$ from the left by ${\overline{W}}^{\dagger }$ and then decompose it into After multiplying both sides by $\overline{W}$ from the left, the above equation can be rewritten as where $\widetilde{{\mathsf{M}}}\equiv \widetilde{{\mathsf{P}}}{\overline{W}}^{\dagger }{\partial }_{\tau }\overline{W}\widetilde{{\mathsf{P}}}$ and ${\widetilde{{\mathsf{M}}}}_{\perp }\equiv ({\mathsf{I}}-\widetilde{{\mathsf{P}}}){\overline{W}}^{\dagger }{\partial }_{\tau }\overline{W}\widetilde{{\mathsf{P}}}$ are two block matrices and are schematically shown in figure 1. The matrix elements of $\widetilde{{\mathsf{M}}}$ can be easily read out from equation (39) as In order to extract the matrix elements of ${\widetilde{{\mathsf{M}}}}_{\perp }$, we make use of the fact that d_j = 0 for j  >  S. As a result, equation (39) leads to After plugging it into equation (42), the ITEOM for the effective squeezed modes becomes which, for the k-th squeezed mode, takes the form Finally, making use of equations (16) and (33), one may easily find the ITEOM for the coherent mode Equations (40), (45), and (46) form the working equations for EGST, which we shall refer to as the effective equations in the rest of this article. To solve these equations, we start with a set of initial states, φ(r) and ${\{{d}_{j},{\bar{\varphi }}_{j}({\boldsymbol{r}})\}}_{j=1}^{S}$. We then simultaneously evolve the effective equations until a steady state is obtained. Here, to fix the total number of atoms, the chemical potential is updated accordingly at every time step. Furthermore, the number of squeezed modes S is chosen such that the result converges when S is increased.

We now discuss the computational complexity of the effective equations. In the coordinate representation, if L is the number of spatial grids, φ, $\bar{\varphi },$ and η are all vectors with L elements, and ${ \mathcal G }$, ${ \mathcal F }$, ${ \mathcal E },$ and Δ are L × L matrices. As can be seen, the computational complexities of the effective equations and the full Gaussian-state equations (16) and  (21) are O(S × L) and O(L²), respectively. Therefore, the effective-equation approach can be much more efficient when S is small.

As a demonstration of EGST, here we revisit the self-bound droplets of dipolar condensates. In the second-quantized form, the Hamiltonian of the system takes the form where h₀ = - ℏ²∇²/(2M) - μ with M being the mass of the atoms and μ the chemical potential. The two-body interaction potential is consisting of the contact and the dipolar interactions, with a_s being the s-wave scattering length, μ₀ the vacuum permeability, μ the magnetic dipole moment of the atom, and θ_r the polar angle of r. For convenience, we introduce ϵ_dd = a_dd/a_s to measure the relative strength of the dipolar interaction with respect to the contact interaction, where a_dd = μ₀μ²M/(12πℏ²) is the dipolar length. For the three-body interaction, g₃ represents the three-body coupling constant. After applying Wick's theorem, we obtain [32] where ${ \mathcal G }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })={\sum }_{j=1}^{S}{\sinh }^{2}{d}_{j}{\bar{\varphi }}_{j}^{* }({{\boldsymbol{r}}}^{{\prime} }){\bar{\varphi }}_{j}({\boldsymbol{r}})$ and ${ \mathcal F }({\boldsymbol{r}},{{\boldsymbol{r}}}^{{\prime} })\,={\sum }_{j=1}^{S}\sinh {d}_{j}\cosh {d}_{j}{\bar{\varphi }}_{j}({{\boldsymbol{r}}}^{{\prime} }){\bar{\varphi }}_{j}({\boldsymbol{r}})$. With these quantities specified, one can then numerically evolve the effective equations (40), (45), and (46) to obtain the ground-state wave function. The total energy can then be calculated as E = E_kin + E_2B + E_3B, where and are, respectively, the kinetic, two-body interaction, and three-body interaction energies. In addition, we can further compute the total density n(r) = n_c(r) + n_s(r), where n_c(r) = ∣φ(r)∣² and ${n}_{s}({\boldsymbol{r}})={\sum }_{j=1}^{S}{\sinh }^{2}{d}_{j}| {\bar{\varphi }}_{j}({\boldsymbol{r}}){| }^{2}$ are, respectively, the densities of the coherent and squeezed atoms. We note that all densities are axially symmetric due to the cylindrical symmetry of the dipolar interaction.

As a concrete example, we consider a condensate of N = 10³ dysprosium atoms, each possessing a magnetic dipole moment of μ = 10μ_B, with μ_B being the Bohr magneton. Since the scattering length a_s is tunable via Feshbach resonance, we treat a_s, or equivalently ϵ_dd, as a free parameter. More specifically, we shall focus on the regime with ϵ_dd > 1.23, because the condensate becomes self-bound only in this regime [32]. For the three-body coupling strength, we adopt the value used in Reference [32], i.e., g₃ = 6.73 × 10^-40ℏm⁶/s. Previously, the full GST showed that, by increasing ϵ_dd, the condensate falls into the SVS, SCS, and CS phases sequentially [32]. The phase boundaries are located at ϵ_dd = 1.33 and 1.79. Finally, for convenience, we shall refer to the results obtained from the effective theory as approximate results and those from the full GST as exact results.

To validate EGST, we consider the case with ϵ_dd = 2. In figure 2(a), we plot the S dependence of the total energy E and the coherent atom fraction f_c. Next, in figure 2(b), we plot the S dependence of the condensate radial size ${\sigma }_{\rho }={\left[\int {\rm{d}}{\boldsymbol{r}}{\rho }^{2}n({\boldsymbol{r}})\right]}^{1/2}$ and the peak density n_peak. Here the peak density n_peak is the maximal value of n(r), and we find under all circumstances that n_peak = n(0). Finally, figure 2(c) shows the S dependence of the occupation numbers in the first four squeezed modes. As can be seen, all physical quantities converge to the corresponding exact values when S is roughly over S_conv = 50, which seems quite large. This, however, is not completely unexpected, since, for the given parameter, the system lies deep inside the CS phase which involves many squeezed modes. As shall be shown, the situation can be greatly improved in the SVS and SCS phases.

For a more detailed comparison, in figure 3, we plot the approximate densities of the coherent and squeezed atoms along both radial and axial directions for ϵ_dd = 2 and different values of S. As a comparison, we also plot the exact densities in the corresponding figures. As can be seen, all approximate densities roughly converge to the exact ones for S = 40, which confirms the observation in figure 2.

Finally, to see the convergence in other quantum phases, we plot in figure 4 the error of the total energy, i.e., ΔE = E_approximate - E_exact as a function of ϵ_dd for various S's. As expected, ΔE monotonically decreases with the increase of S throughout the interaction regime. In addition, for a given S, ΔE monotonically increases with ϵ_dd. More specifically, in the SVS phase, it is found that the solutions roughly converge for S = 10. However, across the whole SCS phase, S has to be as large as 40 to ensure the convergence of the solution. The required cutoff becomes even larger in the CS phase, where, as shown in figure 2, S_conv increases to 50. Therefore, substantial efficiency can be gained from the effective equations mainly in the SVS and SCS phases.

In conclusion, we have developed an effective-equation approach for GST by deriving the imaginary-time evolution equations for the squeezed modes and the occupation numbers. We have validated these effective equations using a dipolar droplet system. The results obtained from the effective equations converge to those of the full GST when the number of squeezed modes included in the calculation is sufficiently large, confirming the accuracy of the approach. Since the effective equations involve only the retained squeezed modes, the approach can be numerically more efficient than the full GST. In particular, for the SVS and SCS phases of dipolar droplets, in which the number of notably occupied squeezed modes is small, the effective equation approach has a substantial advantage.