Bosonic Particle-Correlated States: A Nonperturbative Treatment Beyond Mean Field

We consider a natural generalization of the product ansatz for Bose-Einstein condensates; the particle-correlated state of $N=l\times n$ identical particles 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$ particles 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$. We pay special attention to the pure-state case for $l=2$, where the many-body state is constructed from a two-particle pure state. These paired wave functions were introduced by Leggett [Rev. Mod. Phys. $\mathbf{73}$, 307 (2001)] as a particle-number-conserving version of the Bogoliubov approximation. For large $N$, we derive few-particle reduced density matrices (correlation functions) for these wave functions. To test the efficacy of our theory, we solve the two-site Bose-Hubbard model by minimizing the energy using the two-particle reduced density matrices that we derived analytically. We find that the relative errors of the ground state energy are within $10^{-5}$ for $N=1000$ particles over the entire range from a single condensate to a Mott insulator.