Bosonic particle-correlated states: A nonperturbative treatment beyond mean field

Many useful properties of dilute Bose gases at ultralow temperature are predicted precisely by the (mean-field) product-state Ansatz, in which all particles are in the same quantum state. Yet, in situations where particle-particle correlations become important, the product Ansatz fails. To include correlations nonperturbatively, we consider a new set of states: the particle-correlated state of $N=l\times n$ bosons is derived by symmetrizing the $n$-fold product of an $l$-particle quantum state. Quantum correlations of the $l$-particle state “spread out” to any subset of the $N$ bosons by symmetrization. The particle-correlated states can be simulated efficiently for large $N$, because their parameter spaces, which depend on $l$, do not grow with $n$. Here we formulate and develop in great detail the pure-state case for $l=2$, where the many-body state is constructed from a two-particle pure state. These paired wave functions, which we call pair-correlated states (PCS), were introduced by A. J. Leggett [Rev. Mod. Phys. 73, 307 (2001)] as a particle-number-conserving version of the Bogoliubov approximation. We present an iterative algorithm that solves for the reduced (marginal) density matrices (RDMs), i.e., the correlation functions, associated with PCS in time $O\left(N\right)$. The RDMs can also be derived from the normalization factor of PCS, which is derived analytically in the large-$N$ limit. To test the efficacy of PCS, we analyze the ground state of the two-site Bose-Hubbard model by minimizing the energy of the PCS state, both in its exact form and in its large-$N$ approximate form, and comparing with the exact ground state. For $N=1000$, the relative errors of the ground-state energy for both cases are within ${10}^{-5}$ over the entire parameter region from a single condensate to a Mott insulator. We present numerical results that suggest that PCS might be useful for describing the dynamics in the strongly interacting regime.

Figure 1

The three eigenvalues of the single-particle RDM (normalized to unity) plotted as a function of the number of particles $N$ for the case ${\lambda }_1^2=0.5$, ${\lambda }_2^2=0.3$, ${\lambda }_3^2=0.2$. As $N$ becomes large, the eigenvalue corresponding to the biggest $\lambda$ approaches 1 while the other eigenvalues become negligible.

Figure 2

The three eigenvalues of the 1RDM (normalized to unity) plotted as a function of the number of particles $N$ for the case ${\lambda }_1^2=0.41$, ${\lambda }_2^2=0.39$, ${\lambda }_3^2=0.2$. The solid lines are the eigenvalues calculated by setting ${\lambda }_3=0$ and keeping ${\lambda }_1$ and ${\lambda }_2$ unchanged except for renormalizing; they conform pretty well to the other results for large $N$.

Figure 3

Eigenvalues of the 1RDM (here normalized to 1) as a function of the number of particles. The coefficients ${\lambda }_1$ and ${\lambda }_2$ are fixed; i.e., the parameter $s_{-}$ grows linearly in $N$. The validity of our approximation in the large-$N$ limit is confirmed by the numerically determined eigenvalues.

Figure 4

(a) Population imbalance $\mathrm{\Delta }$ of the Schmidt modes for the exact ground-state solution (black lines, left axis), given by Eq. (7.16), as a function of $U/\left(nJ\right)$ and for number of bosons $N=50$ (dotted lines), 100 (dot-dashed lines), 500 (dashed lines), 1000 (solid lines). Note that by plotting $\mathrm{\Delta }/N$, we remove the trivial dependence on particle number coming from the 1RDM normalization and, as a result, all cases lie on top of each other. This plot evidences the transition between the superfluid/Mott-insulator regimes of the BH model. For weak interactions $U/\left(nJ\right)\ll 1$, the ground state corresponds to a superfluid with a near-unity occupation of a single Schmidt mode ($\mathrm{\Delta }/N\sim {\varrho }_1/N\sim 1$), whose spatial wave function is symmetrically delocalized between the two sites. As interactions increase, so does the depletion of the dominant Schmidt mode; for strong interactions $U/\left(nJ\right)\gg 1$, the superfluid fragments into uncorrelated components on the two sites, with both 1RDM eigenvalues approaching the same order of magnitude (i.e., $\mathrm{\Delta }/N\ll 1$). In the limit $U/\left(nJ\right)\to \infty$, $\mathrm{\Delta }/N\to 0$ and the ground state becomes the Mott-insulator state $|n,n\rangle$, a product state of $n=N/2$ particles at each site. This is demonstrated on the right axis of (a) by the fidelity $|\langle n,n|{\mathrm{\Psi }}_{\mathrm{ex}}\rangle {|}^2$ (orange lines). For reference, we also plot the population imbalance ${\mathrm{\Delta }}_{\mathrm{bog}}$ given by the conventional Bogoliubov approximation for $N=1000$ [thin, green (gray) line, left axis] (see Appendix pp1 for details). As long as interactions remain weak, i.e., $U/\left(nJ\right)\lesssim 1$, the Bogoliubov approximation provides a good description of the condensate depletion. In the transition region from a superfluid to a Mott insulator, $U/\left(nJ\right)\sim 1$, the Bogoliubov approximation breaks down, and as interactions become sufficiently strong and the system transitions to a Mott insulator, the Bogoliubov approximation fails entirely, unable even to proceed past $U/\left(nJ\right)\simeq 64$. (b) Comparison of the exact population imbalance $\mathrm{\Delta }$ (circles) against the PCS and APCS predictions (7.26) and (B3) [blue and light blue (gray and light gray) lines] for $N=50$ (dotted lines), 100 (dot-dashed lines), 500 (dashed lines), 1000 (solid lines). These plots confirm that the ground state is close to a single condensate when $U/\left(nJ\right)\lesssim 1/n$. (c) Relative error between the exact expressions and the PCS and APCS predictions [blue and light blue (gray and light gray) lines] for $N=50,100,500,$ and 1000; plotted for comparison are the Bogoliubov predictions ${\mathrm{\Delta }}_{\mathrm{bog}}$ of Eq. (A12) for the same values of $N$ [thin, green (gray) lines]. The performance of the exact PCS approximation is remarkably good, nearly matching the exact solution over the whole parameter space and for all particle numbers considered. As expected, the error of the exact PCS goes to zero away from the transition region and increases for $U/\left(nJ\right)\sim 1$. In this transition region, the APCS nearly matches the exact PCS, with the accuracy of the large-$N$ approximation increasing with $N$. The fluctuating relative errors of the exact PCS in the Mott-insulator regime come from the numerical minimization of Eq. (B8), whose minimum gets shallower with increasing $U/\left(nJ\right)$. The Bogoliubov approximation breaks down through the transition region and fails entirely in the Mott-insulator regime. For $U/\left(nJ\right)\lesssim {10}^{-1}$, the exact PCS approximation is generally more accurate than the Bogoliubov approximation, whereas APCS is less accurate, reflecting the fact that APCS cannot model precisely a single-mode condensate.

Figure 5

Relative error of the PCS [dark red (dark gray) lines] and APCS [light coral (light gray) lines] approximations to the exact ground-state energy of the two-site BH model as a function of $U/\left(nJ\right)$, for (a) $N=50$ (dotted lines), (b) 100 (dot-dashed lines), (c) 500 (dashed lines), and (d) 1000 (solid lines). Plotted for comparison are the relative errors of the energy of the Bogoliubov ground state [green (gray) lines], calculated using Eq. (A7). The agreement of the (A)PCS approximation is very good over the whole parameter space, with the accuracy of the approximation increasing with $N$. As $U/\left(nJ\right)$ becomes smaller than 1, the accuracy of the APCS energy is degraded, reflecting the inability of the large-$N$ approximation to describe a single-mode condensate. For $U/\left(nJ\right)\gtrsim 1$, the APCS predictions are quite accurate even for $N=50$. The exact PCS prediction is significantly more accurate than the conventional Bogoliubov approximation for all parameters considered, providing a compelling demonstration of the potential of the PCS formalism. Not only does it outperform the Bogoliubov approximation for weak interactions; it also provides a means for an accurate description across the transition into the Mott-insulator phase. The APCS prediction, on the other hand, does not outperform the Bogoliubov approximation in the weak-interaction regime, because the approximations made for the APCS are inconsistent with giving an accurate description of a single condensate mode. Yet, around $U/\left(nJ\right)\sim 1$, APCS becomes more accurate than Bogoliubov and is able to capture successfully the transition to the Mott-insulator phase.

Figure 6

Measures of the error of the PCS (darker colors) and APCS (lighter colors) 2RDMs, ${\rho }_{\mathrm{pcs}}^{\left(2\right)}$ and ${\rho }_{\mathrm{apcs}}^{\left(2\right)}$, relative to the exact 2RDM ${\rho }^{\left(2\right)}$, as quantified by the trace distance (7.32) [blue (upper) lines] and infidelity (7.33) [red (lower) lines] as a function of $U/\left(nJ\right)$, for $N=50$ (dotted lines), 100 (dot-dashed lines), 500 (dashed lines), 1000 (solid lines). Both measures vanish if and only if the 2RDMs being compared are equal. A nonzero value quantifies the overall error (bounded by 1) of the (A)PCS 2RDM. Both error metrics reveal a very good agreement between ${\rho }_{\left(\mathrm{a}\right)\mathrm{pcs}}^{\left(2\right)}$ and ${\rho }^{\left(2\right)}$ over the entire parameter space, from small ($N=50$) to large ($N=1000$) particle numbers. For ${\rho }_{\mathrm{pcs}}^{\left(2\right)}$, the agreement is nearly exact in the limits $U/\left(nJ\right)\to 0$ and $U/\left(nJ\right)\to \infty$. The PCS accuracy is worse near the phase transition, $U/\left(nJ\right)\sim 1$, where the trace distance (infidelity) shows an imprecision of the order of ${10}^{-2}$ (${10}^{-3}$). The agreement between ${\rho }_{\mathrm{apcs}}^{\left(2\right)}$ and ${\rho }^{\left(2\right)}$ is nearly as good near the transition, but saturates to a nonzero value for $U/\left(nJ\right)\ll 1$ and $U/\left(nJ\right)\gg 1$. Not surprisingly, the accuracy of the APCS approximation increases with increasing particle number.

Figure 7

(a) and (b) Comparison of the PCS (darker colors) and APCS (lighter colors) 2RDM parameters ${\gamma }_{\left(\mathrm{a}\right)\mathrm{pcs}}$ and ${\delta }_{\left(\mathrm{a}\right)\mathrm{pcs}}$ to the exact ones $\gamma$ and $\delta$ [see Eqs. (7.22) and (7.22)], with their respective relative errors (c) and (d), as a function of $U/\left(nJ\right)$ for $N=50$ (dotted), 100 (dot-dash), 500 (dashed), 1000 solid. The behavior here is generally consistent with that of the error measures plotted in Fig. 6; the exact PCS predictions well approximate the exact 2RDM parameters, with both lines lying on top of each other.

Figure 8

Time evolution of the two-site Bose-Hubbard model of $N=1000$ atoms. Initially, all the atoms sit in the ground state of the noninteracting Hamiltonian, and then an on-site interaction with strength $U=0.8J$ is suddenly turned on. The plots show the trace distances (normalized to 1) of the exact 2RDM to the closest 2RDMs given by PCS, the double-Fock state (DFS), and the Gross-Pitaevskii state (GPS) as functions of evolution time (in units of $\hslash /J$).