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

TL;DR

This paper introduces particle-correlated states (PCS) for dilute Bose gases, providing a nonperturbative method that captures correlations beyond mean-field theory and efficiently models large systems.

Contribution

The paper develops a detailed formulation and efficient algorithms for pair-correlated states (PCS), extending beyond mean-field approximations to include particle correlations nonperturbatively.

Findings

PCS accurately predicts ground state energies within 10^{-5} for N=1000 bosons.

The iterative algorithm for RDMs operates in linear time with respect to N.

Numerical results indicate PCS's potential for modeling strongly interacting Bose gases.

Abstract

Many useful properties of dilute Bose gases at ultra-low 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. 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. ${\bf 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(N)$. 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=1\,000$, 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.

Equations

$$\displaystyle\big{\lvert}\,\varPsi_{\mathrm{Fock}}\,\big{\rangle}=\prod_{j=1}^{\nu}\frac{1}{\sqrt{N_{j}!}}\big{(}a^{\dagger}_{j}\big{)}^{N_{j}}\big{\lvert}\,\mathrm{vac}\,\big{\rangle}\;,$$

$${\Psi }_{\mathrm{Jast}}({\mathbf{x}}_1,\dots ,{\mathbf{x}}_N)\sim \prod _{\begin{array}{c}j,k=1\\ j<k\end{array}}^{N}f({\mathbf{x}}_j-{\mathbf{x}}_k),$$

$$\displaystyle\rho\hskip 0.4pt\big{(}\mathbf{x}^{(l)}_{1},\ldots,\mathbf{x}^{(l)}_{n}\,;\,\mathbf{y}^{(l)}_{1},\ldots,\mathbf{y}^{(l)}_{n}\big{)}=\frac{\mathcal{P}_{S}\,\sigma(\mathbf{x}^{(l)}_{1},\,\mathbf{y}^{(l)}_{1})\otimes\cdots\otimes\sigma(\mathbf{x}^{(l)}_{n},\,\mathbf{y}^{(l)}_{n})\,\mathcal{P}_{S}}{\operatorname{tr}\!\Big{(}\mathcal{P}_{S}\,\sigma(\mathbf{x}^{(l)}_{1},\,\mathbf{y}^{(l)}_{1})\otimes\cdots\otimes\sigma(\mathbf{x}^{(l)}_{n},\,\mathbf{y}^{(l)}_{n})\,\mathcal{P}_{S}\Big{)}}\;.$$

$$\displaystyle\mathcal{P}_{S}\,\big{\lvert}\,\psi_{1},\ldots,\psi_{N}\,\big{\rangle}=\frac{1}{N!}\sum_{\pi\in S_{N}}\,\big{\lvert}\,\psi_{\pi(1)},\ldots,\psi_{\pi(N)}\,\big{\rangle}\;,$$

$${\Psi }_{\mathrm{pcs}}$$

$${\Psi }^{(2)}({\mathbf{x}}_1,{\mathbf{x}}_2)=\sum _{j=1}^{\nu }{\lambda }_j{\psi }_j({\mathbf{x}}_1){\psi }_j({\mathbf{x}}_2),$$

$$\sum _{j=1}^{\nu }{\lambda }_j^2=1.$$

$${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_{\nu }.$$

$$\displaystyle\mathcal{A}^{\dagger}\equiv\int\varPsi^{(2)}(\mathbf{x}_{1},\,\mathbf{x}_{2})\uppsi^{\dagger}(\mathbf{x}_{1})\uppsi^{\dagger}(\mathbf{x}_{2})\,d\mathbf{x}_{1}\,d\mathbf{x}_{2}=\sum_{j=1}^{\nu}\lambda_{j}\hskip 0.4pt\big{(}a^{\dagger}_{j}\big{)}^{2}\;.$$

$$a_j^{†}=\int {\psi }_j(\mathbf{x}){\text{uppsi}}^{†}(\mathbf{x})d\mathbf{x}$$

$$\lvert\,\varPsi^{(2)}\,\rangle=\frac{1}{\sqrt{2}}\,\mathcal{A}^{\dagger}\,\lvert\,\mathrm{vac}\,\rangle=\frac{1}{\sqrt{2}}\,\sum_{j=1}^{\nu}\lambda_{j}\hskip 0.4pt\big{(}a^{\dagger}_{j}\big{)}^{2}\,\lvert\,\mathrm{vac}\,\rangle\;,$$

$$|{\Psi }^{(2)}⟩=\frac{1}{\sqrt{2}}\sum _{j,k=1}^{\mu }{\Lambda }_{jk}{b}_k^{†}{b}_j^{†}|\mathrm{vac}⟩.$$

$$U\Lambda\,U^{T}=\operatorname{diag}\big{(}\lambda_{1},\ldots,\lambda_{\mu}\big{)}\;.$$

$$\lvert\,\varPsi^{(2)}\,\rangle=\frac{1}{\sqrt{2}}\,\sum_{j=1}^{\nu}\lambda_{j}\hskip 0.4pt\big{(}a^{\dagger}_{j}\big{)}^{2}\,\lvert\,\mathrm{vac}\,\rangle\;.$$

$$\displaystyle\begin{split}\Psi^{\dagger}(\mathbf{x}_{1})\cdots\Psi^{\dagger}(\mathbf{x}_{N})\,\big{\lvert}\,\mathrm{vac}\,\big{\rangle}&=\sqrt{N!}\;\mathcal{P}_{S}\big{\lvert}\,\mathbf{x}_{1},\ldots,\mathbf{x}_{N}\,\big{\rangle}=\frac{1}{\sqrt{N!}}\sum_{\pi\in S_{N}}\big{\lvert}\,\mathbf{x}_{\pi(1)},\ldots,\mathbf{x}_{\pi(N)}\,\big{\rangle}\;,\end{split}$$

$$\displaystyle\big{(}\mathcal{A}^{\dagger}\big{)}^{n}\,\big{\lvert}\,\mathrm{vac}\,\big{\rangle}$$

$$\displaystyle\mathcal{P}_{S}\Bigl{(}\varPsi^{(2)}(\mathbf{x}_{1},\,\mathbf{x}_{2})\cdots\varPsi^{(2)}(\mathbf{x}_{2n-1},\,\mathbf{x}_{2n})\Bigr{)}$$

$$\frac{1}{\sqrt{N!}}\,\big{\langle}\,\mathbf{x}_{1},\ldots,\mathbf{x}_{2n}\,\big{\rvert}\big{(}\mathcal{A}^{\dagger}\big{)}^{n}\big{\lvert}\,\mathrm{vac}\,\big{\rangle}=\mathcal{P}_{S}\Bigl{(}\varPsi^{(2)}(\mathbf{x}_{1},\,\mathbf{x}_{2})\cdots\varPsi^{(2)}(\mathbf{x}_{2n-1},\,\mathbf{x}_{2n})\Bigr{)}\;,$$

$$\displaystyle\lvert\,\varPsi_{\mathrm{pcs}}\,\rangle\equiv\frac{1}{\sqrt{\mathsf{N}}}\,\big{(}\mathcal{A}^{\dagger}\big{)}^{n}\,\lvert\,\mathrm{vac}\,\rangle=\frac{1}{\sqrt{\mathsf{N}}}\bigg{(}\sum_{j=1}^{\nu}\lambda_{j}\hskip 0.4pt\big{(}a^{\dagger}_{j}\big{)}^{\!2}\bigg{)}^{n}\,\lvert\,\mathrm{vac}\,\rangle\;,$$

$$\displaystyle\varPsi_{\mathrm{pcs}}\propto\,\sum_{j}\lambda_{j}\,\psi_{j}(\mathbf{x}_{1})\,\mathcal{P}_{S}\Bigl{(}\psi_{j}(\mathbf{x}_{2})\,\varPsi^{(2)}(\mathbf{x}_{3},\,\mathbf{x}_{4})\cdots\varPsi^{(2)}(\mathbf{x}_{2n-1},\,\mathbf{x}_{2n})\Bigr{)}\;.$$

$$\displaystyle\begin{split}\varPsi_{\mathrm{pcs}}&\,\propto\,\varPsi^{(2)}(\mathbf{x}_{1},\,\mathbf{x}_{2})\,\mathcal{P}_{S}\Bigl{(}\varPsi^{(2)}(\mathbf{x}_{3},\,\mathbf{x}_{4})\cdots\varPsi^{(2)}(\mathbf{x}_{2n-1},\,\mathbf{x}_{2n})\Bigr{)}\\ &\qquad+(N-2)\sum_{j,k}\lambda_{\smash{j}}\lambda_{\smash{k}}\,\psi_{j}(\mathbf{x}_{1})\psi_{k}(\mathbf{x}_{2})\,\mathcal{P}_{S}\Bigl{(}\psi_{j}(\mathbf{x}_{3})\psi_{k}(\mathbf{x}_{4})\,\varPsi^{(2)}(\mathbf{x}_{5},\,\mathbf{x}_{6})\cdots\varPsi^{(2)}(\mathbf{x}_{2n-1},\,\mathbf{x}_{2n})\Bigr{)}\;.\end{split}$$

$$\displaystyle\big{\lvert}\,\varPsi_{\mathrm{VB}}\,\big{\rangle}=\frac{1}{\sqrt{\mathcal{M}}}\,\exp\!\Big{(}\,\lambda_{0}\hskip 0.4pta_{0}^{\dagger}\hskip 0.4pta_{0}^{\dagger}+2\sum_{k}\lambda_{k}\hskip 0.4ptb^{\dagger}_{k}b^{\dagger}_{-k}\Big{)}\,\big{\lvert}\,\mathrm{vac}\,\big{\rangle}\;,$$

$$\displaystyle\begin{split}\big{\lvert}\,\varPsi_{\mathrm{VB}}\,\big{\rangle}&=\frac{1}{\sqrt{\mathcal{M}}}\,e^{\lambda_{0}\hskip 0.4pta_{0}^{\dagger}\hskip 0.4pta_{0}^{\dagger}}\,\big{\lvert}\,\mathrm{vac}\,\big{\rangle}\\ &=\frac{1}{\sqrt{\mathcal{M}}}\,\sum_{n=0}^{\infty}\frac{\lambda_{0}^{n}}{n!}\,\sqrt{(2n)!}\;\big{\lvert}\,2n\,\big{\rangle}_{0}\otimes\big{\lvert}\,\mathrm{vac}\,\big{\rangle}_{\perp}\\ &\simeq\frac{1}{\sqrt{\mathcal{M}}}\,\sum_{n=0}^{\infty}\frac{\;(2\lambda_{0})^{n}}{\sqrt[4]{\pi n}}\>\big{\lvert}\,2n\,\big{\rangle}_{0}\otimes\big{\lvert}\,\mathrm{vac}\,\big{\rangle}_{\perp}\;,\end{split}$$

$$\displaystyle\big{\lvert}\,\varPsi_{\mathrm{Legg}}\,\big{\rangle}\propto\Big{(}\,a_{0}^{\dagger}\hskip 0.4pta_{0}^{\dagger}+2\sum_{k}\lambda_{k}b^{\dagger}_{k}b^{\dagger}_{-k}\Big{)}^{\!N/2}\,\big{\lvert}\,\mathrm{vac}\,\big{\rangle}\;,$$

$$a_{k0}^{†}=\frac{1}{\sqrt{2}}(b_k^{†}+b_{-k}^{†}),a_{k1}^{†}=-\frac{i}{\sqrt{2}}(b_k^{†}-b_{-k}^{†}),$$

$$\displaystyle\lvert\,\varPsi_{\mathrm{pcs}}\,\rangle\propto\Big{(}\,a_{0}^{\dagger}a_{0}^{\dagger}+\sum_{m>0}^{\nu-1}\lambda_{m}a^{\dagger}_{m}a^{\dagger}_{m}\Big{)}^{N/2}\,\lvert\,\mathrm{vac}\,\rangle\;,$$

$$|{\Psi }_{\mathrm{gpgs}}⟩$$

$$\tilde{b}_{j}^{\dagger}=b_{j}^{\dagger}a_{0}/\sqrt{N}\quad\mbox{for $j=1,2,\ldots,\nu-1$,}$$

$$\displaystyle\mathcal{H}_{\mathrm{ncb}}=\sum_{j,k=1}^{\nu-1}H_{jk}\tilde{b}_{j}^{\dagger}\tilde{b}_{k}+\frac{1}{2}\Big{(}H_{jk}^{\prime}\tilde{b}_{j}^{\dagger}\tilde{b}_{k}^{\dagger}+\mathrm{H.c.}\Big{)}\;,$$

$${\mathcal{B}}^{†}{\mathcal{H}}_{\mathrm{ncb}}\mathcal{B}=\sum _{k=1}^{\nu -1}{ϵ}_k{\tilde{b}}_k^{†}{\tilde{b}}_k,$$

Full text

I Introduction

The mean-field Gross-Pitaevskii (GP) theory Gross (1961); Pitaevskii (1961) successfully predicts many useful properties of weakly interacting Bose gases at ultra-low temperatures, yet a number of interesting many-particle phenomena, such as quantum phase transitions from superfluids to Mott insulators Jaksch et al. (1998); Greiner et al. (2002a), cannot be explained by the mean-field GP approach. This is because particle-particle correlations are neglected when one approximates the state of the system by the GP Ansatz, i.e., by a product of identical single-particle states at all times. The Bogoliubov approximation Bogoliubov (1947); Fetter (1972); Huang (1987) is an attempt to include particle correlations by perturbing the GP mean-field solution with collective excitations. Although useful for analyzing the stability and validity of the GP solution, the Bogoliubov approximation is only valid when depletion of the condensate mode is quite small, i.e., when the single-particle reduced density matrix has only one dominant eigenvalue Penrose and Onsager (1956).

To tackle the case where the condensate is fragmented into two parts, many authors have adopted the double-Fock state (sometimes also called a twin-Fock state); more generally, for a condensate fragmented into many parts, one can use the many-Fock state,

[TABLE]

where $\nu$ is the number of fragments, $\lvert\,\mathrm{vac}\,\rangle$ is the vacuum state, and $a_{j}^{\dagger}$ creates a particle in the single-particle state for the $j$th fragment. This approach thus uses a product-state Ansatz for each of the fragments. Using this Ansatz, Streltsov et al. Streltsov et al. (2004) argued that fragmentation of the ground state of a condensate only happens when the total number of particles is finite; Mueller et al. Mueller et al. (2006) showed that as degeneracies multiply, so do the varieties of fragmentation; and Alon et al. Alon et al. (2007) generalized the GP equation to include the several single-particle wave functions for the various fragments. In contrast, fragmented condensates are also treated by evolving the single-particle reduced density matrix in an approximation that includes the back-reaction from Bogoliubov excitations Vardi and Anglin (2001); Tikhonenkov et al. (2007).

A powerful idea for including particle correlations in many-particle bosonic or fermionic wave functions is to construct them from two-particle states. For instance, the Jastrow wave function of $N$ particles is a product of two-particle states of all $N(N-1)/2$ pairs,

[TABLE]

with bosons (fermions) corresponding to $f$ having even (odd) parity. Many famous wave functions are of Jastrow type, e.g., the Laughlin wave function Laughlin (1983) and the Gutzwiller wave function Gutzwiller (1963); Rokhsar and Kotliar (1991). The Jastrow wave function has found wide application in strongly interacting systems: It is used in variational quantum Monte Carlo as a trial wave function Umrigar et al. (1988); it is used to show that the single-particle reduced (marginal) density matrix of ${}^{4}\mathrm{He}$ is an extensive quantity, thus demonstrating that Bose-Einstein condensation (BEC) underlies superfluidity Reatto (1969); and it is used to investigate the effect of the interatomic correlations and the accuracy of the GP equation Fabrocini and Polls (1999); DuBois and Glyde (2001); Cowell et al. (2002). The Jastrow wave function, however, suffers from increasing demand for computational power for large numbers of bosons and from the requirement of using quantum Monte Carlo.

Here we propose a new set of states, which constitute a natural generalization of the GP product-state Ansatz,

[TABLE]

Here \mathbf{x}^{(l)}_{j}=\big{(}\mathbf{x}_{j,1},\ldots,\mathbf{x}_{j,l}\big{)} and \mathbf{y}^{(l)}_{j}=\big{(}\mathbf{y}_{j,1},\ldots,\mathbf{y}_{j,l}\big{)} denote the coördinates of blocks of $l$ particles, $\sigma$ is an arbitrary state (density matrix) of $l$ particles, and $\mathcal{P}_{S}$ is the projection operator onto the symmetric subspace of all the $N=l\times n$ particles,

[TABLE]

where the sum is over the permutations $\pi$ in the symmetric group $S_{N}$. The state (3) is derived by symmetrizing the $n$-fold tensor product of the $l$-particle state $\sigma$; we call the resulting state a bosonic particle-correlated state (BPCS). The $l$-particle states $\sigma$ can be restricted to symmetrized states, or they can be left arbitrary, with $\mathcal{P}_{S}$ taking care of the symmetrization when the BPCS state is constructed. Note that we can extend (3) to the case of fermions simply by substituting the anti-symmetrizing operator $\mathcal{P}_{A}$ for the symmetrizing operator $\mathcal{P}_{S}$. The resulting state can be called a fermionic particle-correlated state (FPCS). As a consequence of symmetrization, the quantum correlations existing in the $l$-boson state $\sigma$ “spread out” to any subset of the $N=l\times n$ bosons. Moreover, the parameter space of the BPCS does not grow with $n$; it remains the same as that of the bosonic states for $l$ particles.

This article is devoted to the case that $\sigma$ is pure and $l=2$, which we refer to simply as PCS, where now PCS can be read as pair-correlated state. Despite being constructed from a two-particle wave function, PCS is different from Jastrow’s wave function.

This article is organized as follows. In Sec. II, we express PCS, i.e., the case that $\sigma$ is pure and $l=2$, in second-quantized form, and we discuss its introduction by Leggett Leggett (2001) to implement a particle-number-conserving version of the Bogoliubov approximation. In Sec. III, we prove that PCS indeed reproduces the number-conserving Bogoliubov approximation as a special case and, hence, encompasses weakly interacting Bose gases with small condensate depletion. In Sec. IV, we show how to calculate the diagonal elements of reduced density matrices (RDMs) exactly from the PCS normalization factor. In Sec. V, we derive approximations for the RDMs that are valid in the case of a large number of bosons; we discuss exactly solvable examples of these approximations in Sec. V.2. In Sec. VI, we show that the off-diagonal elements of the RDMs can be determined exactly from the diagonal elements. In particular, we show in Sec. VI.1 how to calculate exactly all the matrix elements of the two-particle RDM (2RDM), we give the large-$N$ limit of the 2RDM for two Schmidt orbitals in terms of modified Bessel functions in Sec. VI.2, and in Sec. VI.3, we find the off-diagonal elements of all RDMs in the large-$N$ limit. The analytical calculation of few-particle RDMs of PCS allows the direct evaluation of physical observables in various regimes. An appealing feature of PCS is that they can represent quantum states with or without off-diagonal long-range order Yang (1962). For example, both the superfluid and the Mott-insulating phases in the two-site Bose-Hubbard model can be described by the PCS Ansatz. An interesting question is whether the PCS Ansatz can faithfully interpolate between the two phases. We show in Sec. VII that the answer is yes: The PCS description, both in its exact form and in its large-$N$ approximate form, provides a remarkably good account of the ground state of the two-site Bose-Hubbard model across the entire parameter region from superfluid to Mott insulator. Moreover, we present numerical results that suggest that PCS might be useful for describing the dynamics in the strongly interacting regime. We conclude in Sec. VIII.

This paper is mainly devoted to developing the mathematical formalism for manipulating and using PCS and, in particular, to calculating the reduced (marginal) density matrices (RDMs)—these are the correlation functions—associated with PCS and to investigating the PCS ground state of the two-site Bose-Hubbard model. The results we present at the end of Sec. VII on the dynamics of PCS for the two-site Bose-Hubbard model assume, as is always true, that the single-particle spatial wave functions (orbitals) are held fixed. A subsequent article will explore the power of PCS dynamics when the spatial mode functions are included in the dynamics Jiang and Caves . This investigation will be based on deriving time-dependent equations for the PCS Ansatz by evolving a state (initially in PCS form) for an infinitesimal time and then projecting it back to the PCS manifold. In subsequent papers, we also plan to discuss the PCS ground state for more general problems, e.g., the Bose-Hubbard model for more than two sites and fragmented spin-orbital coupled BECs in a trapping potential.

The material in this paper is taken mainly from Jiang’s doctoral dissertation at the University of New Mexico Jiang (2014).

II Pair-Correlated States

II.1 Second quantization of pair-correlated states

It is difficult to do any calculation with the form (3), because of the need for an explicit symmetrization. This motivates going to a second-quantized picture, where the symmetrization is taken care of automatically. For a pure PCS with $l=2$, the PCS is specified by a two-boson wave function $\varPsi^{(2)}(\mathbf{x}_{1},\,\mathbf{x}_{2})$, and the PCS wave function is given by

[TABLE]

Such a PCS can be regarded as constructed by a mapping of the two-boson Hilbert space into a submanifold of the Hilbert space of $2n$ bosons. Note that the two-boson wave function always has a Schmidt decomposition of the form Paškauskas and You (2001)

[TABLE]

where the single-particle wave functions $\{\psi_{j}(\mathbf{x})\,|\,j=1,2,\ldots,\nu\}$ form orthonormal Schmidt bases for the two particles, $\nu$ is the Schmidt rank, and the (real and positive) $\lambda_{j}$s are the Schmidt coefficients. They satisfy the normalization condition

[TABLE]

The coincidence of the Schmidt bases of the two bosons is a consequence of the symmetry of the wave function. Throughout this paper we order the $\lambda_{j}$s as

[TABLE]

Note also that we use the terms single-particle states, modes, and orbitals interchangeably throughout the paper, both for the single-particle states and for their wave functions.

To get into the second-quantized picture, we find it convenient to introduce a pair creation operator,

[TABLE]

Here $\uppsi^{\dagger}(\mathbf{x})$ is the field operator that creates a particle at $\mathbf{x}$, and

[TABLE]

is the operator that creates a particle in the single-particle state with wave function $\psi_{j}(\mathbf{x})$. The two-particle state,

[TABLE]

has the wave function $\langle\,\mathbf{x}_{1},\mathbf{x}_{2}\,|\,\varPsi^{(2)}\,\rangle=\varPsi^{(2)}(\mathbf{x}_{1},\,\mathbf{x}_{2})$ of Eq. (6).

It is instructive to derive the Schmidt decomposition directly in the second-quantized picture, where an arbitrary two-boson state is given by

[TABLE]

Here $b^{\dagger}_{j}$ creates a particle in the $j$th single-particle state, and $\Lambda_{jk}=\Lambda_{kj}$ is a symmetric matrix [if $\Lambda$ is not symmetric, we can always make it so by redefining $\Lambda\rightarrow(\Lambda+\Lambda^{T})/2$, without changing $\lvert\,\varPsi^{(2)}\,\rangle$], satisfying $\mbox{tr}(\Lambda\Lambda^{\dagger})=1$ to make \big{\lvert}\,\varPsi^{(2)}\,\big{\rangle} normalized. The Autonne-Takagi factorization theorem (see Corollary 4.4.4(c) of Horn and Johnson Horn and Johnson (2013)) says that any complex symmetric matrix $\Lambda$ can be diagonalized by a unitary matrix $U$:

[TABLE]

The $\lambda_{j}$s are real and nonnegative and normalized according to Eq. (7); they are the singular values of $\Lambda$, and the diagonalization (13) is a special case of the singular-value decomposition, specialized to symmetric matrices. By introducing a new set of creation operators, $a_{j}^{\dagger}=\sum_{k=1}^{\nu}U_{jk}^{*}\,b_{k}^{\dagger}$, we bring the state (12) into Schmidt form:

[TABLE]

The number of nonzero singular values is the rank $\nu\leq\mu$ of $\Lambda$. We only need the nonzero singular values in Eq. (11), and the sum is written so as to recognize this.

We are now prepared to give the second-quantized version of the PCS. We begin with

[TABLE]

where the sum is over the permutations $\pi$ in the symmetric group $S_{N}$. This gives immediately that

[TABLE]

where

[TABLE]

Equation (16) is the relation between the first- and second-quantized pictures; it can be written in the equivalent form

[TABLE]

which states that one gets the PCS state (5) by applying the pair creation operator $n$ times to the vacuum.

The state in Eq. (18) is not normalized. The properly normalized PCS state for $N=2n$ particles is given by

[TABLE]

where $\mathsf{N}$ is a normalization factor that plays an important role in our consideration of PCS. The PCS and this normalization factor are functions of $n$ (or $N$) and the $\lambda_{j}$s, but for brevity, we generally allow this functional dependence to remain implicit; when we want to be explicit, we denote the normalization factor as $\mathsf{N}_{\smash{\vec{\lambda}},n}$. Henceforth, we abandon the normalization (7) of the Schmidt coefficients, requiring only that the $\lambda_{j}$s be real and positive; we can do this because an overall scaling of the $\lambda_{j}$s is automatically absorbed into $\mathsf{N}$.

The second-quantized form (19) is convenient for calculations, but we can build some intuition by considering relative-state decompositions in the position basis of the first-quantized form (5) of the PCS. The relative-state decomposition of the particle $\mathbf{x}_{1}$, relative to all the other particles, is

[TABLE]

This relative-state decomposition reveals a Schmidt decomposition, and the Schmidt basis of the particle $\mathbf{x}_{1}$ consists of all the single-particle wave functions $\psi_{j}(\mathbf{x}_{1})$. This means that the single-particle reduced density matrix (1RDM) is diagonal in the basis of the wave functions $\psi_{j}(\mathbf{x}_{1})$. The Schmidt coefficients of the decomposition (20), i.e., the square roots of the eigenvalues of the 1RDM, are not, however, given by the $\lambda_{j}$s, because the norms of the relative states of the other particles are different for different values of $j$.

More interestingly, we have the following relative-state decomposition for particles $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$,

[TABLE]

What this shows is that the two-particle reduced density matrix (2RDM) comes from two sorts of terms: In the first term, the two particles $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ can be perfectly correlated, whereas in the second, they are only partially correlated. To determine the pairwise quantum correlations in the PCS, we need to find the 2RDM, and to do that, it turns out to be useful to investigate thoroughly the functional dependence of the normalization factor $\mathsf{N}$, which plays a role akin to a partition function.

Before getting to that in Sec. IV, however, we detour into some historical remarks in Sec. II.2 and into showing how the Bogoliubov approximation arises from the PCS in Sec. III.

II.2 Pairing wave functions

Inspired by the BCS wave function proposed by Bardeen, Cooper, and Schrieffer for superconductivity Bardeen et al. (1957a, b), Valatin and Butler Valatin and Butler (1958) introduced a similar pairing wave function for bosons,

[TABLE]

where $\mathcal{M}$ is a normalization factor. This state, with a quadratic form of the creation operators in the exponential, is very different from a coherent state. A coherent state has a Poissonian number distribution, peaked around $N=\lvert\alpha\rvert^{2}$, whereas the Valatin-Butler state (22) satisfies an exponential number distribution. Consequently, it might not be suited to situations where the total number of particles is conserved. An easy way to see the exponential number distribution is by setting $\lambda_{k}=0$ for all $k\neq 0$, which gives

[TABLE]

where the subscript $\perp$ denotes the modes labeled by $k\neq 0$, which are orthogonal to the zero mode, and we use Stirling’s formula in the last step. It turns out this exponential distribution is valid even when more than one $\lambda_{k}$ is nonzero. Unlike the pairing wave function for fermions, which is always normalizable, the Valatin-Butler wave function (22) can only be normalized when $\lvert\lambda_{k}\rvert<1/2$ for all $k\in\{0,1,\ldots\}$. The Valatin-Butler wave function has been used to investigate the transition from a single condensate to a multimode condensate Coniglio and Marinaro (1967); Evans and Imry (1969); Noziéres and Saint James (1982).

A number-conserving version of the Valatin-Butler wave function,

[TABLE]

was introduced by Leggett Leggett (2001, 2003). Leggett’s state is the $n$th term, where $n=N/2$, in the expansion of the exponential in the Valatin-Butler state. It is a special case of the BPCS with $\sigma$ pure and $l=2$; i.e., it is a PCS. Indeed, by introducing modes with creation operators,

[TABLE]

the Leggett state takes on the standard PCS form, but with each pair of modes for $k\geq 1$ sharing a Schmidt coefficient $\lambda_{k}$. It is worth mentioning that the ground state of a spin-1 Bose gas with an antiferromagnetic interaction takes a form similar to Eq. (24) Ho and Yip (2000); the fermion version of this state was used by Leggett et al. to treat the BEC-BCS crossover problem Leggett (1980); Leggett and Zhang (2012) and by others to discuss the composite boson problem Combescot et al. (2003); Law (2005); Combescot et al. (2008); Tichy et al. (2014).

Dziarmaga and Sacha Dziarmaga and Sacha (2003) generalized Leggett’s wavefunction to the inhomogeneous case (i.e., translational symmetry is broken) while retaining a similar pair-correlated form,

[TABLE]

where $a^{\dagger}_{0}$ and $a^{\dagger}_{m}$ are the creation operators for the condensate mode and the modes orthogonal to the condensate wavefunction; the orthogonal modes and the corresponding real numbers $\lambda_{m}$ are derived from the Bogoliubov Hamiltonian by using a singular-value decomposition. Later, Dziarmaga and Sacha generalized their results to include time-dependent evolutions Dziarmaga and Sacha (2006), where they showed that the state of the system retains the same structure if it starts in a Bogoliubov vacuum state. Although derived from a different perspective, the wavefunction (26) is essentially identical to the PCS introduced in Eq. (19). We derive analytical expressions for the physical quantities (particularly, the 1RDM and 2RDM) associated with the PCS in the large-$N$ limit when more than one mode is macroscopically occupied, i.e., $\lambda_{m}\sim 1+O(1/N)$ for several modes, so the Bogoliubov approximation fails. These results enable us to analyze a PCS with constant computational cost (independent of $N$), whereas a naïve full numerical simulation of the same quantum state requires a computational resource of order $N^{2\nu-2}$. Moreover, in App. C, we present an iterative algorithm to calculate the RDMs of PCS exactly with computational cost $O(N)$.

We discuss the single-particle, two-particle, and higher-order RDMs of the $N$-particle PCS in Secs. IV–VI and present PCS results for the two-site Bose-Hubbard model in Sec. VII. We plan to present a set of coupled equations to determine the Schmidt coefficients and orbitals for the PCS ground state, as well as to formulate a generalized GP-type equation for the dynamics of the two-particle PCS wave function, in a subsequent paper Jiang and Caves .

III PCS and the Bogoliubov Approximation

When the Schmidt rank of the two-particle state (11) is 1, the PCS (19) is a product of single-particle states and reduces to the product-state Ansatz that is the basis of the mean-field GP equation. In the case that the Schmidt rank is greater than 1, but there is a single dominant Schmidt coefficient, Dziarmaga and Sacha Dziarmaga and Sacha (2003) showed that the PCS reproduces the particle-number-conserving ($N$-conserving) Bogoliubov approximation Girardeau and Arnowitt (1959); Gardiner (1997); Castin and Dum (1998). This result validates using PCS for BECs with small depletion, and it provides a different way to formulate the $N$-conserving Bogoliubov approximation. We give an alternative demonstration of this equivalence in this section.

The Bogoliubov approximation fails when depletion from a single condensate mode is large and, hence, the picture of quasi-particle excitations on top of a condensate becomes invalid. The PCS with more than one dominant Schmidt coefficient (a multimode condensate) allows for correlations beyond the weak correlations present in the Bogoliubov approximation. In Sec. VII, we discuss the superfluid-Mott insulator transition in the two-site Bose-Hubbard model using a large-$N$ expansion of the PCS states. In this expansion about a multiple-mode PCS, we find that terms beyond Bogoliubov order must be retained to get an accurate description of the transition.

Consider a Bose-Einstein condensate of $N$ particles in a trapping potential. We restrict our considerations here to the ground state of the condensate, but similar considerations apply to mean-field dynamics and accompanying Bogoliubov excitations. The Gross-Pitaevskii (GP) ground state takes the form

[TABLE]

where the subscript $\perp$ denotes all the modes that are orthogonal to the condensate mode. The essence of the $N$-conserving Bogoliubov approximation is to perturb about the state (27) by introducing the operators

[TABLE]

where $a_{0}$ is the annihilation operator of the condensate mode and $b_{j}^{\dagger}$ is the creation operator of the $j$th orthogonal mode. Even though the operators $\tilde{b}_{j}^{\dagger}$ and $\tilde{b}_{j}$ are not exactly creation and annihilation operators, they satisfy the canonical commutation relations approximately in the limit of large $N$ and small depletion (small number of particles not in the condensate mode). Indeed, for a condensate with small depletion, $a_{0}$ is of order $\sqrt{N}$ and the $\tilde{b}_{j}$s are of order $1$; the modification of the ground-state energy in the Bogoliubov approximation is also of order $1$. When we do operator algebra at Bogoliubov order, we can regard $a_{0}$ and $a_{0}^{\dagger}$ and the $\tilde{b}_{j}$s and $\tilde{b}_{j}^{\dagger}$s as making up a canonical set of creation and annihilation operators and neglect the small corrections to the canonical commutation relations, all of which are order $1/\sqrt{N}$ or smaller.

The $N$-conserving Bogoliubov Hamiltonian, with the mean field removed, is a quadratic function of $\tilde{b}_{j}^{\dagger}$ and $\tilde{b}_{j}$,

[TABLE]

where $H^{\dagger}=H$ and $(H^{\prime})^{T}=H^{\prime}$. This Hamiltonian can be diagonalized by a Bogoliubov transformation,

[TABLE]

where $\mathcal{B}$ is the Gaussian unitary operator corresponding to a symplectic transformation of the operators $\tilde{b}_{j}^{\dagger}$ and $\tilde{b}_{j}$ for $j=1,2,\ldots,\nu-1$. We can use the Bloch-Messiah reduction theorem Bloch and Messiah (1962) to decompose the Gaussian unitary $\mathcal{B}$ as

[TABLE]

where $\mathcal{V}^{\dagger}$ and $\mathcal{U}$ are multiport beamsplitters (which preserve the number of particles and thus have no effect on the vacuum), and the operators in the product are single-mode squeeze operators,

[TABLE]

The ground state of the Bogoliubov Hamiltonian (29) thus takes the form

[TABLE]

Note that we have replaced $\mathcal{V}^{\dagger}$ with $\mathcal{U}^{\dagger}$ in the final form. Neither $\mathcal{U}^{\dagger}$ nor $\mathcal{V}^{\dagger}$ has any effect on the GP ground state $\lvert\,N\,\rangle_{0}\otimes\lvert\,\mathrm{vac}\,\rangle_{\perp}$, so we can include either or neither of them on the right side of the product of the squeeze operators.

The action of the multiport beamsplitter $\mathcal{U}$ is given by

[TABLE]

where $U$ is the unitary matrix that specifies the multiport beamsplitter. The $a_{j}$s are the annihilation operators for a new set of modes orthogonal to the condensate mode, in terms of which we can write the Bogoliubov ground state (33) as

[TABLE]

Now we apply the “quasi-normal-ordered” factored form of the squeeze operator Perelomov (1977); Hollenhorst (1979); Schumaker (1986),

[TABLE]

to the squeeze operators in Eq. (35), with the result

[TABLE]

Note that we can always make the coefficients $\mathord{-}\tanh\gamma_{k}$ in Eq. (37) positive by redefining the phase of $\tilde{b}_{k}$ to make $\gamma_{k}$ negative.

On the other hand, if we separate out a single, dominant Schmidt coefficient $\lambda_{0}$, the PCS (19) of $N=2n$ particles takes the form

[TABLE]

We have made no approximations as yet, but we now use $N!/(N-M)!\simeq N^{M}$, good when $M\ll N$, which holds in our case when $N=2n$ is large and the depletion $M=2m$ is small, to write

[TABLE]

Plugging this into Eq. (38) yields

[TABLE]

where the final approximation is good when $N$ is large. Equations (37) and (40) can be made the same by choosing

[TABLE]

Hence, as promised, the PCS (19) encompasses the Bogoliubov approximation to the BEC ground state.

Another way to display the same result is to notice that if $\mathcal{P}_{N}$ is the projection operator onto the $N$-particle sector and $\lvert\,\alpha\,\rangle_{0}$ is a coherent state for the condensate mode, with $\alpha=\sqrt{N}$, then $\lvert\,N\,\rangle_{0}\otimes\lvert\,\mathrm{vac}\,\rangle_{\perp}\propto\mathcal{P}_{N}\lvert\,\alpha\,\rangle_{0}\otimes\lvert\,\mathrm{vac}\,\rangle_{\perp}$. Since $\tilde{a}_{k}^{\dagger}$ is number conserving, we can write the approximate PCS of Eq. (40) as

[TABLE]

The state before application of the projection into the $N$-particle sector coincides with the “extended catalytic state” that we have introduced for a treatment of the number-conserving Bogoliubov approximation Jiang and Caves (2016).

Now that we have the squeeze operator and its quasi-normal-ordered form available, we can display some general relations among the states we have introduced. In particular, a tensor product of single-mode squeezed states can be written as

[TABLE]

where $\lambda_{j}=-\tanh\gamma_{j}$. A product of single-mode squeezed vacua generates PCS as the terms in its expansion into number sectors, thus, as noted above, making a product of single-mode squeezed vacua an extended catalytic state for a PCS. The Valatin-Butler state (22) is seen to be a product of single-mode squeezed vacua once it is written in terms of the decoupled modes (25), with the modes except for $k=0$ coming in pairs that share the same amount of squeezing.

IV Normalization Factor

The importance of the normalization factor $\mathsf{N}$ in the manipulation and analysis of PCS is analogous to the utility of the partition function in statistical physics. By taking derivatives of the normalization factor with respect to the Schmidt coefficients $\lambda_{j}$, one can calculate the reduced (marginal) density matrices (RDMs) in the Schmidt basis, and these, given the Schmidt-basis wave functions (orbitals), give all the physical observables.

Generally, for any density operator $\rho$ in the symmetric subspace, the $q$-particle RDM $\rho^{(q)}$ has the following matrix elements:

[TABLE]

We normalize the $q$RDM to $N!/(N-q)!=N(N-1)\cdots(N-q+1)$ unless stated otherwise; normalized in this way, the RDMs are clearly identically equal to normally ordered correlation functions. For a PCS, we default to the Schmidt basis for the matrix elements of the $q$RDM:

[TABLE]

In the second-quantized picture, the normalization factor $\mathsf{N}$ introduced in Eq. (19) takes the form

[TABLE]

where we insert a basis of coherent states $\lvert\,\vec{\alpha}\,\rangle=\lvert\,\alpha_{1}\,\rangle\otimes\cdots\otimes\lvert\,\alpha_{\nu}\,\rangle$ for the Schmidt modes. Expanding the monomial in Eq. (47), we have

[TABLE]

where $u_{j}$($v_{j}$) are nonnegative integers satisfying $\sum_{j=1}^{\nu}u_{j}=\sum_{j=1}^{\nu}v_{j}=n$. Although the sum in Eq. (48) appears to require exponential time to evaluate, it can be evaluated in polynomial time by using an iterative algorithm. It is still, however, computationally demanding for large $N$, in addition to being unintuitive.

It is useful to note here that in the case of two Schmidt modes ($\nu=2$), Eq. (48) reduces to

[TABLE]

Here

[TABLE]

is a Gauss hypergeometric function [often denoted ${}_{2}F_{1}(a,b;c;z)$, but missing the two subscripts here for brevity and clarity], which for the particular parameters in Eq. (49) has a series that terminates and is thus a polynomial. Notice that \mathsf{N}/\Gamma(2n+1)=\lambda_{1}^{2n}F\big{(}\frac{1}{2},-n;\frac{1}{2}-n;\lambda_{2}^{2}/\lambda_{1}^{2}\big{)}=\lambda_{2}^{2n}F\big{(}\frac{1}{2},-n;\frac{1}{2}-n;\lambda_{1}^{2}/\lambda_{2}^{2}\big{)}, reflecting the symmetry under exchange of labels of the two Schmidt modes. We apply the results of this section to this expression in App. B to find the exact 1RDM and 2RDM for Schmidt rank $\nu=2$.

To make progress in interpreting and evaluating $\mathsf{N}$, we now take what might be construed as a backward step by writing $\mathsf{N}$ in a different integral form. To do so, we use $(2u_{j})!/u_{j}!=2^{u_{j}}(2u_{j}-1)!!$ to write

[TABLE]

and we then put this into Eq. (48) to give

[TABLE]

We use the above expression to evaluate the normalization factor in the large-$N$ limit in the next section. For now, however, we employ Eq. (52) to derive a generating function,

[TABLE]

An alternate way to obtain the generating function (53) is to evaluate the quantity $\langle\,\mathrm{vac}\,\rvert e^{\sqrt{\tau}\mathcal{A}/2}e^{\sqrt{\tau}\mathcal{A}^{\dagger}/2}\lvert\,\mathrm{vac}\,\rangle$ using the “quasi-normal-ordered” factored form of the squeeze operator found in Eq. (36). Indeed, we can read off the result directly from Eq. (44),

[TABLE]

which is equivalent to Eq. (53).

The normalization factor can be used to find expressions for the diagonal elements of RDMs. As a first example, we are able to represent, using Wick’s theorem, the $j$th diagonal element of the 1RDM $\rho^{(1)}$ in terms of the normalization factor:

[TABLE]

In the second line of Eq. (55), the term and its complex conjugate, which correspond to contracting $a^{\dagger}_{j}$ and $a_{j}$ with the pair annihilation and creation operators, are equal. In addition, by using Wick’s theorem, it is not hard to prove that all the off-diagonal elements of $\rho^{(1)}$ in the Schmidt basis are zero; therefore, the normalization factor and its first derivative determine the 1RDM.

We can apply the same thinking more generally to the $q$RDM matrix elements (46). Indeed, Wick’s theorem [or just thinking hard about the form of Eq. (46)] shows that $\rho_{k_{1}\cdots k_{q},\,j_{1}\cdots j_{q}}^{(q)}$ is nonzero only when for all $j=1,2,\ldots,\nu$, the number of times index $j$ occurs in the matrix element,

[TABLE]

is even. We use this result in Sec. VI.3 to help find the off-diagonal matrix elements of $\rho^{(q)}$ in terms of the diagonal elements.

Returning to diagonal matrix elements, we have the following result for the diagonal elements of the $q$RDM $\rho^{(q)}$:

[TABLE]

This result is most easily proved by mathematical induction. We already have that it holds for $q=1$, so to show that it holds for all positive integers $q$ is that if it holds for $q$, it is satisfied for $q+1$. The inductive hypothesis is thus that

[TABLE]

By taking derivatives with respect to $\lambda_{j_{q+1}}$ of both sides of Eq. (58), we have

[TABLE]

which is the required result. We note that for $q>1$, the $q$RDM is generally not diagonalized in the Schmidt basis. In Sec. VI.3, we show how to construct the entire $q$RDM using only the diagonal elements calculated by the above method.

Although we can avail ourselves of the power of Wick’s theorem to derive the results of this section, the content of Wick’s theorem here and elsewhere in this paper is contained in two commutators, $[a_{j},(\mathcal{A}^{\dagger})^{n}]=2n\lambda_{j}a_{j}^{\dagger}(\mathcal{A}^{\dagger})^{n-1}$ and $[\mathcal{A}^{n},a_{j}^{\dagger}]=2n\lambda_{j}\mathcal{A}^{n-1}a_{j}$, and the following relations, which follow immediately from the commutators:

[TABLE]

V Limit of Large Particle Number

V.1 Diagonal matrix elements of qRDMs in the large-N limit

Often there are thousands to millions of atoms in a BEC, and it is informative to have available results that are only valid in the large-$N$ limit. In this section, we discuss how to derive an asymptotic form of the normalization factor for large $N$. To get the desired analytical results, terms that are of order $1/N$ smaller than the leading terms are neglected. The depletion predicted by the Bogoliubov approximation is of order $1$; i.e., it modifies the 1RDM by order $1$, which is $1/N$ times smaller than the 1RDM itself. This means that the large-$N$ results derived in this section do not encompass the Bogoliubov approximation; instead we should think of them as a generalization of the mean-field GP description to a multimode condensate.

Before the analytical calculation, let us look first at some numerical examples to get some intuition. In Fig. 1, the eigenvalues of the 1RDM (normalized to unity instead of to $N$) are plotted as a function of the number $N$ of particles for $\lambda_{1}^{2}=0.5$, $\lambda_{2}^{2}=0.3$, and $\lambda_{3}^{2}=0.2$. For $N=2$, the eigenvalues of the 1RDM are, of course, equal to the square of the $\lambda_{j}$s, but as $N$ gets larger, the eigenvalues become farther apart. Eventually, the biggest eigenvalue approaches 1, leaving the other two eigenvalues negligible; thus, as far as the 1RDM can tell, the PCS becomes an uninteresting product state for large $N$.

In the second numerical example, plotted in Fig. 2, we consider the situation where there are two $\lambda_{j}$s that are nearly degenerate, $\lambda_{1}^{2}=0.41$ and $\lambda_{2}^{2}=0.39$, and a third smaller value, $\lambda_{3}^{2}=0.2$. When $N$ is of the order $1/(\lambda_{1}-\lambda_{2})=40$ or larger, the third eigenvalue has died out, and only the two biggest eigenvalues play much of a role in determining the 1RDM. Our numerics suggest that in the large-$N$ limit, only those $\lambda_{j}$s that are within order $1/N$ of $\lambda_{1}$ survive.

Generally, we speculate that one only needs to keep those $\lambda_{j}$s that are within $1/N$ of $\lambda_{1}$, the largest eigenvalue given our ordering convention (8). The other $\lambda_{j}$s can be omitted without affecting the PCS; i.e., the relevant low-order RDMs are not affected. This speculation is supported not only by the above numerical results, but also by the following analytical results for the large-$N$ limit. The dominant $\lambda_{j}$s, being very close to each other, can be rescaled (i.e., they are no longer normalized to 1) and parametrized as

[TABLE]

where the $s_{j}$s are real parameters of order unity. Notice that

[TABLE]

Because of the rescaling, all the $\lambda_{j}$s are very close to $1$, and their differences are of order $1/N$. Putting Eq. (62) into Eq. (52), we manipulate the normalization factor $\mathsf{N}$ through the following sequence of steps:

[TABLE]

where $d\hskip 0.4pt\Omega$ denotes the area element on the unit $(\nu-1)$-dimensional sphere defined by $\lvert\vec{y}\hskip 0.4pt\rvert=1$. The only approximation here is to replace, in the fourth line of Eq. (64), the power function by the exponential function. For each low-degree monomial of $s_{j}$s, the error in its expansion coefficient as a result of this replacement is of order $1/n$; such error only becomes substantial when the degree of the monomial approaches $n$. This is an excellent approximation for our purpose of calculating low-order RDMs, because the high-degree monomials only affect high-order RDMs.

Now denote the Gaussian integral in Eq. (64) by

[TABLE]

Since the area of the unit $(\nu-1)$-dimensional sphere is $2\pi^{\nu/2}/\Gamma(\nu/2)$, we have $\Upsilon(0)=1/\Gamma(\nu/2)$. In terms of $\Upsilon(\vec{s})$, the normalization factor takes the form

[TABLE]

The significance of this expression is that the dependencies on $n$ and on $\vec{s}$ (or $\vec{\lambda}$) factorize. According to Eq. (57), the diagonal elements of the RDMs can now be expressed approximately, with errors of order $1/N$, as

[TABLE]

The function \Upsilon\big{(}\vec{s}\,\big{)} determines the diagonal elements of the $q$RDMs in the large-$N$ limit, with all the $n$ dependence removed to the factor $(2n)^{q}$. Relative to exact expressions like Eq. (48), the complexity of evaluating \Upsilon\big{(}\vec{s}\,\big{)} is dramatically reduced because of the removal of the $n$ dependence.

As an example, the $k$th eigenvalue of the 1RDM in the large-$N$ limit takes the form

[TABLE]

Consider the case where $\lambda_{\nu}$, the smallest eigenvalue given our ordering convention (8), satisfies $1-\lambda_{\nu}\gg 1/2n$ and all the other $\lambda_{j}$s are within roughly $1/2n$ of $1$. In this situation $s_{\nu}$ is a negative number with large magnitude, and all the other $s_{j}$s are of order unity. Because of the exponential function in Eq. (68), the magnitude of $y_{\nu}$ must be very small to contribute to the integral; this tells us that $\rho_{\nu\nu}^{(1)}\ll 2n$. On the other hand, we have $\sum_{j=1}^{\nu-1}y_{j}^{2}\simeq 1$ for the other dimensions, so $\rho_{jj}^{(1)}$ for $j=1,2,\ldots,\nu-1$ can be calculated by neglecting the last dimension $y_{\nu}$; in effect, the integral in Eq. (68) is reduced to an integral over a hypersphere of one less dimension. This argument can be easily generalized to higher-order RDMs, and it confirms our speculation that we need only keep those $\lambda_{j}$s that are within $1/N$ of $\lambda_{1}$.

Although we have already made life easier by introducing \Upsilon\big{(}\vec{s}\,\big{)}, it is still a difficult task to evaluate the Gaussian integral (65) over the hypersphere. Fortunately, we can reduce the expression for \Upsilon\big{(}\vec{s}\,) to a single-variable integral. To do so, notice that

[TABLE]

where we do the substitution $\chi=r^{2}/2$ and where $\tau>s_{1}$ for convergence ($s_{1}$ is the largest of the $s_{j}$s). Because Eq. (69) is the Laplace transform of the function \chi^{\nu/2-1}\,\Upsilon\big{(}\chi\hskip 0.4pt\vec{s}\,\big{)}, we have

[TABLE]

where $\mathcal{L}^{-1}$ stands for the inverse Laplace transformation and the real parameter $\Delta>s_{1}$ for convergence. We have thus succeeded in reducing the high-dimensional integral (65) to the one-dimensional integral (70).

For numerical calculations, one might find a straightforward series expansion of the function $\Upsilon(\vec{s}\,)$ to be useful:

[TABLE]

The integral in Eq. (71) can be manipulated by a change of variables into the form

[TABLE]

where

[TABLE]

is the multivariable Beta function [notice that $(-1)!!=1$]. Putting this back into Eq. (71) gives

[TABLE]

It is worth pointing out that we can also represent the function \Upsilon\big{(}\chi\hskip 0.4pt\vec{s}\,\big{)} as a convolution. This representation is particularly useful in Sec. V.2, when we discuss exactly solvable examples. Note that

[TABLE]

where

[TABLE]

Putting Eq. (75) into Eq. (69), we have

[TABLE]

where $*$ stands for the convolution,

[TABLE]

Doing the inverse Laplace transformation, we have

[TABLE]

We now have four representations of \Upsilon\big{(}\hskip 0.4pt\vec{s}\,\big{)}, expressed in Eqs. (65), (70), (74), and (79). All of these turn out to be useful, and we use whichever is most convenient.

We turn now to an exploration of relations among the RDMs (in the large-$N$ limit) that can be derived and expressed through the function \Upsilon\big{(}\vec{s}\,\big{)}. First, from the definition (65) or from the Laplace transform (70), we have

[TABLE]

where $\vec{1}=(1,1,\ldots,1)^{T}$ and $\delta$ is a $c$-number. The only effect of adding a constant $\delta$ to all the $s_{j}$s, which changes the normalization of the $\lambda_{j}^{2}$s by $\nu\delta/n$, is to change \Upsilon\big{(}\vec{s}\,\big{)} and, hence, $\mathsf{N}_{\vec{s},n}$ by multiplying by a factor $e^{\delta}$. This trivial fact implies that

[TABLE]

which applied to Eq. (67), confirms the normalization condition for the 1RDM:

[TABLE]

Equation (81) can be generalized to

[TABLE]

which corresponds to the partial-trace condition on the higher RDMs:

[TABLE]

Notice that we have the factor $2n$, instead of $2n-q$, in Eq (84); this is a consequence of and illustrates the approximations we have used, which are fine for large $n$ and $q$RDMs such that $q\ll N$.

A summary of these considerations is that $\sum_{j=1}^{\nu}\partial/\partial s_{j}$ is the unit operator on \Upsilon\big{(}\vec{s}\,\big{)} and its derivatives. Thus, letting

[TABLE]

we can say that \Upsilon\big{(}\vec{s}\,\big{)} has the dependence $e^{s_{+}}$ on the sum of the PCS parameters; as a consequence, the RDMs do not depend on $s_{+}$, and the differential operator $\sum_{j=1}^{\nu}\partial/\partial s_{j}$ gives zero when acting on RDMs.

More powerful relations for the RDMs can be derived using the Laplace form (70),

[TABLE]

For $s_{j}\neq s_{k}$, we find

[TABLE]

Equations (67) and (87) together give the following relation between the single and two-particle RDMs,

[TABLE]

which we rederive in its exact form, using Wick’s theorem, in Sec. VI.1 [see Eq. (128)]. We can also write $\rho^{(2)}_{j\hskip 0.4ptk,\,j\hskip 0.4ptk}$ in terms of derivatives of the 1RDM,

[TABLE]

which can be used to evaluate $\rho^{(2)}_{j\hskip 0.4ptk,\,j\hskip 0.4ptk}$ when $s_{j}=s_{k}$ or even when $j=k$. Once we know $\rho^{(2)}_{j\hskip 0.4ptk,\,j\hskip 0.4ptk}$ for all $k\neq j$, an alternative way to find $\rho^{(2)}_{j\hskip 0.4ptj,\,j\hskip 0.4ptj}$ is by using the marginalization condition,

[TABLE]

Our conclusion is that the diagonal elements of the 2RDM can be calculated from the diagonal elements of the 1RDM. We consider the off-diagonal elements of the 2RDM in Sec. VI.1.

In view of the results in this section, we should squarely address the question of whether the 1RDM encodes all the information about a PCS. Clear at the start of this discussion is that a PCS is determined by the single-particle Schmidt states $\lvert\,\psi_{j}\,\rangle$, the Schmidt coefficients $\lambda_{j}$, and the number of particles, $N=2n$. Yet the entire structure of RDMs, when written in the Schmidt basis, has exactly the same form, regardless of the Schmidt basis; the relations among the Schmidt-basis RDMs depend only on $n$ and the Schmidt coefficients. The 1RDM is diagonal in the Schmidt basis, and its eigenvalues encode the Schmidt coefficients, so we can determine all the other RDMs, in the Schmidt basis, in terms of the 1RDM.

Does this mean that the 1RDM determines all the physical properties of a PCS? The answer to this is a resounding no. Given a 1RDM in some arbitrary basis, one can diagonalize it, extract the Schmidt coefficients (knowing $n$), and determine the Schmidt orbitals up to a phase for each orbital (equivalently, the phase of the corresponding creation operator). The phase of one orbital (or the associated creation operator) can be fixed by using the overall phase freedom of the PCS, but the PCS depends on the remaining phases, so the 1RDM is clearly not sufficient to determine the PCS. Nonetheless, the form of the RDMs, in the Schmidt basis, is independent of these phases; the phases show up in the RDMs, or correlation coefficients, when they are written in other bases. This is particularly clear when the $q$RDM is written in the position basis:

[TABLE]

Even when one takes into account that for PCS the off-diagonal elements vanish unless each index $j$ occurs an even number of times, this still allows the possibility that the phases of the orbitals appear in the position-basis $q$RDM.

V.2 Examples

In the following, we give some exactly solvable examples in the large-$N$ limit, which include the case of Schmidt rank $\nu=2$ and the cases where the $s_{j}$ coefficients are either totally or pair-wise degenerate.

V.2.1 PCS with Schmidt rank two

For the $\nu=2$ case, we notice that

[TABLE]

where

[TABLE]

is the zeroth-order modified Bessel function. Putting Eq. (92) into Eq. (79), we have

[TABLE]

where

[TABLE]

Notice that if we add $\delta=-(s_{1}+s_{2})/2=-s_{+}$ to both $s_{1}$ and $s_{2}$, as in Eq. (80), we would remove $s_{+}$ from \Upsilon\big{(}s_{1},s_{2}\big{)}.

It is now straightforward to calculate the 1RDM using Eq. (67),

[TABLE]

Similarly, we have

[TABLE]

These equations can be compared with numerical evaluation of the eigenvalues of the 1RDM. As shown in Fig. 3, the two conform quite well in the large-$N$ limit.

To examine particle-particle correlations in the $\nu=2$ PCS, we calculate the 2RDM. By putting Eq. (94) into Eq. (67), we have

[TABLE]

similarly, we have

[TABLE]

It is straightforward to check that these 2RDM elements marginalize correctly, i.e., $\rho_{11,11}^{(2)}+\rho_{12,12}^{(2)}=2n\,\rho_{11}^{(1)}$ and $\rho_{12,12}^{(2)}+\rho_{22,22}^{(2)}=2n\,\rho_{22}^{(1)}$, where we neglect the difference between $2n$ and $2n-1$ in the large-$N$ limit. The relation (88) can be verified by using the recursion relations of the Bessel functions,

[TABLE]

V.2.2 Totally degenerate coefficients

In the totally degenerate case ($s_{1}=s_{2}=\cdots=s_{\nu}$), all of the eigenvalues of the 1RDM are the same and thus are determined by the normalization condition,

[TABLE]

Similarly, we have that the 2RDM matrix elements $\rho_{jj,\,jj}^{(2)}$ are the same for all $j$ and the matrix elements $\rho_{j\hskip 0.4ptk,\,j\hskip 0.4ptk}^{(2)}$ are the same for all $j\neq k$; moreover, putting the Laplace form (86) into Eq. (67) implies that

[TABLE]

We can now use the marginalization condition (84) to determine that

[TABLE]

In this situation, plugging Eq. (74) into Eq. (67) gives us the diagonal matrix elements of the RDMs of all orders,

[TABLE]

where $q_{j}=\sum_{k=1}^{q}\delta_{j,j_{k}}$ is the number of times $j$ appears in the list of single-particle states, $j_{1},\ldots,j_{q}$.

The degenerate case illustrates that when the system possesses special symmetries, the expressions for matrix elements can sometimes be solved exactly.

V.2.3 Pair-degenerate coefficients

Another example occurs when the $s_{j}$s come in degenerate pairs, i.e., $s_{j}=s_{j+1}$ for odd $j$. Note that Eq. (92) gives the following for $s_{1}=s_{2}$:

[TABLE]

Using Eq. (79) and (106), we have

[TABLE]

where $\mathrm{odd}=\{1,3,\ldots,\nu-1\}$. The convolutions in Eq. (107) correspond to the sum of many exponential distributions (hypoexponential distribution); results from probability theory can then be used to derive the final form. Another way of deriving Eq. (107) is to use the residue theorem.

VI Higher-Order Reduced Density Matrices

VI.1 2RDMs

In preceding sections, we focused on the diagonal elements of RDMs. In this section we first show how to represent the off-diagonal elements of the 2RDM in terms of its diagonal elements and thus to relate the entire 2RDM to the 1RDM. We do this exactly, without making the large-$N$ approximation in this section, except where we explicitly introduce that approximation to illustrate how the exact results simplify when $N$ is large. We then go on to showing how to find the off-diagonal elements of $q$RDMs from the diagonal elements in the large-$N$ limit.

In the Schmidt basis, the 2RDM reads

[TABLE]

All the matrix elements of $\rho^{(2)}$ are real, and $\rho^{(2)}$ satisfies the normalization condition

[TABLE]

Using Wick’s theorem (or doing some hard thinking about the definition of the PCS), one finds that the 2RDM has the form

[TABLE]

The matrices $\xi$ and $\xi^{\prime}$ are real because none of the creation and annihilation operators introduces any complex numbers into Eq. (108); they are symmetric as a consequence of the Hermiticity of the 2RDM.

There are two instructive ways to write the 2RDM, which avoid the blizzard of Kronecker deltas in Eq. (110). The first is to write out $\rho^{(2)}$ as an operator in the Schmidt basis:

[TABLE]

A second instructive way to state the content of Eq. (110) is that the only nonzero matrix elements of the 2RDM are the diagonal elements,

[TABLE]

and the following specific off-diagonal elements,

[TABLE]

These latter are only off-diagonal when $j\neq k$; all the formulas coincide for $\rho_{jj,jj}^{(2)}$. Note that the off-diagonal elements in Eq. (114) are equal to the diagonal elements (112) because of the symmetry under the exchange of the two particles. The relations (114) can be inverted in the following way:

[TABLE]

The diagonal elements of $\xi$ and $\xi^{\prime}$ cannot be determined separately, but we can always choose, for example, $\xi_{jj}=\xi^{\prime}_{jj}$.

To proceed, we derive an exact result for how the 2RDM changes under exchange of $j_{1}$ and $k_{1}$,

[TABLE]

If we make explicit that $\rho^{(1)}$ is diagonal, this exchange takes the form

[TABLE]

This expression contains information about $O(n)$ contributions to the 2RDM, whose size is generally $O(n^{2})$, whereas the 1RDM is $O(n)$. In the large-$N$ limit, where $\lambda_{j}/\lambda_{k}=1+O(1/n)$ [see Eq. (62)], the expression (117) simply says that $\rho_{k_{1}k_{2},\,j_{1}j_{2}}^{(2)}$ is symmetric at leading order under the exchange of $j_{1}$ and $k_{1}$. The useful content of Eq. (117) comes from specializing to $j_{1}=j_{2}=j$ and $k_{1}=k_{2}=k$:

[TABLE]

This expression relates the two kinds of off-diagonal elements of the 2RDM, and since $\rho_{jk,kj}^{(2)}=\rho_{jk,jk}^{(2)}$, it also allows us to relate all the off-diagonal elements of the 2RDM to the 1RDM and diagonal elements of the 2RDM:

[TABLE]

To see what is going on more clearly, it is useful to manipulate Eq. (118) into an equivalent form in which the Schmidt coefficients no longer multiply the 2RDM matrix elements. To do this, we use the condition $\rho^{(2)}_{kk,jj}=\rho^{(2)}_{jj,kk}$ to write

[TABLE]

Combining Eqs. (119) and (120), we have

[TABLE]

In the large-$N$ limit, we have, from Eq. (62), (\lambda_{j}-\lambda_{k})\big{/}(\lambda_{j}+\lambda_{k})=(s_{j}-s_{k})/4n+O(1/n^{2}) and thus

[TABLE]

The 2RDM matrix elements on the left are $O(n^{2})$, whereas the two terms on the right-hand side are $O(n^{1})$ and $O(n^{0})$.

We see that to leading order, $O(n^{2})$, in the large-$N$ limit, we have $\rho^{(2)}_{kk,jj}=\rho^{(2)}_{jk,jk}$, confirming our assertion above that to leading order, $\rho_{k_{1}k_{2},\,j_{1}j_{2}}^{(2)}$ is symmetric under exchange of $j_{1}$ and $k_{1}$. Indeed, using Eqs. (113) and (114), we have that $\xi_{jk}-\xi_{jk}^{\prime}=\rho^{(2)}_{kk,jj}-\rho^{(2)}_{jk,jk}$, so to leading order in $n$, we can conclude that $\xi$ and $\xi^{\prime}$ are the same, and a single matrix $\xi$ determines the 2RDM. At leading order, Eq. (110) simplifies to

[TABLE]

Notice that from Eq. (112), to this level of approximation, $\xi$—and thus the entire 2RDM—is determined by the diagonal elements of $\rho^{(2)}$, and from Eqs. (88) and (89), the diagonal components of the 2RDM can be related to the 1RDM,

[TABLE]

Thus, to leading order in the large-$N$ limit, we have determined the 2RDM in terms of the 1RDM. Notice that Eq. (124) implies that to leading order, all the matrix elements of $\xi$ are nonnegative. Since from Eq. (113) we also have $\rho^{(2)}_{kk,jj}=\xi_{jk}$ for $j\neq k$, we can conclude that the off-diagonal elements $\rho^{(2)}_{kk,jj}$ are nonnegative to leading order. We can make this last conclusion secure to all orders by considering the exact relation (121), which shows that the off-diagonal elements $\rho^{(2)}_{kk,jj}$ are nonnegative to all orders.

We must use the leading-order, large-$N$ approximation with great care, however, because the small terms that are neglected in this approximation are crucial to determining the PCS ground state precisely in those circumstances where the PCS ground state is different from ground state obtained from the product-state (GP) Ansatz. The first term on the right-hand side of Eq. (122)—this is the first correction to the leading-order behavior—is already a Bogoliubov-order correction that is not captured by the large-$N$ limit. The second term on the right-hand side of Eq. (122) is thus a post-Bogoliubov correction. We realize this is not a standard Bogoliubov expansion, because we are expanding about our large-$N$ PCS approximation, not about a single-mode condensate. Nonetheless, we believe that referring to terms in this expansion as Bogoliubov order and post-Bogoliubov order is both useful and informative, so we adopt this terminology in what follows. As we show in Sec. VII, this second term, of order $n^{0}$, although two orders smaller than the leading term in the 2RDM, is crucial in capturing the superfluid-insulator transition in the Bose-Hubbard model.

To complete our program of finding the exact relation of the 2RDM to the 1RDM, we need to find the exact analog of the large-$N$ relation (88) between the diagonal elements of $\rho^{(2)}$ and those of $\rho^{(1)}$. Thus we rederive this relation here using Wick’s theorem, without making the large-$N$ approximation. By doing contractions of the single annihilation and creation operators with the pair creation and annihilation operators, we have for $j\neq k$,

[TABLE]

Switching the roles of $j$ and $k$, we have

[TABLE]

Multiplying Eqs. (126) and (125) by $\lambda_{j}^{2}$ and $\lambda_{k}^{2}$, respectively, and subtracting the results gives

[TABLE]

which leads to

[TABLE]

In the large-$N$ limit, the first term on the right-hand side of Eq. (128) dominates the second term and reproduces Eq. (88). If $\lambda_{j}=\lambda_{k}$, we can fall back on the analog of Eq. (89):

[TABLE]

These relations allow us to determine $\rho_{j\hskip 0.4ptk,\,j\hskip 0.4ptk}^{(2)}$, for $j\neq k$, from the 1RDM, and then $\rho_{j\hskip 0.4ptj,\,j\hskip 0.4ptj}^{(2)}$ is determined by the marginalization condition,

[TABLE]

Equation (121) can then be used to calculate the off-diagonal elements $\rho^{(2)}_{kk,jj}$:

[TABLE]

Equations (128)–(131) provide a recipe for calculating the 2RDM from the 1RDM.

We close this subsection by discussing a subtle point regarding formulas like Eqs. (122) and (128), which involve rewriting exact results that depend on the Schmidt coefficients $\lambda_{j}=\sqrt{1+s_{j}/n}$ in terms of the PCS parameters $s_{j}$. Rescaling the Schmidt coefficients by a factor $r$ rescales the normalization factor $\mathsf{N}$ by a factor $r^{2n}$; this factor disappears from all the RDMs, however, and this property means that we can choose any normalization for the Schmidt coefficients. We also know that in the large-$N$ limit, adding a constant $\delta$ to all the PCS parameters $s_{j}$, rescales $\Upsilon(\vec{s})$ by a factor $e^{\delta}$ and thus does not affect the large-$N$ RDMs. What we have to be careful about in exact expressions is that adding $\delta$ to all the PCS parameters is not a rescaling of the Schmidt coefficients. Instead, we rescale the Schmidt coefficients by a factor $r=\sqrt{1+\delta/n}$ by defining new PCS parameters $s^{\prime}_{j}=(1+\delta/n)s_{j}+\delta$ and thus new Schmidt coefficients $\lambda^{\prime}_{j}=\lambda_{j}\sqrt{1+\delta/n}$. In the large-$N$ limit, the term $\delta/n$ in the primed PCS parameters should be ignored, as it is of the size of terms neglected in the large-$N$ limit; this takes us back to simply adding a constant to the PCS parameter. Exact expressions for RDMs are invariant under the rescaling by $r=\sqrt{1+\delta/n}$, but because the rescaling depends on $n$, when one expands in powers of $1/n$ by introducing the PCS parameters, the terms in the expansion are generally not invariant under the rescaling. In Eq. (122), this does not cause a problem, because the rescaling only introduces terms of higher order than those kept; thus we can derive the Bogoliubov and post-Bogoliubov corrections in Eq. (122). In contrast, in Eq. (128), the rescaling mixes the leading-order term with the Bogoliubov correction; this shows up as a difficulty in identifying the Bogoliubov correction to the large-$N$ contribution.

VI.2 2RDM examples in the large-N limit

In the case $\nu=2$, we can solve for the 2RDM exactly in the large-$N$ limit. Plugging the results (98), (99), and (100) into Eq. (123), we have the following expression for the 2RDM in the Schmidt basis, $\{\,\lvert\,11\,\rangle\,$, $\,\lvert\,12\,\rangle\,$,$\,\lvert\,21\,\rangle\,$,$\,\lvert\,22\,\rangle\,\}$,

[TABLE]

where $I_{0}$, $I_{1}$, and $I_{2}$ are the zeroth, first, and second order modified Bessel functions with argument $s_{-}=(s_{1}-s_{2})/2$. Equivalently, we can write $\rho^{(2)}$ in the Pauli basis,

[TABLE]

The two-particle state (133) is not entangled; i.e., it has zero concurrence. It is known that all pairwise entanglement vanishes in large bosonic systems due to the monogamy of entanglement Koashi et al. (2000).

Another case that can be solved analytically in the large-$N$ limit is the totally degenerate case, i.e., $s_{1}=s_{2}=\cdots=s_{\nu}$. Using Eqs. (104) and (123), we have

[TABLE]

VI.3 Off-diagonal elements of qRDMs in the large-N limit

We already know how to calculate the diagonal elements of the $q$RDMs in the large-$N$ limit, by using Eq. (67), but we have not yet discussed how to derive the off-diagonal elements for $q>2$. It turns out that the methods we developed for the 2RDM are also useful for the $q$RDM, and in the large-$N$ limit, the result is very simple.

First, we establish how the $q$RDM matrix element (46) changes under exchange of $j_{1}$ and $k_{1}$ to leading order, i.e., to $O(n^{q})$:

[TABLE]

In the large-$N$ limit, $\lambda_{j_{1}}/\lambda_{k_{1}}=1+O(1/n)$, so at leading order, we have

[TABLE]

Because of the particle-exchange symmetry, $\rho_{k_{1}\cdots k_{q},\,j_{1}\cdots j_{q}}^{(q)}$ is invariant under any permutation of $\{j_{1},\ldots,j_{q}\}$ and any permutation of $\{k_{1},\ldots,k_{q}\}$. Thus Eq. (136) implies that $\rho_{k_{1}\cdots k_{q},\,j_{1}\cdots j_{q}}^{(q)}$ remains the same at leading order under any permutation of all the indices, $\{j_{1},\ldots,j_{q},\,k_{1},\ldots,k_{q}\}$.

Start now with any nonzero matrix element $\rho_{k_{1}k_{2}\cdots k_{q},\,j_{1}j_{2}\cdots j_{q}}^{(q)}$. Recalling from Eq. (56) that any particular index occurs an even number of times, those occurrences can be permuted to positions where that index is paired up on the row and column sides of the matrix element, thus giving a diagonal matrix element. Thus the original nonzero matrix element is equal at leading order to a diagonal matrix element, which can be obtained from Eq. (67). This relates any off-diagonal element of $\rho^{(q)}$ to a corresponding diagonal element at leading order, i.e., $O(n^{q})$. Exact relations, not assuming that $N$ is large, but of increasing complexity as $q$ increases, could be worked out using the techniques developed in this section.

VII Two-Site Bose-Hubbard Model

To illustrate the utility of the PCS formalism, we investigate the performance of the PCS Ansatz in approximating ground-state properties and RDMs of the two-site Bose-Hubbard (BH) model, given by the Hamiltonian

[TABLE]

which describes the hopping of $N$ interacting bosons between two identical lattice sites. Here $b_{k}(b_{k}^{\dagger})$ represents the annihilation (creation) operator for a boson on site $k=1,2$, $J$ sets the strength of tunneling between sites, and $U>0$ determines the strength of the on-site repulsive interaction between two particles. Despite its simplicity, the two-site BH Hamiltonian (137) describes the rich physics of a quantum phase transition from a superfluid to a Mott insulator Fisher et al. (1989) with increasing ratio $U/(nJ)$, for integer-filling $n$. In this section we focus on situations where $n\gg 1$. For very weak interactions, $(nU)/J\ll 1$, the ground state of the BH model is a superfluid with almost all particles occupying a single (condensate) mode that is symmetrically delocalized between the two sites; this regime is well described by GP mean-field theory. Depletion of the single condensate mode increases as the interaction strength increases, with the transition from a superfluid to a Mott insulator occurring when $U/(nJ)\sim 1$. For strong interactions, $U/(nJ)\gg 1$, the condensate fragments into two uncorrelated components at the two sites. Notice that in terms of powers of $1/n$, there is another regime of relatively weak interaction strength, $U/J\sim 1$, between the transition to mean-field GP at $(nU)/J\sim 1$ and the transition to a Mott insulator at $U/(nJ)\sim 1$. In this intermediate regime, the ground state is still close to a single-mode condensate; we discuss how the PCS formalism treats this intermediate regime in the following. Even though the PCS Ansatz and the number-conserving Bogoliubov approximation share the same form (26), the latter only works for small corrections about a single-mode condensate and fails when there are two or more well-populated modes. The PCS approach is nonperturbative and does not have such a constraint.

Because the size of the Hilbert space grows linearly with $N$, we can directly diagonalize Eq. (137) numerically for large system sizes ($N\gtrsim 100$) and, hence, benchmark the PCS Ansatz against the exact BH ground state for a wide range of parameters. Specifically, we focus on the comparison of the ground-state energy and on the reconstruction of the one- and two-particle RDMs for the ground state. As appropriate, we include comparisons with the ground state of the conventional Bogoliubov approximation to the two-site BH Hamiltonian Paraoanu et al. (2001); the predictions of the conventional Bogoliubov approximation are developed in App. A. Lastly, we briefly explore the potential of the PCS Ansatz by showing that it can faithfully represent the time evolution of the two-site BH model Jiang (2014), which will be discussed in detail elsewhere.

We start here with the exact $N$-particle ground state $\lvert\,\varPsi_{\mathrm{ex}}\,\rangle$ and ground-state energy $E$ for parameters ranging across the transition region from very weak interactions, $(nU)/J\ll 1$, to strong interactions, $U/(nJ)\gg 1$. We then calculate the exact one- and two-particle RDMs,

[TABLE]

Note that the RDMs (138) and (139) are written in the physical basis that defines sites 1 and 2 and, in turn, determine directly important physical properties, such as the on-site number fluctuation,

[TABLE]

the (normalized) second-order site correlation functions,

[TABLE]

and the ground-state energy,

[TABLE]

Before going further, it is useful to spell out the symmetries of the ground-state 2RDM $\sigma^{(2)}$. The 2RDM is, of course, Hermitian. Time-reversal invariance of the Hamiltonian (137) implies that $\sigma^{(2)}$ has real matrix elements; together with Hermiticity, this means that $\sigma^{(2)}$ is a real, symmetric matrix. The boson symmetry implies that the two row labels or the two column labels can be interchanged without changing the value of the matrix elements. These symmetries imply that the 2RDM can be parameterized as

[TABLE]

where all the parameters are real and $\delta=(\delta_{1}+\delta_{2})/2$. We use the basis ordering $\{\lvert\,11\,\rangle,\lvert\,12\,\rangle,\lvert\,21\,\rangle,\lvert\,22\,\rangle\}$. Normalization implies that

[TABLE]

Notice that the boson symmetry means that the matrix (143) is symmetric under interchange of the middle two rows or the middle two columns. There is one further symmetry we can apply. The Bose-Hubbard Hamiltonian (137) is symmetric under interchange of the site labels. The ground state inherits this symmetry, which implies that the matrix $\sigma^{(2)}$ is invariant under an inversion through its center. The one further identification this gives us is that $\sigma^{(2)}_{11,\,11}=\sigma^{(2)}_{22,\,22}$, i.e., $\delta_{1}=\delta_{2}$, which leaves the 2RDM in the form

[TABLE]

The 1RDM in the site basis takes the form

[TABLE]

and the second-order correlation functions (141) and the ground-state energy (142) become

[TABLE]

The connection to the PCS Ansatz is made in the Schmidt basis, in which the 1RDM has the diagonal form

[TABLE]

This diagonal form is obtained via the transformation

[TABLE]

Here $e^{i\mu}$ is the one phase (after removing the global phase freedom) that is left arbitrary in the two Schmidt orbitals. We show below that this phase is fixed (to within an irrelevant sign) by choosing it to minimize the PCS ground-state energy, and this leads to the choice $\mu=\pi/2$, i.e., $e^{i\mu}=i$.

Since $\operatorname{tr}\rho^{(1)}=N$, Eq. (150) is uniquely defined by the population imbalance of the Schmidt modes,

[TABLE]

which we plot in Fig. 4(a) for the exact ground state as a function of $U/(nJ)$. It is important to note that $\Delta$ contains the same information as the off-site correlation function $\sigma^{(1)}_{12}$. For weak interactions, $U/(nJ)\ll 1$, the ground state of the BH model is a superfluid with almost all particles occupying the (condensate) Schmidt mode 1 ($\Delta\sim\varrho_{1}\sim N$), which is symmetrically delocalized between the two sites. As $U/(nJ)$ increases, the dominant Schmidt mode becomes depleted, and the transition to a Mott insulator occurs around $U/(nJ)\sim 1$. For strong interactions, $U/(nJ)\gg 1$, the two Schmidt modes are nearly equally occupied, i.e., $\Delta\sim 0$, and the superfluid fragments into two uncorrelated components at the two sites. In the limit $U/(nJ)\rightarrow\infty$, the ground state becomes the Mott-insulator state $\lvert\,n,n\,\rangle$, i.e., a product state of $n$ particles at each site, as shown by the near-unity fidelity $|\langle\,n,n\,|\,\varPsi_{\rm ex}\,\rangle|^{2}$ in Fig. 4(a) (right axis).

Using the transformation (151), we can write the Hamiltonian (137) as

[TABLE]

In this basis the expectation value of the ground-state energy takes the form

[TABLE]

where \rho^{(2)}_{k_{1}k_{2},\,j_{1}j_{2}}=\big{\langle}\,\varPsi_{\mathrm{ex}}\,\big{\rvert}\,a_{j_{1}}^{\dagger}a_{j_{2}}^{\dagger}a_{k_{2}}a_{k_{1}}\,\big{\lvert}\,\varPsi_{\mathrm{ex}}\,\big{\rangle}. We remind the reader that the subscripts on $\rho$ index Schmidt orbitals, as in all previous sections of the paper, and not sites, as in the subscripts on $\sigma$. Now we use the PCS Ansatz to determine a particular Schmidt basis. By taking the expectation value of Eq. (153) over the PCS Ansatz, one obtains an expression of the same form as Eq. (154), but with the values of the matrix elements given by the PCS formalism. Because the Schmidt-basis matrix elements $\rho^{(2)}_{22,\,11}$ and $\rho^{(2)}_{11,\,22}$ in the PCS formalism are always nonnegative [see discussion following Eq. (124)], the PCS ground-state energy is minimized, as promised, by the phase choice $e^{i\mu}=i$. This specifies the Schmidt basis for the PCS formalism. In the Schmidt basis, the ground-state energy then has the form

[TABLE]

The transition from superfluid to Mott insulator is described by the competition between the first and second terms on the right-hand side of Eq. (155). In the Schmidt basis, the two-site Bose-Hubbard Hamiltonian in the Bogoliubov approximation now becomes

[TABLE]

We can transform the 2RDM (145) to the Schmidt basis using the transformation (151) with $e^{i\mu}=i$, obtaining

[TABLE]

The exact ground-state 2RDM (157) has the same general form (110) as a PCS state. The only use of the PCS Ansatz in deriving the ground-state energy $(\ref{eq:pcs_energy_initial})$ and the 2RDM $(\ref{eq:exact_2rdm_Sch})$ is to determine a particular Schmidt basis. Once that determination is made, the ground-state energy and the 2RDM have the given forms, both for the exact results and for the PCS Ansatz; what changes between the two are the values of the matrix elements in the expressions. The four parameters in the 2RDM are given by Eqs. (144) and (152) and by

[TABLE]

For the present case of Schmidt rank $\nu=2$, we have calculated the PCS normalization factor exactly in Eq. (49), and we use that expression to derive exact PCS RDMs in App. B. In App. C, we show how to evaluate the PCS 1RDM for $\nu\geq 2$ in time $O(n)$ based on a simple iterative method, but we do not use that method here. We center our discussion in this section around the use of the large-$N$ results derived in the previous sections to determine the RDMs approximately. We relegate the calculation of the exact PCS to App. B, but include its predictions in our analysis in order to investigate how well the PCS Ansatz does before any approximations are made to it. To avoid confusion when presenting both predictions, we refer for the remainder of this section—and in Fig. 4 and App. B—to the approximate PCS Ansatz as APCS and reserve PCS for results of the PCS Ansatz without approximations.

We expect the exact PCS to do well for weak interactions, $U/(nJ)\ll 1$, where we get a single condensate by choosing Schmidt coefficients $\lambda_{1}=1$ and $\lambda_{2}=0$; the resulting pair creation operator, $\mathcal{A}^{\dagger}=(a_{1}^{\dagger})^{2}=(b_{1}^{\dagger}+b_{2}^{\dagger})^{2}/2$, creates pairs of particles in the mode $a_{1}$. We also expect the exact PCS to do well for strong interactions, $U/(nJ)\gg 1$, where we get $n$ particles at each site by choosing $\lambda_{1}^{2}=\lambda_{2}^{2}=1/2$; the resulting pair creation operator, $\mathcal{A}^{\dagger}=[(a_{1}^{\dagger})^{2}+(a_{2}^{\dagger})^{2}]/\sqrt{2}=\sqrt{2}b_{1}^{\dagger}b_{2}^{\dagger}$, creates particles pairwise at the two sites. The purpose of the numerics for the exact PCS is thus to test how well the PCS Ansatz does in the transition region, $U/(nJ)\sim 1$, and what we find below is that it does a credible job throughout all of $U/(nJ)$ parameter space, even for relatively small particle numbers.

The situation is different for the approximate PCS. Recall from our previous discussions of two Schmidt modes that all physical quantities in the APCS are functions of the parameter $s_{-}=(s_{1}-s_{2})/2$ and that the approximations made in the large-$N$ limit require that $s_{-}\ll n$. The purpose of the APCS predictions is really to see how well the APCS approximates the exact PCS. In the transition region between a superfluid and a Mott insulator, $U/(nJ)\sim 1$, we expect $s_{-}\sim 1$, and in the strong-interaction Mott regime, $U/(nJ)\gg 1$, we expect $s_{-}\ll 1$, so throughout these two regimes, we expect the APCS to match the exact PCS well even for fairly small particle numbers. In the mean-field GP regime of very weak interactions, $(nU)/J\gtrsim 1$, the APCS attempts to model a single condensate by having $s_{-}$ formally become of order $N$ or larger; this, however, contradicts the definition of PCS parameters in Eq. (62). We cannot expect the APCS predictions to do well or even to be physical when $s_{-}$ gets as big as or bigger than $n$. The reason for this problem, of course, is that the approximations made in the APCS are inconsistent with giving an accurate description of a single condensate mode. What we find below, however, is that there is an intermediate regime of weak interactions, $U/J\sim 1$, where $s_{-}\sim\sqrt{N}$, and in this regime, the APCS does a credible job of modeling a single-mode condensate, with corrections of order $1/\sqrt{N}$. The numerics are roughly consistent with all of this, except that the APCS does worse than expected in the strong-interaction Mott regime, for reasons that we discuss further below.

Comparison of the 2RDM (157) with the large-$N$ approximation (132) gives the following APCS values:

[TABLE]

Here $s_{-}=(s_{1}-s_{2})/2$ is the APCS parameter, and $\chi_{j}(s_{-})=I_{j}(s_{-})/I_{0}(s_{-})$, with $I_{j}$ denoting the $j$th modified Bessel function.

To find the APCS predictions for the ground-state energy, however, we need to go beyond the parameters (159). A glance at the 2RDM (157) and at the ground-state energy (155) shows that if $\delta=0$, then there is no competition between tunneling and the on-site interaction, so we do not capture the transition from a single-mode condensate to a Mott insulator. To model this transition, we hybridize the large-$N$ results with the exact PCS relations derived in Sec. VI.1. We start by using Eq. (122) to get a more accurate value of $\delta$,

[TABLE]

This then reduces the energy (155) to

[TABLE]

We reiterate here some terminology that we employed in previous sections. The large-$N$ approximation that we developed in Sec. V.1 is a generalization of the mean-field description to multimode condensates. Corrections to the large-$N$ approximation, such as those derived for exact 2RDMs in Sec. VI.1, occur as increasing powers of $1/N$. Thus, as we pointed out in Sec. VI.1, we refer to the first term, $-N$, on the right-hand side of Eq. (160) as a Bogoliubov-order correction to the large-$N$ 2RDM and the second term as a post-Bogoliubov correction. It is the post-Bogoliubov correction that gives rise to the quantity $Us_{-}/(4N)$ in Eq. (161). In the transition region, where the on-site interaction competes with the tunneling splitting, the competition is expressed by the two quantities $Us_{-}/(4N)$ and $J$, which are of the same order in the transition region. This stands in marked contrast with the mean-field GP regime of a single condensate mode, which is nominally defined by $(nU)/J\ll 1$ [$U/(nJ)\ll 1/n^{2}$]; i.e., the magnitude of the tunneling splitting, $nJ$, far exceeds the on-site interaction energy, $Un^{2}$. Because the Bogoliubov correction in Eq. (160) is a constant, however, the single-condensate behavior of the mean-field regime extends, at least approximately, to values of $J$ one order smaller in $1/N$, i.e., to the intermediate regime of weak interaction strength, $U/J\sim 1$, as is evident from the plots in Fig. 4(a).

Using the relations (96) and (97), we have

[TABLE]

The APCS energy can thus be written as a function of the parameter $s_{-}$ alone,

[TABLE]

The APCS ground-state energy for the two-site Bose-Hubbard model is found by minimizing this expression with respect to $s_{-}$; the optimal value of $s_{-}$ is determined by

[TABLE]

The optimal value of $s_{-}$ (solved numerically) is then inserted into the parameters (159) to give the APCS 1RDM and 2RDM and into Eqs. (162) and (163) to find the APCS population difference and ground-state energy. We stress that Eq. (160) has a precision of $O(1/N)$ and Eq. (163) has a precision of $O(1)$, both two orders beyond the precision of the large-$N$ approximation discussed in Sec. V.1 and one order beyond the Bogoliubov corrections. The ability of the PCS Ansatz to include selectively such crucial post-Bogoliubov corrections is the primary strength of the formalism, but does indicate that using the formalism requires care if one wants to find physically meaningful results.

We can use the large-argument asymptotic expansions of the modified Bessel functions to find that for $s_{-}\gg 1$, the ground-state value determined by Eq. (164) is

[TABLE]

In the single-condensate, mean-field regime, where $(Un)/J\lesssim 1$, this gives $s_{-}\gtrsim n$, which lies outside the region of validity of the large-$N$ approximation, reflecting the expected fact that the large-$N$ approximation cannot accommodate a single-mode condensate. In the transition region, where $U/(nJ)\sim 1$, this gives $s_{-}\sim 1$, as expected, but meaning that the asymptotic expansions used to find Eq. (165) are not valid. The sweet spot for Eq. (165) is in the intermediate regime of weak interaction strength, where $U/J\sim 1$, which gives $s_{-}\sim\sqrt{N}$; the population difference and ground-state energy in this intermediate regime are

[TABLE]

both of which are close to the values for a single-mode condensate. Notice that the limits of validity of the approximate APCS predictions (166) and (167), i.e., $1\ll s_{-}=\sqrt{2NJ/U}\ll N$, are equivalent to being in the intermediate regime of weak interaction strength, $1/N\ll U/J\ll N$.

We can also investigate the behavior of the APCS in the Mott regime by considering the power-series expansions of the modified Bessel functions for small $s_{-}$. What one obtains are corrections to the Mott population difference and the Mott ground-state energy that are nominally as small as or smaller than the terms neglected in Eqs. (162) and (163). This suggests that though the APCS gives an essentially correct description of the Mott regime for large $N$, it is missing details that are contained in the exact ground state and in the exact PCS ground state; the plots in Fig. 4(c) are consistent with this suggestion. Moreover, the difficulty in finding the APCS ground-state zero of Eq. (163) is a hint that it might be difficult to find the exact PCS ground-state energy in the Mott-insulator regime from Eq. (190). Indeed, we find that Eq. (190) has a very shallow minimum in the Mott regime; this accounts for the noisy data for the exact PCS in the Mott-insulator regime in Fig. 4(c).

In Fig. 4(b), we plot the APCS population difference $\Delta_{\mathrm{apcs}}$ from Eq. (162), the exact PCS population difference $\Delta_{\mathrm{pcs}}$ from Eq. (185), and also the exact solution $\Delta$ as functions of $U/(nJ)$ for $N=50,100,500,1\,000$. The corresponding relative error is shown in Fig. 4(c), which gives a first indication of the accuracy of the PCS formalism. The agreement of the PCS approximation is remarkably good over the whole parameter space and for all particle numbers considered. In fact, the APCS gives good predictions across the entire parameter space, improving with increasing $N$, but fairly good even for $N=50$. In the transition region, the APCS matches the exact PCS for $N\gtrsim 500$, but is not as good as the exact PCS, for any value of $N$, in the regimes of weak or strong interaction strength $U$.

In Fig. 5, we compare the APCS ground-state energy (163) and the PCS ground-state energy (190) to the exact ground-state energy, by computing their relative error as a function of $U/(nJ)$, for $N=50,100,500,1\,000$. The error of the PCS approximation is impressively small for all parameters considered. Note that the largest discrepancy (of only a few percent) of the APCS approximation happens for $U/(nJ)\sim 0$ and $N=50$, a regime where we do not expect the APCS to do well, because it cannot describe a single-mode condensate, and where, in any case, the ground state, as a single-mode condensate, is trivial and the PCS Ansatz is unnecessary. Just as for the population imbalance, the APCS matches the exact PCS in the transition region for $N\gtrsim 500$, but is not as good as the exact PCS, for any value of $N$, in the regimes of weak or strong interaction strength $U$.

We turn our attention now to the PCS construction of the two-particle RDM for the ground state. In the Schmidt basis (151), the exact 2RDM $\rho^{(2)}$ and the PCS 2RDM, in its exact form $\rho_{\mathrm{pcs}}^{(2)}$, given in App. B, and its approximate, large-$N$ form $\rho_{\mathrm{apcs}}^{(2)}$, constructed from the parameters (159), have the general form (157). Thus we can assess the accuracy of the PCS 2RDMs by comparing them directly to the exact 2RDM. First, to serve as overall error measures, we compute the trace distance and the infidelity between the PCS RDMs, $\rho_{\mathrm{pcs}}^{(2)}$ and $\rho_{\mathrm{apcs}}^{(2)}$, and the exact 2RDM $\rho^{(2)}$. These measures are given by

[TABLE]

with all density matrices now normalized to unity. Both quantities vanish if and only if $\rho_{\mathrm{(a)pcs}}^{(2)}=\rho^{(2)}$; hence a nonzero value quantifies the overall error (bounded by 1) of the PCS 2RDM. Figure 6 shows the very good agreement between the PCS 2RDMs and the exact 2RDM, as measured by the metrics (168) and (169), over the entire $U/(nJ)$ parameter space, from small ($N=50$) to large ($N=1000$) particle numbers. As shown by both metrics, $\rho_{\mathrm{pcs}}^{(2)}$ is very close to $\rho^{(2)}$ for all parameters considered. It is worth pointing out that in the extreme limits of zero and hard-core interactions [$U/(nJ)\rightarrow 0$ and $U/(nJ)\rightarrow\infty$], the exact PCS Ansatz gives a nearly exact prediction, but reaches its worst point in the transition region $U/(nJ)\sim 1$. In the transition region, the APCS performs just as well, with the trace distance (infidelity) showing imprecision of the order of $10^{-2}$ ($10^{-3}$); this region of good agreement gets wider as $N$ increases, but the APCS deviates from the PCS predictions for weak and strong interactions.

To understand the origin of the disagreements among the 2RDMs revealed in Fig. 6, we compare the parameters that define the exact and (A)PCS 2RDM matrix elements [see Eqs. (145) and (157)]. Note that for this purpose, we revert to setting the normalization of all 2RDMs equal to $N(N-1)$, as opposed to the 1 in Fig. 6. The PCS values of these parameters are given in Eqs. (185), (192), and (196), and the APCS values by Eqs. (159). The physical role played by these parameters is clear from looking at how they express the correlations in the site-label basis. From Eq. (152) the parameter $\beta$ expresses the first-order correlation function between the two sites; since it also gives the population difference of the Schmidt modes, we have already considered it in Fig. 4, so do not consider it further now. As one sees from Eqs. (158), the parameters $\gamma$ and $\delta$ give the second-order site correlation functions. We show in Fig. 7 a direct comparison between the parameters that define the exact and the PCS 2RDM matrix elements, together with their respective relative errors. We already have some evidence of the behavior of $\delta$ from its role in the ground-state energy (149), where it describes how the interaction strength $U$ competes with the tunneling $J$ through the transition from the mean-field regime to a Mott insulator. Indeed, it was the need for this contribution that motivated us to develop the more accurate, hybrid estimate of $\delta$ in Eq. (160). In the plots in Fig. 7, we use the more accurate, hybrid estimate,

[TABLE]

as the APCS result, in place of the large-$N$ estimate of zero from Eq. (159d). That $\gamma$ has the large-$N$ behavior proportional to $N^{2}$, whereas the leading-order behavior of $\delta$ is the Bogoliubov correction proportional to $N$ is immediately evident from the difference in scale of the ordinates in Fig. 7.

Finally, as a teaser to motivate using the PCS formalism to tackle many-body dynamics, we briefly investigate how faithfully PCS can represent the time evolution of the two-site BH model (137) compared to other ansätze. We consider a system of $N=1\,000$ particles initially in the ground state of the Hamiltonian (137) with $U/J=0$; this ground state is a single condensate mode that is symmetrically delocalized between the two sites. Then we suddenly quench the system by turning on an on-site repulsion $U/J=0.8$ [$U/(nJ)=0.0016$], which is in the intermediate regime of weak interactions. We calculate the resulting evolution of the exact state vector numerically using the Crank-Nicolson method; the exact 2RDM is derived from the state vector. We then optimize the parameter spaces of PCS, double-Fock state (DFS), and Gross-Pitaevskii state (GPS), i.e., a single-mode condensate, to minimize their trace distance to the exact 2RDM. The errors of the 2RDMs, as measured by the trace distances normalized to 1, given by PCS, DFS, and GPS, are plotted in Fig. 8 as functions of evolution time (in units of $\hbar/J$). The oscillation of the error given by GPS is a consequence of the collapse and revival of phase Greiner et al. (2002b); Schachenmayer et al. (2011); i.e., the purity of the 1RDM oscillates. These numerical results suggest that PCS might be useful in describing the dynamics in the strongly interacting regime.

VIII Conclusion and Future Work

In this paper we introduce the bosonic particle-correlated state (BPCS), a state of $N=l\times n$ bosons that is derived by symmetrizing the $n$-fold tensor product of an arbitrary $l$-boson (pure or mixed) state $\sigma^{(l)}$. We analyze in detail the pure-state case for $l=2$, i.e., $\sigma^{(2)}=\lvert\,\varPsi^{(2)}\,\rangle\langle\,\varPsi^{(2)}\,\rvert$, which we call a pair-correlated state (PCS). When there is just one Schmidt coefficient in the two-particle wave function, PCS reproduces the mean-field, Gross-Pitaevskii description of a single-mode condensate. When there is one dominant Schmidt coefficient in $\lvert\,\varPsi^{(2)}\,\rangle$, the leading-order corrections to just one Schmidt mode reproduce the particle-number-conserving Bogoliubov approximation. The allure of the PCS Ansatz lies in the case of many Schmidt coefficients of nearly the same size; in this situation, PCS describes a fragmented state that has large two-particle quantum correlations. At leading order in the large-$N$ limit, this PCS description is a sort of mean-field theory for a multimode condensate; the corrections to leading order describe Bogoliubov and higher-order modifications to the behavior of the multimode condensate. We provide methods for calculating the one- and two-particle reduced density matrices of PCS; from these RDMs come the predictions for physical observables such as mode populations and correlation functions.

As a test of the PCS Ansatz, we consider the two-site Bose-Hubbard model for the case of two identical lattice sites and large particle numbers. We calculate the ground-state energy of this model and analyze in detail the one- and two-particle reduced density matrices of the ground state. By comparing to exact results, we find that the PCS description, both in its exact form and in its large-$N$ approximate form, provides a remarkably good account of these ground-state properties. The main lesson of this analysis concerns the transition from a mean-field, single-mode condensate, in the regime where tunneling between sites dominates the on-site interaction, to a Mott insulator, in the regime where the on-site interaction dominates. This transition is well described in the large-$N$ limit by the PCS formalism. For this large-$N$ approximate PCS to give an accurate modeling, however, requires including in the approximate PCS ground-state energy not just Bogoliubov-order corrections, but also the first correction beyond Bogoliubov order. As a teaser, we also present results for how well the PCS Ansatz can match the dynamics of the symmetric two-site Bose-Hubbard model when it is suddenly quenched from being a single-mode condensate; our results suggest that the PCS formalism has the potential to describe this process quite well.

The success of the PCS formalism for the symmetric, two-site Bose-Hubbard model motivates further work. One question involves the ground state of an asymmetric Bose-Hubbard model, where comparison of Eqs. (143) and (145) shows that in the asymmetric case, the exact ground-state two-particle RDM does not have PCS form. Does the PCS formalism continue to provide a good description in this situation? A second question involves generalizing to the ground state of multi-site Bose-Hubbard Hamiltonians. The PCS description can incorporate exactly the strong-interaction Mott-insulator state of disconnected condensates at each site in the case of two sites, but it cannot do so for three or more sites. Does the PCS formalism continue to provide an accurate description of the transition from single condensate to Mott insulator for three or more sites? The question here might be put more generally as whether the correlations built into the two-particle RDM in the PCS description are the dominant correlations when there can be correlations among more than two sites.

The ultimate goal of our work on the PCS Ansatz is to analyze the dynamics of a Bose-Einstein condensate as one or more barriers is raised and lowered within it. In this situation, one must model how the “sites” within the BEC are evolving, i.e., how some appropriate set of single-particle wave functions are changing in time, and at the same time model how these sites are populated. The PCS formalism has the potential to capture this situation by providing a prescription for how the Schmidt orbitals are changing in time and by using the evolving two-particle PCS correlations to describe how the Schmidt orbitals are populated. Developing this description of PCS dynamics is the main goal of our future work.

Acknowledgements.

This work was supported by National Science Foundation Grants No. PHY-1314763 and No. PHY-1521016 and by Office of Naval Research Grant No. N00014-15-1-2167. A.B.T. acknowledges support from AFOSR Grant No. FA9550-12-1-0057. Z.J. acknowledges support from the NASA Advanced Exploration Systems program and NASA Ames Research Center.

Appendix A Bogoliubov Approximation to the Two-Site Bose-Hubbard Model

We start by considering the two-site Bose-Hubbard Hamiltonian in the Schmidt basis, which is given by Eq. (156). In the absence of interactions ($U=0$), the ground state of $\mathcal{H}_{\mathrm{tbh}}$ consists of a condensate, with all particles occupying the Schmidt mode $a_{1}=(b_{1}+b_{2})/\sqrt{2}$; the ground-state energy is $-JN$.

To derive the Bogoliubov approximation to $\mathcal{H}_{\mathrm{tbh}}$, we replace $a_{1}$ and $a_{1}^{\dagger}$ with $\sqrt{N}-a_{2}^{\dagger}a_{2}/2\sqrt{N}$ (so that $a_{1}^{\dagger}a_{1}$ is replaced by $N-a_{2}^{\dagger}a_{2}$) and keep terms up to quadratic order in $a_{2}$ and $a_{2}^{\dagger}$. This gives the Bogoliubov Hamiltonian, which in this case, is just a Hamiltonian for mode $a_{2}$:

[TABLE]

The Bogoliubov Hamiltonian (171) can be diagonalized by introducing the (squeezed) bosonic operators

[TABLE]

from which we get

[TABLE]

With the choice

[TABLE]

Eq. (171) takes the form

[TABLE]

Thus the Bogoliubov ground state has energy

[TABLE]

and a population in mode $a_{2}$ equal to

[TABLE]

The population imbalance of the Schmidt modes is given by

[TABLE]

The Bogoliubov approximation provides a very good account in the regime of very weak interactions, $NU/J\ll 1$, where the interaction energy is much smaller than the tunneling energy, but in this regime, the Bogoliubov corrections are really too tiny to worry about. The approximation should remain valid over a larger range of interaction strengths. Indeed, it should give a reasonably good account as long as the depletion is small, i.e., $\langle a_{2}^{\dagger}a_{2}\rangle\ll N$; notice that this region of validity includes regions where $NU/J\gg 1$. Taking the Bogoliubov approximation at face value for all values of interaction strength and assuming $N\gg 1$ and $NU/J\gg 1$, we can put the Bogoliubov results (177)–(179) in the simpler, approximate forms

[TABLE]

As noted above, the Bogoliubov approximation should be valid as long as $\langle a_{2}^{\dagger}a_{2}\rangle/N\ll 1$, i.e., $\sqrt{U/(NJ)}\ll 1$, which means that we should be able to rely on the expressions (180)–(182) as long as $NU/J\gg 1$ and $\sqrt{U/(NJ)}\ll 1$. This is the intermediate regime of weak interaction strength delimited by $1/N\ll U/J\ll N$. Notice now that Eqs. (180)–(182) are identical to the APCS predictions contained in Eqs. (166) and (167), which apply in this same intermediate regime of weak interaction strength. [In Figs. 4(c) and 5, the differences between APCS and the Bogoliubov approximation in this intermediate regime are accounted for by the fact that the differences are as small as the terms neglected in getting to Eqs. (166) and (167) and Eqs. (180)–(182).]

The crucial distinction between APCS and the Bogoliubov approximation is that APCS, by including post-Bogoliubov corrections, successfully navigates the transition from a superfluid to a Mott insulator as the interaction strength increases through $U/(NJ)\sim 1$, whereas the Bogoliubov approximation continues to make the predictions  (180)–(182) as $U$ increases, failing to notice the transition. The Bogoliubov approximation runs completely off the rails when the interaction strength increases to the point that $\langle a_{2}^{\dagger}a_{2}\rangle=N$, i.e., $U/(NJ)=32$, and cannot proceed to stronger interactions (the Bogoliubov expressions yield values for stronger interaction strengths, but they are unphysical). All these conclusions are satisfyingly consistent with our contention that the intermediate regime of weak interaction strength is where Bogoliubov corrections become important and to get beyond this regime and through the transition to a Mott insulator requires including post-Bogoliubov corrections, which the PCS formalism can accommodate.

Appendix B Exact PCS 2RDM for Schmidt Rank Two

For Schmidt rank $\nu=2$, we have in Eq. (49) the exact PCS normalization factor $\mathsf{N}$ in terms of a Gauss hypergeometric function. From this, we get, using Eq. (55),

[TABLE]

where

[TABLE]

denotes a ratio of contiguous Gauss hypergeometric functions. To get this result, we use, when taking the derivatives, the property

[TABLE]

which comes from the general identity,

[TABLE]

The population difference (185) is sufficient to calculate the PCS ground-state energy. In particular, from Eq. (121), we have

[TABLE]

and plugging this into the expression (149) for the ground-state energy, we get

[TABLE]

Equations (189) and (190) are the exact PCS analogues of the large-$N$, approximate expressions (160) and  (163). We find the exact PCS ground-state energy by minimizing $E_{\mathrm{pcs}}$ with respect to $z=\lambda_{2}^{2}/\lambda_{1}^{2}$; this minimum value of $z$ is then inserted in the expressions in this appendix to find the PCS ground-state 1RDM and 2RDM.

One can develop power-series expansions for the population difference (185) and the ground-state energy (190) in the single-condensate, mean-field regime, $\lambda_{2}^{2}/\lambda_{1}^{2}\ll 1$, and in the Mott regime, $|1-\lambda_{2}^{2}/\lambda_{1}^{2}|\ll 1$; these are, however, of limited value since they only provide small corrections to the mean-field and Mott limiting behaviors, corrections valid only in the extreme mean-field regime and the extreme Mott regime.

The normalization factor (49) also allows us to calculate the nonzero matrix elements of the PCS 2RDM, $\rho^{(2)}$. Using Eq. (57), we find the following diagonal elements:

[TABLE]

When taking the derivatives, we use the property

[TABLE]

which comes from applying a relation between contiguous hypergeometric functions to the general identity,

[TABLE]

There is, however, an easier way to find these 2RDM matrix elements than taking derivatives of the normalization factor. The second-to-last formulas on the right in Eqs. (191)–(193) express the procedure developed in Sec. VI.1 for determining the diagonal elements of the 2RDM from the 1RDM. In particular, these are Eqs. (128) and (130), and they can be used directly to find these 2RDM diagonal elements from the 1RDM. The final formula on the right relates the 2RDM matrix elements to the parameters in the Eq. (157).

Because $\rho_{12,12}^{(2)}=\rho_{21,21}^{(2)}=\rho_{12,21}^{(2)}=\rho_{21,12}^{(2)}$, the only nontrivial off-diagonal matrix elements remaining are $\rho^{(2)}_{11,22}=\rho^{(2)}_{22,11}$, which can be determined from Eq. (131) to be

[TABLE]

Notice that in accordance with our general conclusions, the entire 2RDM in the Schmidt basis is determined by the 1RDM and, in particular, in this case of $\nu=2$, by the population imbalance $\Delta_{\mathrm{pcs}}$ or, equivalently, by the hypergeometric ratio $\mathcal{F}$.

Appendix C Iterative Relations

In App. B, we demonstrated that the exact 2RDM of the PCS Ansatz for $\nu=2$ can be expressed using hypergeometric functions. It is very challenging, if not impossible, to generalize this analytical result to $\nu>2$. In this appendix, we derive an exact relation between the PCS 1RDMs of $2n$ and $2n+2$ particles using Wick’s theorem. This leads to an iterative algorithm that gives the exact 1RDM for $\nu\geq 2$ in time $O(N)$. Using Eqs. (128)–(131), one can calculate the 2RDM in time $O(1)$ from the 1RDM. In the large-$N$ limit, this iterative algorithm for the 1RDMs can be turned into a differential equation; a solution based on aseries expansion is given here. This procedure reduces the time complexity of the iterative algorithm for the 1RDMs to $O(1)$, at a price of introducing relative errors of order $O(1/N)$.

In the Schmidt basis, we have

[TABLE]

This relation implies that $\rho^{(1)}(n+1)$ is diagonalized in the Schmidt basis provided that $\rho^{(1)}(n)$ is diagonalized. Since $\rho^{(1)}(1)$ is diagonalized, mathematical induction allows us to conclude that $\rho^{(1)}(n)$ is diagonalized.

That $\rho^{(1)}(n)$ is diagonalized in the Schmidt basis can also be seen directly from the first line of Eq. (197): \big{(}\mathcal{A}^{\dagger}\big{)}^{n}\big{\lvert}\,\mathrm{vac}\,\big{\rangle} is a superposition of Fock states that have an even number of particles in each of the Schmidt single-particle states, so in the Fock-state superposition for a_{k}\,\big{(}\mathcal{A}^{\dagger}\big{)}^{n}\big{\lvert}\,\mathrm{vac}\,\big{\rangle}, the single-particle state $k$ always has an odd number of particles; thus a_{k}\,\big{(}\mathcal{A}^{\dagger}\big{)}^{n}\big{\lvert}\,\mathrm{vac}\,\big{\rangle} is orthogonal to a_{j}\,\big{(}\mathcal{A}^{\dagger}\big{)}^{n}\big{\lvert}\,\mathrm{vac}\,\big{\rangle} unless $j=k$.

For the diagonal elements $\rho^{(1)}_{jj}(n)\equiv\varrho_{j}(n)$, we have

[TABLE]

Using the normalization condition, $\sum_{j}\varrho_{j}(n)=2n$, we can write the ratio of the normalization factors, $\mathsf{N}_{n}/\mathsf{N}_{n+1}$, as

[TABLE]

Plugging this into Eq. (198), we get an iterative relation in which the normalization factors do not appear. We can also defer normalization to the end of the process, instead of imposing it at each iteration. The new iterative equation then reads

[TABLE]

where we replace the 1 with $(1/2n)\,\sum_{k}\varrho_{k}(n)$ to deal with the fact that $\varrho_{j}(n)$ is not normalized. The number of steps required in each iteration is proportional to $\nu$ and is independent of $n$, so the time complexity of calculating $\rho^{(1)}(n)$ using Eq. (198) or (200) is $O(\nu n)$. Because the PCS 2RDM can be calculated from the 1RDM (up to relative phase factors of the orbitals), searching for the lowest energy state in the PCS submanifold is exponentially faster than a brute-force approach considering the entire Hilbert space of $2n$ bosons.

We now let $n^{\prime}$ denote the iteration variable, and we suppose that we wish to iterate from $n^{\prime}=1$ to $n^{\prime}=n$. For sufficiently large $n$, all the probability concentrates on the dominant eigenvalues. To get useful results, we assume, as in Sec. V, that the differences of the $\lambda_{j}$s are small, and we use the parametrization (62) for the endpoint of the iteration, i.e., $\lambda_{j}^{2}=1+s_{j}/n$. To turn the iterative equation into a continuous differential equation, we introduce the parameter $\tau=n^{\prime}/n$. The new iterative equation then reads as follows:

[TABLE]

We now can write

[TABLE]

Taking the limit $n\rightarrow\infty$ yields the differential equation,

[TABLE]

where $\bar{\varrho}_{j}(\tau)=\varrho_{j}(n\tau)$. The problem with Eq. (203) is that it diverges at small $\tau$ due to the factor $1/2\tau$ unless $\bar{\varrho}_{j}(0)=0$ for all $j=1,2,\ldots,\nu$. This divergence is a consequence of our decision to defer the normalization to the end; one remedy is to modify the differential equation to

[TABLE]

where the extra term, which only introduces an overall factor, keeps $\bar{\varrho}_{j}$ from diverging for the initial condition $\bar{\varrho}_{j}(0)=1$, for $j=1,2,\ldots,\nu$. By mapping the iterative relation Eq. (200) to the differential equation (204), the PCS 1RDM can be solved approximately within time independent of $N$.

For the case $\nu=2$, we assume $s_{2}=-s_{1}$ without loss of generality. By introducing new variables, $\bar{\varrho}_{\pm}=(\bar{\varrho}_{1}\pm\bar{\varrho}_{2})/2$, we can write Eq. (204) as

[TABLE]

where $s_{-}=(s_{1}-s_{2})/2=s_{1}$. Taking derivatives with respect to $\tau$ on both sides of Eqs. (205) and (206) and manipulating the results, we have the decoupled second-order differential equations for $\bar{\varrho}_{\pm}$,

[TABLE]

These differential equations are solved by the zeroth- and first-order modified Bessel functions $I_{0}(\tau s_{-})$ and $I_{1}(\tau s_{-})$, respectively. Therefore, we recover our former results (96) and (97) using an entirely different approach.

For the general case $\nu\geq 2$, the solution to the linear differential (204) can be expressed using a Taylor series in $\tau$. This can be done most conveniently by introducing the transition matrix $T(\tau)$, such that

[TABLE]

where $\bar{\varrho}_{j}(\tau)$, the diagonal elements of the 1RDM, form a vector. Using the initial condition $\bar{\varrho}_{j}(0)=1$, for $j=1,2,\ldots,\nu$, we have

[TABLE]

It is convenient to choose the initial value of the transition matrix as $T(0)=M/\nu$ [any matrix that stabilizes the vector $(1,\,1,\ldots,1)^{T}$ suffices], where $M_{jk}=1$ (for $j,k=1,2,\ldots,\nu$) is the matrix of ones. The equation that governs the evolution of $T(\tau)$ can be derived from Eq. (204),

[TABLE]

where $S=\operatorname{diag}(s_{1},s_{2},\ldots,s_{\nu})$, and 𝟙 is the identity matrix. Suppose the transition matrix has a Taylor expansion,

[TABLE]

where the matrices $T_{\ell}$, for $\ell=1,2,\ldots$, are to be determined. Putting Eq. (212) into Eq. (211), we get

[TABLE]

The term of order $\tau^{-1}$ on the right hand side of Eq. (213) disappears with the choice $T(0)=M/\nu$. Comparing the coefficients of the terms $\tau^{\ell-1}$ on both sides of Eq. (213), we have

[TABLE]

By inverting the matrix \ell\hskip 0.4pt\mbox{\textbb{1}}+\big{(}\nu\mbox{\textbb{1}}-M\big{)}/2, we have

[TABLE]

The $\ell$th order matrix $T_{\ell}$ can be solved iteratively by applying the above relation $\ell$ times to $T_{0}=M/\nu$. To simply this procedure, we notice that

[TABLE]

holds for any diagonal matrix $D$, including powers of $S$. Therefore, the solution to Eq. (215) takes the form $T_{\ell}=P_{\ell}(S)\hskip 0.4ptM/\nu$, where $P_{\ell}(S)$ is some polynomial of $S$ of order $\ell$, e.g., $P_{1}(S)=S/(1+\nu/2)$ and $P_{2}(S)=[S^{2}+\operatorname{tr}(S^{2})/4]/[1+\nu/2)(2+\nu/2)]$. The matrix norm of $T_{\ell}$ begins to fall quickly after $\ell>s_{1}$, which gives an estimate on how many terms are needed in the expansion to get a desired precision.

In this appendix, we discussed an iterative algorithm to calculate the 1RDMs of the PCS Ansatz exactly in time $O(N)$. For large $N$, the iterative steps can be approximately mapped to a differential equation. This differential equation gives the same result as those derived previously for $\nu=2$. For $\nu\geq 2$, the solution to the differential equations can be expressed as a Taylor expansion. Compared to our other approaches, the iterative approach is particularly suitable for numerical evaluations of the PCS Ansatz.