Solvable Model of a Generic Trapped Mixture of Interacting Bosons: Reduced Density Matrices and Proof of Bose-Einstein Condensation

TL;DR

This paper presents an exactly solvable model of a two-species bosonic mixture in a harmonic trap, analytically proving Bose-Einstein condensation and properties of reduced density matrices in the infinite-particle limit.

Contribution

It provides an analytical proof of Bose-Einstein condensation and reduced density matrix properties for a solvable two-species bosonic mixture model.

Findings

Both species are 100% condensed in the infinite-particle limit.

Inter-species reduced density matrix is separable and product form.

Energy and densities match Gross-Pitaevskii predictions in the limit.

Abstract

A mixture of two kinds of identical bosons, species $1$ and species $2$, held in a harmonic potential and interacting by harmonic intra-species and inter-species particle-particle interactions is discussed. To prove Bose-Einstein condensation of the mixture three steps are needed. First, we integrate the all-particle density matrix, employing a four-parameter matrix recurrence relations, down to the lowest-order intra-species and inter-species reduced density matrices of the mixture. Second, the coupled Gross-Pitaevskii (mean-field) equations of the mixture are solved analytically. Third, we analyze the mixture's reduced density matrices in the limit of an infinite number of particles of both species $1$ and $2$ (when the interaction parameters, i.e., the products of the number of particles times the intra-species and inter-species interaction strengths, are held fixed) and prove that: (i) Both species $1$ and $2$ are 100\% condensed; (ii) The inter-species reduced density matrix per particle is separable and given by the product of the intra-species reduced density matrices per particle; and (iii) The mixture's energy per particle, and reduced density matrices and densities per particle all coincide with the Gross-Pitaevskii quantities. Finally, when the infinite-particle limit is taken with respect to, say, species $1$ only (with interaction parameters held fixed) we prove that: (iv) Only species $1$ is 100\% condensed and its reduced density matrix and density per particle, as well as the mixture's energy per particle, coincide with the Gross-Pitaevskii quantities of species $1$ alone; and (v) The inter-species reduced density matrix per particle is nonetheless separable and given by the product of the intra-species reduced density matrices per particle.

Equations

$$\hat{H}({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{N_1},{\mathbf{y}}_1,\dots ,{\mathbf{y}}_{N_2})=$$

$$=\sum _{j=1}^{N_1}\left(-\frac{1}{2m_1}\frac{{\partial }^2}{\partial {\mathbf{x}}_j^2}+\frac{1}{2}{m}_1{\omega }_1^2{\mathbf{x}}_j^2\right)+\sum _{j=1}^{N_2}\left(-\frac{1}{2m_2}\frac{{\partial }^2}{\partial {\mathbf{y}}_j^2}+\frac{1}{2}{m}_2{\omega }_2^2{\mathbf{y}}_j^2\right)+$$

$$\displaystyle\quad+\lambda_{1}\sum_{1\leq j<k}^{N_{1}}({\bf x}_{j}-{\bf x}_{k})^{2}+\lambda_{2}\sum_{1\leq j<k}^{N_{2}}({\bf y}_{j}-{\bf y}_{k})^{2}+\lambda_{12}\sum_{j=1}^{N_{1}}\sum_{k=1}^{N_{2}}({\bf x}_{j}-{\bf y}_{k})^{2}.\$$

$${\mathbf{Q}}_k=\frac{1}{\sqrt{k(k+1)}}\sum _{j=1}^k({\mathbf{x}}_{k+1}-{\mathbf{x}}_j),1\le k\le N_1-1,$$

$${\mathbf{Q}}_{N_1-1+k}=\frac{1}{\sqrt{k(k+1)}}\sum _{j=1}^k({\mathbf{y}}_{k+1}-{\mathbf{y}}_j),1\le k\le N_2-1,$$

$${\mathbf{Q}}_{N-1}=\sqrt{\frac{N_2}{N_1}}\sum _{j=1}^{N_1}{\mathbf{x}}_j-\sqrt{\frac{N_1}{N_2}}\sum _{j=1}^{N_2}{\mathbf{y}}_j,$$

$$\displaystyle{\bf Q}_{N}=\frac{m_{1}}{M}\sum_{j=1}^{N_{1}}{\bf x}_{j}+\frac{m_{2}}{M}\sum_{j=1}^{N_{2}}{\bf y}_{j}.\$$

$$\hat{H}({\mathbf{Q}}_1,\dots ,{\mathbf{Q}}_N)=\sum _{k=1}^{N_1-1}\left(-\frac{1}{2m_1}\frac{{\partial }^2}{\partial {\mathbf{Q}}_k^2}+\frac{1}{2}{m}_1{\Omega }_1^2{\mathbf{Q}}_k^2\right)+$$

$$+\sum _{k=N_1}^{N-2}\left(-\frac{1}{2m_2}\frac{{\partial }^2}{\partial {\mathbf{Q}}_k^2}+\frac{1}{2}{m}_2{\Omega }_2^2{\mathbf{Q}}_k^2\right)+\left(-\frac{1}{2M_{12}}\frac{{\partial }^2}{\partial {\mathbf{Q}}_{N-1}^2}+\frac{1}{2}{M}_{12}{\Omega }_{12}^2{\mathbf{Q}}_{N-1}^2\right)+$$

$$\displaystyle\quad+\left(-\frac{1}{2M}\frac{\partial^{2}}{\partial{\bf Q}_{N}^{2}}+\frac{1}{2}M\omega^{2}{\bf Q}_{N}^{2}\right),\qquad M_{12}=\frac{m_{1}m_{2}}{M},\ M=N_{1}m_{1}+N_{2}m_{2}.\$$

$${\Omega }_1=\sqrt{{\omega }^2+\frac{2}{m_1}(N_1{\lambda }_1+N_2{\lambda }_{12})},{\Omega }_2=\sqrt{{\omega }^2+\frac{2}{m_2}(N_2{\lambda }_2+N_1{\lambda }_{12})},$$

$$\displaystyle\Omega_{12}=\sqrt{\omega^{2}+2\left(\frac{N_{1}}{m_{2}}+\frac{N_{2}}{m_{1}}\right)\lambda_{12}}=\sqrt{\omega^{2}+\frac{2\lambda_{12}}{M_{12}}},\qquad\omega.\$$

$${\Omega }_{12}^2={\omega }^2+\frac{2{\lambda }_{12}}{M_{12}}>0⟹{\lambda }_{12}>-\frac{M_{12}{\omega }^2}{2},$$

$${\Omega }_1^2={\omega }^2+\frac{2}{m_1}(N_1{\lambda }_1+N_2{\lambda }_{12})>0⟹{\lambda }_1>-m_1\frac{N_2}{N_1}{\lambda }_{12}-\frac{m_1{\omega }^2}{2N_1},$$

$$\displaystyle\Omega_{2}^{2}=\omega^{2}+\frac{2}{m_{2}}(N_{2}\lambda_{2}+N_{1}\lambda_{12})>0\quad\Longrightarrow\quad\lambda_{2}>-m_{2}\frac{N_{1}}{N_{2}}\lambda_{12}-\frac{m_{2}\omega^{2}}{2N_{2}}.\$$

$$\Psi ({\mathbf{Q}}_1,\dots ,{\mathbf{Q}}_N)={\left(\frac{m_1{\Omega }_1}{\pi }\right)}^{\frac{3(N_1-1)}{4}}{\left(\frac{m_2{\Omega }_2}{\pi }\right)}^{\frac{3(N_2-1)}{4}}{\left(\frac{M_{12}{\Omega }_{12}}{\pi }\right)}^{\frac{3}{4}}{\left(\frac{M\omega }{\pi }\right)}^{\frac{3}{4}}\times$$

$$\displaystyle\quad\times e^{-\frac{1}{2}\left(m_{1}\Omega_{1}\sum_{k=1}^{N_{1}-1}{\bf Q}_{k}^{2}+m_{2}\Omega_{2}\sum_{k=N_{1}}^{N-2}{\bf Q}_{k}^{2}+M_{12}\Omega_{12}{\bf Q}_{N-1}^{2}+M\omega{\bf Q}_{N}^{2}\right)},\$$

$$\displaystyle E=\frac{3}{2}\left[(N_{1}-1)\Omega_{1}+(N_{2}-1)\Omega_{2}+\Omega_{12}+\omega\right]=\frac{3}{2}\Bigg{[}(N_{1}-1)\sqrt{\omega^{2}+\frac{2}{m_{1}}(N_{1}\lambda_{1}+N_{2}\lambda_{12})}+$$

$$\displaystyle\quad+(N_{2}-1)\sqrt{\omega^{2}+\frac{2}{m_{2}}(N_{2}\lambda_{2}+N_{1}\lambda_{12})}+\sqrt{\omega^{2}+\frac{2\lambda_{12}}{M_{12}}}+\omega\Bigg{]}\$$

$$\Psi ({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{N_1},{\mathbf{y}}_1,\dots ,{\mathbf{y}}_{N_2})={\left(\frac{m_1{\Omega }_1}{\pi }\right)}^{\frac{3(N_1-1)}{4}}{\left(\frac{m_2{\Omega }_2}{\pi }\right)}^{\frac{3(N_2-1)}{4}}{\left(\frac{M_{12}{\Omega }_{12}}{\pi }\right)}^{\frac{3}{4}}{\left(\frac{M\omega }{\pi }\right)}^{\frac{3}{4}}\times$$

$$\displaystyle\quad\times e^{-\frac{\alpha_{1}}{2}\sum_{j=1}^{N_{1}}{\bf x}_{j}^{2}-\beta_{1}\sum_{1\leq j<k}^{N_{1}}{\bf x}_{j}\cdot{\bf x}_{k}}\times e^{-\frac{\alpha_{2}}{2}\sum_{j=1}^{N_{2}}{\bf y}_{j}^{2}-\beta_{2}\sum_{1\leq j<k}^{N_{2}}{\bf y}_{j}\cdot{\bf y}_{k}}\times e^{+\gamma\sum_{j=1}^{N_{1}}\sum_{k=1}^{N_{2}}{\bf x}_{j}\cdot{\bf y}_{k}},\$$

$${\alpha }_1=m_1\left[{\Omega }_1\left(1-\frac{1}{N_1}\right)+(m_2{N}_2{\Omega }_{12}+m_1{N}_1\omega )\frac{1}{M{N}_1}\right]=m_1{\Omega }_1+{\beta }_1,$$

$${\beta }_1=m_1\left[-{\Omega }_1\frac{1}{N_1}+(m_2{N}_2{\Omega }_{12}+m_1{N}_1\omega )\frac{1}{M{N}_1}\right],$$

$${\alpha }_2=m_2\left[{\Omega }_2\left(1-\frac{1}{N_2}\right)+(m_1{N}_1{\Omega }_{12}+m_2{N}_2\omega )\frac{1}{M{N}_2}\right]=m_2{\Omega }_2+{\beta }_2,$$

$${\beta }_2=m_2\left[-{\Omega }_2\frac{1}{N_2}+(m_1{N}_1{\Omega }_{12}+m_2{N}_2\omega )\frac{1}{M{N}_2}\right],$$

$$\displaystyle\gamma=\frac{m_{1}m_{2}}{M}(\Omega_{12}-\omega)=M_{12}(\Omega_{12}-\omega).\$$

$$\Psi ({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{N_1},{\mathbf{y}}_1,\dots ,{\mathbf{y}}_{N_2}){\Psi }^{\ast }({\mathbf{x}}_1^{\prime },\dots ,{\mathbf{x}}_{N_1}^{\prime },{\mathbf{y}}_1^{\prime },\dots ,{\mathbf{y}}_{N_2}^{\prime })=$$

$${\left(\frac{m_1{\Omega }_1}{\pi }\right)}^{\frac{3(N_1-1)}{2}}{\left(\frac{m_2{\Omega }_2}{\pi }\right)}^{\frac{3(N_2-1)}{2}}{\left(\frac{M_{12}{\Omega }_{12}}{\pi }\right)}^{\frac{3}{2}}{\left(\frac{M\omega }{\pi }\right)}^{\frac{3}{2}}\times$$

$$\times e^{-\frac{{\alpha }_1}{2}{\sum }_{j=1}^{N_1}({\mathbf{x}}_j^2+{{\mathbf{x}}_j^{\prime }}^2)-{\beta }_1{\sum }_{1\le j<k}^{N_1}({\mathbf{x}}_j\cdot {\mathbf{x}}_k+{\mathbf{x}}_j^{\prime }\cdot {\mathbf{x}}_k^{\prime })}\times$$

$$\times e^{-\frac{{\alpha }_2}{2}{\sum }_{j=1}^{N_2}({\mathbf{y}}_j^2+{{\mathbf{y}}_j^{\prime }}^2)-{\beta }_2{\sum }_{1\le j<k}^{N_2}({\mathbf{y}}_j\cdot {\mathbf{y}}_k+{\mathbf{y}}_j^{\prime }\cdot {\mathbf{y}}_k^{\prime })}\times$$

$$\displaystyle\quad\times e^{+\gamma\sum_{j=1}^{N_{1}}\sum_{k=1}^{N_{2}}({\bf x}_{j}\cdot{\bf y}_{k}+{\bf x}^{\prime}_{j}\cdot{\bf y}^{\prime}_{k})},\$$

$${\rho }_1(\mathbf{x},{\mathbf{x}}^{\prime })=N_1\int d{\mathbf{x}}_2\cdots d{\mathbf{x}}_{N_1}d{\mathbf{y}}_1\cdots d{\mathbf{y}}_{N_2}\Psi (\mathbf{x},{\mathbf{x}}_2,\dots ,{\mathbf{x}}_{N_1},{\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_{N_2})\times$$

$$\times {\Psi }^{\ast }({\mathbf{x}}^{\prime },{\mathbf{x}}_2,\dots ,{\mathbf{x}}_{N_1},{\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_{N_2}),$$

$${\rho }_2(\mathbf{y},{\mathbf{y}}^{\prime })=N_2\int d{\mathbf{x}}_1\cdots d{\mathbf{x}}_{N_1}d{\mathbf{y}}_2\cdots d{\mathbf{y}}_{N_2}\Psi ({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_{N_1},\mathbf{y},{\mathbf{y}}_2,\dots ,{\mathbf{y}}_{N_2})\times$$

$$\displaystyle\qquad\qquad\qquad\times\Psi^{\ast}({\bf x}_{1},{\bf x}_{2},\ldots,{\bf x}_{N_{1}},{\bf y}^{\prime},{\bf y}_{2},\ldots,{\bf y}_{N_{2}}),\$$

$${\rho }_{12}(\mathbf{x},{\mathbf{x}}^{\prime },\mathbf{y},{\mathbf{y}}^{\prime })=N_1{N}_2\int d{\mathbf{x}}_2\cdots d{\mathbf{x}}_{N_1}d{\mathbf{y}}_2\cdots d{\mathbf{y}}_{N_2}\Psi (\mathbf{x},{\mathbf{x}}_2,\dots ,{\mathbf{x}}_{N_1},\mathbf{y},{\mathbf{y}}_2,\dots ,{\mathbf{y}}_{N_2})\times$$

$$\displaystyle\qquad\qquad\qquad\qquad\times\Psi^{\ast}({\bf x}^{\prime},{\bf x}_{2},\ldots,{\bf x}_{N_{1}},{\bf y}^{\prime},{\bf y}_{2},\ldots,{\bf y}_{N_{2}}),\$$

$$F_{N_1,N_2}({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{N_1},{\mathbf{y}}_1,\dots ,{\mathbf{y}}_{N_2};{\mathbf{x}}_1^{\prime },\dots ,{\mathbf{x}}_{N_1}^{\prime },{\mathbf{y}}_1^{\prime },\dots ,{\mathbf{y}}_{N_2}^{\prime };$$

$${\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2,C_{N_1,N_2},C_{N_1,N_2}^{\prime },D_{N_1,N_2},D_{N_1,N_2}^{\prime })=$$

$$=e^{-\frac{{\alpha }_1}{2}{\sum }_{j=1}^{N_1}({\mathbf{x}}_j^2+{{\mathbf{x}}_j^{\prime }}^2)-{\beta }_1{\sum }_{1\le j<k}^{N_1}({\mathbf{x}}_j\cdot {\mathbf{x}}_k+{\mathbf{x}}_j^{\prime }\cdot {\mathbf{x}}_k^{\prime })}{e}^{-\frac{{\alpha }_2}{2}{\sum }_{j=1}^{N_2}({\mathbf{y}}_j^2+{{\mathbf{y}}_j^{\prime }}^2)-{\beta }_2{\sum }_{1\le j<k}^{N_2}({\mathbf{y}}_j\cdot {\mathbf{y}}_k+{\mathbf{y}}_j^{\prime }\cdot {\mathbf{y}}_k^{\prime })}\times$$

$$\times e^{-\frac{1}{4}{C}_{N_1,N_2}({\mathbf{X}}_{N_1}+{\mathbf{X}}_{N_1}^{\prime }{)}^2}{e}^{-\frac{1}{4}{C}_{N_1,N_2}^{\prime }({\mathbf{Y}}_{N_2}+{\mathbf{Y}}_{N_2}^{\prime }{)}^2}\times$$

$$\displaystyle\quad\times e^{+\frac{1}{2}D_{N_{1},N_{2}}({\bf X}_{N_{1}}+{\bf X}^{\prime}_{N_{1}})({\bf Y}_{N_{2}}+{\bf Y}^{\prime}_{N_{2}})}e^{+\frac{1}{2}D^{\prime}_{N_{1},N_{2}}({\bf X}_{N_{1}}-{\bf X}^{\prime}_{N_{1}})({\bf Y}_{N_{2}}-{\bf Y}^{\prime}_{N_{2}})},\$$

$${\mathbf{X}}_{N_1}=\sum _{j=1}^{N_1}{\mathbf{x}}_j,{\mathbf{X}}_{N_1}^{\prime }=\sum _{j=1}^{N_1}{\mathbf{x}}_j^{\prime },{\mathbf{Y}}_{N_2}=\sum _{j=1}^{N_2}{\mathbf{y}}_j,{\mathbf{Y}}_{N_2}^{\prime }=\sum _{j=1}^{N_2}{\mathbf{y}}_j^{\prime }$$

$$C_{N_1,N_2}=0,C_{N_1,N_2}^{\prime }=0,D_{N_1,N_2}=\gamma ,D_{N_1,N_2}^{\prime }=\gamma .$$

$$\begin{array}[]{ccccc}F_{N_{1},N_{2}}&\Longrightarrow&F_{N_{1},1}&\Rightarrow&F_{N_{1},0}\\ \big{\Downarrow}&&\big{\Downarrow}&&\big{\Downarrow}\\ F_{1,N_{2}}&\Longrightarrow&F_{1,1}&&F_{1,0}\\ \Downarrow&&&&\\ F_{0,N_{2}}&\Longrightarrow&F_{0,1}&&\end{array}\ .$$

Full text

I Introduction

Mixtures of Bose-Einstein condensates made of ultra-cold quantum gases and their physical properties have attracted much attention in the past twenty years, see for example B1 ; B2 ; B3 ; B4 ; B5 ; B6 ; B7 ; B8 ; B9 ; B10 ; B11 ; B12 ; B13 ; B14 ; B15 ; B16 ; B17 ; B18 ; B19 ; B20 ; B21 ; B22 ; B23 ; B24 ; B25 ; B26 . The ground state of trapped mixtures has been studied both at the mean-field level (within Gross-Pitaevskii theory) and by different many-body theoretical and numerical tools. Unlike their older sibling – Bose-Einstein condensates made of a single species – for which there have been rigorous results connecting in the infinite-particle limit the many-body and mean-field expressions for the energy per particle, density per particle, and providing proof of $100\%$ Bose-Einstein condensation in the ground state Yngvason_PRA ; Lieb_PRL , there are to the best of our knowledge no such results for trapped mixtures.

To address this topic, we present an exactly-solvable model of a generic mixture of trapped interacting bosons which allows one, for a start, to compute analytically the energy and wave-function of the mixture for any number of particles. From there, and with some effort as we shall see below, intra-species and inter-species reduced density matrices are computed and analyzed. This allows us to get concrete results in the infinite-particle limit (when the products of the number of particles times the intra-species and inter-species interaction strengths are held fixed) on the energy per particle, intra-species and inter-species reduced density matrices per particle, and to prove Bose-Einstein condensate of each of the species. By solving analytically for the ground state at mean-field level (within Gross-Pitaevskii theory), we are also able to compare in the infinite-particle limit the many-body and mean-field results. For a mixture there is a nice twist on the infinite-particle limit, which may be taken with respect to both species or with respect to one of them. The conclusions from both infinite-particle-limit procedures are compared and contrasted.

The model we present for the mixture is that of two species of indistinguishable bosons, species $1$ with $N_{1}$ bosons of mass $m_{1}$ and species $2$ with $N_{2}$ bosons of mass $m_{2}$, held in a harmonic potential of frequency $\omega$ and interacting by harmonic intra-species and inter-species particle-particle interactions of strengths $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{12}$, respectively. This is the harmonic-interaction model for trapped bosonic mixtures. The harmonic-interaction model for identical particles has been widely used for single-species bosons HIM_Cohen ; HIM_Yan ; HIM_Po1 ; HIM_Benchmarks ; HIM_Po0 ; Schilling_FER_BOS ; Ental_FER_BOS , fermions HIM_Po0 ; Schilling_FER_BOS ; Ental_FER_BOS ; Schilling_FER_HIM_WORK ; Axel_MCTDHF_HIM , and un-trapped (i.e., translationally-invariant) bosonic mixtures HIM_MIX_JPA1 ; HIM_MIX_JPA2 ; HIM_MIX_PRE . Most recently, we were able to solve a specific case of the harmonic-interaction model for trapped mixtures, a trapped symmetric mixture for which $N_{1}=N_{2}$, $m_{1}=m_{2}$, and $\lambda_{1}=\lambda_{2}$, and showed that the energy per particle and densities per particle (i.e., the diagonal of reduced density matrices per particle) converge in the infinite-particle limit to their mean-field analogs BB_HIM_SYM . To treat the general mixture with unequal numbers of particles, masses, and interaction strengths, and especially to derive its reduced density matrices and prove Bose-Einstein condensation, it is needed to generalize the analytical techniques developed by Cohen and Lee HIM_Cohen for single-species bosons and extended in BB_HIM_SYM for the symmetric mixture substantially further, which is done below.

The structure of the paper is as follows. In Sec. II we present the harmonic-interaction model for a generic trapped mixture of bosons and discuss its ground-state energy and wave-function. In Sec. III we construct explicitly the lowest-order intra-species and inter-species reduced density matrices of the mixture, solve analytically the two-coupled Gross-Pitaevskii equations, perform the infinite-particle limit in which the many-body and mean-field results are compared, and prove Bose-Einstein condensation. Concluding remarks are provided in Sec. IV. Further details of the derivations are collected in the appendixes.

II The Harmonic-Interaction Model for Trapped Bosonic Mixtures

Consider a mixture of two distinguishable types of identical bosons which we label $1$ and $2$. The bosons are trapped in a three-dimensional isotropic harmonic potential of frequency $\omega$ and interact via harmonic particle-particle interactions. In the present work we only deal with the ground state of the trapped mixture. We treat the case of a generic mixture. Namely, a mixture consisting of $N_{1}$ bosons of type $1$ and mass $m_{1}$ and $N_{2}$ bosons of type $2$ and mass $m_{2}$. The total number of particles is denoted by $N=N_{1}+N_{2}$. Furthermore, the two intra-species interaction strengths are denoted by $\lambda_{1}$ and $\lambda_{2}$, and the inter-species interaction strength by $\lambda_{12}$. Positive values of $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{12}$ mean attractive particle-particle interactions whereas negative values imply repulsive interactions REMARK_INTERACTIONS .

The Hamiltonian of the mixture is then given by ($\hbar=1$)

[TABLE]

Here, the coordinates ${\bf x}_{j}$ denote bosons of type $1$ and ${\bf y}_{k}$ bosons of type $2$. We work in Cartesian coordinates where the vector ${\bf x}=(x_{1},x_{2},x_{3})$ denotes the position of a boson of type $1$ in three dimensions, and $\frac{1}{i}\frac{\partial}{\partial{\bf x}}=\frac{1}{i}\left(\frac{\partial}{\partial x_{1}},\frac{\partial}{\partial x_{2}},\frac{\partial}{\partial x_{3}}\right)$ its momentum. To avoid cumbersome notation, we denote ${\bf x}^{2}\equiv x_{1}^{2}+x_{2}^{2}+x_{3}^{2}$ and $\frac{\partial^{2}}{\partial{\bf x}^{2}}\equiv\frac{\partial^{2}}{\partial x_{1}^{2}}+\frac{\partial^{2}}{\partial x_{2}^{2}}+\frac{\partial^{2}}{\partial x_{3}^{2}}$. Analogous notation is employed for the vector ${\bf y}$ of a boson of type $2$.

To diagonalize the Hamiltonian (II) we transform to the Jacoby coordinates HIM_MIX_PRE

[TABLE]

The meaning of (II) is as follows: The first group of $N_{1}-1$ coordinates are relative coordinates of the bosons of species $1$; the second group of $N_{2}-1$ coordinates are relative coordinates the bosons of species $2$; ${\bf Q}_{N-1}$ can be seen as a relative coordinate between the center-of-mass of the species $1$ and the center-of-mass of the species $2$ bosons; and, finally, ${\bf Q}_{N}$ is the center-of-mass coordinate of all particles in the mixture.

Using (II) and (A)-(A) in Appendix A, the Hamiltonian (II) transforms to the diagonal form

[TABLE]

The transformed Hamiltonian of the mixture (II) is that of $N$ uncoupled harmonic oscillators having the four masses $m_{1}$, $m_{2}$, $M_{12}$ (relative mass of the mixture), and $M$ (total mass of the mixture), and four frequencies

[TABLE]

The multiplicity of the frequencies is $N_{1}-1$, $N_{2}-1$, $1$, and $1$, respectively, corresponding to the types and numbers of Jabobi coordinates REMARK_DEGEN .

The frequencies (II) must be positive in order for a bound solution to exist. This dictates bounds on both the intra-species $\lambda_{1}$ and $\lambda_{2}$ and inter-species $\lambda_{12}$ interactions which are:

[TABLE]

The meaning of these bounds are as follows: The inter-species interaction $\lambda_{12}$ is bounded from below by the frequency of the trap and the relative mass, irrespective of the intra-species interactions $\lambda_{1}$ and $\lambda_{2}$, otherwise the mixture cannot be trapped in the harmonic potential. On the other hand, the intra-species interaction $\lambda_{1}$ is limited by the chosen inter-species interaction $\lambda_{12}$, the numbers of particles, and the mass $m_{1}$, and analogously $\lambda_{2}$.

We can now proceed and prescribe the normalized ground-state wave-function

[TABLE]

along with the ground-state energy

[TABLE]

of the trapped mixture REMARK_NO_TRAP . Using the bounds for $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{12}$ in (II), we obtain that the ground-state energy of the mixture is bound from below by $E>\frac{3}{2}\omega$, which is obtained for $\Omega_{1}\to 0^{+}$, $\Omega_{2}\to 0^{+}$, and $\Omega_{12}\to 0^{+}$. This means that all relative degrees of freedom are marginally bound, and essentially only the center-of-mass degree of freedom is bound in the harmonic trap. The system is then predominantly repulsive. At the other end, when the mixture is predominantly attractive, the energy is unbound from above.

To express the wave-function with respect to the original spatial coordinates we use relations (A) and find

[TABLE]

where the parameters are

[TABLE]

The ground-state wave-function of the mixture (II) is seen to be comprised of a product of a type $1$ boson part, a type $2$ boson part, and a coupling $1$–$2$ part. As required by indistinguishability of identical bosons, $\Psi$ is symmetric to permutation of the coordinates ${\bf x}_{j}$ and ${\bf x}_{k}$ of any two $1$ type bosons, and likewise symmetric to permutation of the coordinates ${\bf y}_{j}$ and ${\bf y}_{k}$ of any two $2$ type bosons, but is not symmetric to permutation of the distinguishable $1$ and $2$ type bosons.

The harmonic-interaction model for a generic trapped mixture of interacting bosons presented above admits a wealth of properties that can all be studied in principle analytically. After all, the ground-state wave-function and energy are given as explicit simple functions of all parameters – masses $m_{1}$ and $m_{2}$, numbers of particles $N_{1}$ and $N_{2}$, interactions strengths $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{12}$, and trapping frequency $\omega$. All quantum properties of the mixture can be computed from its wave-function, although, as we shall see in the next section, not necessarily without some effort.

III Reduced Density Matrices and Proof of Bose-Einstein Condensation

To show that the individual species are Bose-Einstein condensed three steps are required. The first, the computation of the intra-species reduced one-particle density matrices. To discuss the complementary question of separability requires the inter-species reduced two-body density matrix. These are calculated analytically in Sec. III.1 as a function of the masses, interaction strengths, and numbers of particles. The second step is to solve, again analytically, the system at the mean-field level. This means finding the solution to the two-species coupled Gross-Pitaevskii equations. This is performed in Sec. III.2. The third and final step is to perform the infinite-particle limit and to compare and contrast the exact and mean-field solutions in this limit, which is carried out in Sec. III.3. In a mixture one can perform the infinite-particle limit with respect to the two species or with respect to one of the species. We will discuss and compare both infinite-particle-limit procedures below.

III.1 The intra-species and inter-species reduced density matrices

We start from the full $N$-particle density matrix of the mixture,

[TABLE]

which for convenience is here normalized to one, $\int d{\bf x}_{1}\cdots d{\bf x}_{N_{1}}d{\bf y}_{1}\cdots d{\bf y}_{N_{2}}|\Psi|^{2}=1$. The intra-species reduced one-body density matrices of the mixture,

[TABLE]

are given by integrating over $N-1$ coordinates, and the inter-species reduced two-body density matrix

[TABLE]

which is the lowest-order inter-species reduced density matrix of the mixture, is obtained by integrating over $N-2$ coordinates.

To obtain the lowest-order reduced density matrices (III.1) and (III.1) we need to perform multiple integrations of the $N$-particle density matrix of the mixture (III.1). Moreover, the latter contains a coupling $1$–$2$ part because of the inter-species interaction. Thus, without an appropriate construction, the task becomes quickly impractical with increasing $N_{1}$ and $N_{2}$. Below, to perform the integration, we derive a four-parameter matrix-recurrence relations, thus generalizing the one-parameter and two-parameter vector-recurrence relations put forward, respectively, in the cases of the single-species HIM_Cohen and symmetric-mixture BB_HIM_SYM harmonic-interaction models.

To integrate the $N$-particle density matrix we begin by introducing the auxiliary function

[TABLE]

where $\alpha_{1}$, $\beta_{1}$, $\alpha_{2}$, and $\beta_{2}$, are given in (II),

[TABLE]

are vectors REMARK_VECTORS , and $C_{N_{1},N_{2}}$, $C^{\prime}_{N_{1},N_{2}}$, $D_{N_{1},N_{2}}$, and $D^{\prime}_{N_{1},N_{2}}$ are constants to play a further role below. We note that if one were to equate the $N$-particle density (III.1) and the auxiliary function $F_{N_{1},N_{2}}$ then this would imply that

[TABLE]

We will perform the integration of the auxiliary function $F_{N_{1},N_{2}}$ in two steps, first by integrating the ${\bf y}_{N_{2}},{\bf y}^{\prime}_{N_{2}}={\bf y}_{N_{2}}$ variables and then the ${\bf x}_{N_{1}},{\bf x}^{\prime}_{N_{1}}={\bf x}_{N_{1}}$ variables. We shall call the first step the horizontal reduction of the auxiliary function $F_{N_{1},N_{2}}$ and the second step the vertical reduction. The sequences of integrations to be performed below may be written symbolically as

[TABLE]

The end terms of these integrations, $F_{1,1}$ and $F_{1,0}$, $F_{0,1}$, will be connected with the reduced density matrices (III.1) and (III.1), respectively.

Thus, to perform the multiple integration steps in the horizontal reduction (16), we seek for a recurrence relation and write

[TABLE]

where the Gaussian integral $\int d{\bf y}e^{-a{\bf y}^{2}-2b{\bf y}\cdot{\bf X}}=\left(\frac{\pi}{a}\right)^{\frac{3}{2}}e^{\frac{b^{2}{\bf X}^{2}}{a}}$ is used and $\int d{\bf y}_{N_{2}}F_{N_{1},N_{2}}\ldots$ implicitly implies that ${\bf y}^{\prime}_{N_{2}}={\bf y}_{N_{2}}$ is used in the integrand. The equality (III.1) relates the auxiliary function $F_{N_{1},N_{2}}$ with $N_{1}+N_{2}$ coordinates and constants $C_{N_{1},N_{2}}$, $C^{\prime}_{N_{1},N_{2}}$, $D_{N_{1},N_{2}}$, and $D^{\prime}_{N_{1},N_{2}}$, to the auxiliary function $F_{N_{1},N_{2}-1}$ of the same functional form, with $N_{1}+(N_{2}-1)$ coordinates and corresponding constants $C_{N_{1},N_{2}-1}$, $C^{\prime}_{N_{1},N_{2}-1}$, $D_{N_{1},N_{2}-1}$, and $D^{\prime}_{N_{1},N_{2}-1}$ which depend on all constants appearing in $F_{N_{1},N_{2}}$ and read

[TABLE]

Recall that the initial values of $C_{N_{1},N_{2}}$, $C^{\prime}_{N_{1},N_{2}}$, $D_{N_{1},N_{2}}$, and $D^{\prime}_{N_{1},N_{2}}$ are given by (15) and obtained when we equate the auxiliary function $F_{N_{1},N_{2}}$ and the $N$-body density (III.1).

According to (16), the horizontal reduction makes a ‘stopover’ at the auxiliary function

[TABLE]

Its respective constants, $C_{N_{1},1}$, $C^{\prime}_{N_{1},1}$, $D_{N_{1},1}$, and $D^{\prime}_{N_{1},1}$, are required as initial conditions for the vertical reduction en route to evaluate the inter-species reduced density matrix of the mixture (III.1). Note that the dependence of $F_{N_{1},1}$ on $\beta_{2}$ is now implicit, representing the situation that all but the last boson of type $2$ are integrated out. The final result for the constants is

[TABLE]

where the initial conditions (15) have been used, see appendix B for further details.

We can now perform the vertical reduction of $F_{N_{1},1}$ where (III.1) serve as the initial conditions for the constants to be computed. Thus we seek for a recurrence relation

[TABLE]

where

[TABLE]

The relations between the constants in the vertical reduction (III.1) are seen to be analogous to the relations between the constants in the horizontal reduction (III.1), see appendix B for further details.

Thus, combining and interchanging the order of both horizontal and vertical reductions, the integration of $F_{N_{1},N_{2}}$ (16) ends with the auxiliary functions

[TABLE]

Because $F_{N_{1},N_{2}}$ with the initial conditions (15) is proportional to the $N$-particle density matrix (III.1), $F_{1,1}$ is proportional to the inter-species reduced density matrix (III.1), and similarly $F_{1,0}$ and $F_{0,1}$ to the intra-species ones (III.1). The final expressions for the constants in $F_{1,1}$ are

[TABLE]

and for the constants in $F_{1,0}$ and $F_{0,1}$ are

[TABLE]

The explicit dependence of the bottommost constants (III.1) and (III.1) on $C_{N_{1},1}$ in (III.1) and

[TABLE]

is to remind us that both horizontal and vertical reductions have been combined to evaluate the former, see appendix B for further details. All in all, the inter-species (III.1) and intra-species (III.1) reduced density matrices read

[TABLE]

and

[TABLE]

For completeness, we paste below their diagonal parts, i.e., the two-body

[TABLE]

and one-body

[TABLE]

densities. Summarizing, we have derived closed-form expressions for the (lowest-order) intra-species and inter-species reduced density matrices of a generic model of a trapped mixture of interacting bosons.

III.2 Solution of the harmonic-interaction model for trapped mixtures at the mean-field level

At the other end of the exact, many-body treatment of the harmonic-interaction model for trapped mixtures lies the Gross-Pitaevskii, mean-field solution. In the mean-field theory the many-particle wave-function is approximated as a product state, where all the bosons of species $1$ lie in one orbital $\phi_{1}({\bf x})$ and all the bosons of species $2$ lie in another orbital $\phi_{2}({\bf y})$. Thus, the mean-field ansatz for the mixture is the product wave-function

[TABLE]

The Gross-Pitaevskii energy functional of the mixture reads

[TABLE]

where the mean-field interaction parameters are given by $\Lambda_{1}=\lambda_{1}(N_{1}-1)$, $\Lambda_{2}=\lambda_{1}(N_{2}-1)$, $\Lambda_{12}=\lambda_{12}N_{1}$, and $\Lambda_{21}=\lambda_{12}N_{2}$ and satisfy $N_{1}\Lambda_{21}=N_{2}\Lambda_{12}$. We denote hereafter $\varepsilon^{GP}=\frac{E^{GP}}{N}$ as the total mean-field energy of the mixture divided by the total number of particles $N=N_{1}+N_{2}$. Minimizing the energy functional (III.2) with respect to the shapes of the orbitals $\phi_{1}({\bf x})$ and $\phi_{2}({\bf y})$, the two-coupled Gross-Pitaevskii equations of the mixture are derived,

[TABLE]

where $\mu_{1}$ and $\mu_{2}$ are the chemical potentials of the species, see appendix C.

The solution of (III.2) follows a similar strategy as for the single-species HIM_Cohen and symmetric-mixture BB_HIM_SYM harmonic-interaction models. Expanding the interaction terms in (III.2) we find

[TABLE]

where, since $\phi_{1}({\bf x})$ and $\phi_{2}({\bf y})$ are even functions (see below), there are no linear in ${\bf x}$, ${\bf y}$ terms in (III.2). A particular solution of (III.2) are the following (interaction-dressed) Gaussian functions REMARK_GAUSSIAN

[TABLE]

We can now compute the mean-field energy per particle $\varepsilon^{GP}$ which reads

[TABLE]

where $\frac{N_{1}}{N_{2}}=\frac{\Lambda_{12}}{\Lambda_{21}}$ is used. Of course, the many-body energy (II) is always lower than the mean-field energy because of the variational principle.

We now discuss the reduced density matrices at the mean-field level. The Gross-Pitaevskii wave-function reads

[TABLE]

From (III.2) we have

[TABLE]

for the reduced one-body density matrices and

[TABLE]

for the inter-species reduced density matrix. Quite generally and as might have been expected, for mixtures with a finite number of particles, the mean-field reduced density matrices (III.2) and (III.2) differ from their many-body counterparts (III.1) and (III.1). The intra-species reduced density matrices are factorized to products of Gross-Pitaevskii orbitals, and the inter-species reduced density matrix is factorized to product of Gross-Pitaevskii intra-species reduced density matrices. This concludes our derivation of the mean-field solution of the harmonic-interaction model for trapped mixtures.

III.3 The infinite-particle limit

We are now in the position to put together the above two sections, and investigate the energy and reduced density matrices per particle of the mixture at the infinite-particle limit, and how these quantities are connected in this limit with the Gross-Pitaevskii solution of the mixture. Interestingly, for a mixture we can discuss separately two such limits, hereafter referred to as the two-species infinite-particle limit and the one-species infinite-particle limit. In the first, the numbers of particles of both species are taken to infinity whereas in the second the number of particles of one of the species is taken to infinity and the number of particles of the second species remains fixed and finite (in both limits interaction parameters are held fixed, the precise way is discussed below). This is unlike the case of the single-species and symmetric-mixture harmonic-interaction models. We compare and contrast the properties of the mixture in the two limits. We start with the two-species infinite-particle limit.

III.3.1 The two-species infinite-particle limit

In the two-species infinite-particle limit, namely for $N_{1}\to\infty$ and $N_{2}\to\infty$ and holding the interaction parameters $\Lambda_{1}$, $\Lambda_{2}$, $\Lambda_{12}$, and $\Lambda_{21}$ fixed, we find from (II) that

[TABLE]

which establishes the connection between the exact energy per particle and mean-field (Gross-Pitaevskii) energy per particle in this limit for the generic mixture. Like the literature cases of single-species bosons and the symmetric-mixture harmonic-interaction model, the many-body and mean-field solutions coincide in the limit of an infinite number of particles as far as the energy per particle in examined. Note that the two-species limit of an infinite number of particles $N\to\infty$ implies that $N_{1}\to\infty$ and $N_{2}\to\infty$ such the the ratio $\frac{N_{1}}{N_{2}}$ is kept constant. This is since $\frac{N_{1}}{N_{2}}=\frac{\Lambda_{12}}{\Lambda_{21}}$.

To discuss the two-species infinite-particle limit for the reduced density matrices we first have to evaluate the limit of relevant quantities. Thus we find for the frequencies (II)

[TABLE]

for the parameters (II) of the wave-function

[TABLE]

and therefore for the constants (III.1) and (III.1) of the reduced density matrices

[TABLE]

In particular, that $\lim_{N\to\infty}D_{11}=\lim_{N\to\infty}D^{\prime}_{11}=0$ stems from $\lim_{N\to\infty}\gamma=0$ and implies that there is no coupling at the level of the inter-species reduced density matrix, see (III.1), between the species $1$ and $2$ in the two-species infinite-particle limit. Combining the above we find

[TABLE]

for the reduced one-body density matrices per particle and

[TABLE]

for the intra-species reduced density matrix per particle. With this, we have established the 100% condensation of each species in the generic mixture, in the two-species infinite-particle limit. Furthermore, the inter-species reduced density matrix per particle is separable in this limit and given as a product of the intra-species reduced density matrices per particle. Each condensate is described by the Gross-Pitaevskii quantities. This constitutes a generalization for generic mixtures of interacting bosons, at least within the exactly-solvable harmonic-interaction model for trapped mixtures, of what is known in the literature for single-species trapped Bose-Einstein condensates Yngvason_PRA ; Lieb_PRL .

III.3.2 The one-species infinite-particle limit

Let us discuss what happens in (and how to define) the limit of an infinite number of particles of one of the species, say, species $1$. For $N_{1}\to\infty$ the interaction parameters $\Lambda_{1}=\lambda_{1}(N_{1}-1)$ and $\Lambda_{12}=\lambda_{12}N_{1}$ are held fixed by diminishing the interaction strengths $\lambda_{1}$ and $\lambda_{12}$ accordingly. Since the number of particles of the second species $N_{2}$ is finite (and fixed) the interaction parameter $\Lambda_{2}=\lambda_{2}(N_{2}-1)$ is fixed for constant $\lambda_{2}$. However, $\Lambda_{21}=\lambda_{12}N_{2}\to 0$. Thus, we get for the energy per particle in the one-species infinite-particle limit

[TABLE]

When only species $1$ is taken to the infinite-particle limit it naturally becomes dominant over species $2$. The energy per particle is that of species $1$ only, with apparently no contribution or influence from species $2$, and is given by the Gross-Pitaevskii energy per particle of species $1$ alone. What happens then with the reduced density matrices?

To discuss the one-species infinite-particle limit for the reduced density matrices we first have to evaluate with some care the limit of the relevant quantities. Now we find for the frequencies (II)

[TABLE]

for the parameters (II) of the wave-function

[TABLE]

and consequently for the constants (III.1) and (III.1) of the reduced density matrices

[TABLE]

where the last equality stems from $\lim_{N_{1}\to\infty}C^{\prime}_{1,N_{2}}=\lim_{N_{1}\to\infty}C^{\prime}_{0,N_{2}}=0$ and is instrumental in what follows. In particular, that $\lim_{N_{1}\to\infty}D_{1,1}=\lim_{N_{1}\to\infty}D^{\prime}_{1,1}=0$ stems from $\lim_{N_{1}\to\infty}\gamma=0$ and implies that there is no coupling at the level of the inter-species reduced density matrix, see (III.1), between the species $1$ and $2$ also in the one-species infinite-particle limit. Note that the corresponding quantities associated with species $2$ do not vanish in the one-species infinite-particle limit, compare (III.3.2)-(III.3.2) with (III.3.1)-(III.3.1). Now we can prescribe in the one-species infinite-particle limit the single-particle reduced density matrices per particle,

[TABLE]

and the inter-species reduced density matrix per particle,

[TABLE]

We find that species $1$ is described in this limit by the Gross-Pitaevskii quantity whereas species $2$, as might have been expected, is not, implying that species $2$ remains correlated. Interestingly, the inter-species reduced density matrix per particle is separable in this limit, and precisely given by the product of the Gross-Pitaevskii quantity for species $1$ and the correlated quantity for species $2$, see (III.3.2). This is on the account of the last equality in (III.3.2). This concludes our studies of the reduced density matrices of a generic mixture within the harmonic-interaction model for trapped mixtures in the single-species infinite-particle limit.

IV Concluding Remarks

Are the different species in the ground state of a trapped bosonic mixture $100\%$ condensed? Are the many-body and mean-field energies per particle of a trapped bosonic mixture equal at the infinite-particle limit? In the present work we answered these and more questions by treating an exactly-solvable model – the harmonic-interaction model for trapped bosonic mixtures.

From the ground-state wave-function of the mixture we have computed the lowest-order intra-species and inter-species reduced density matrices, by generalizing Cohen and Lee HIM_Cohen recurrence relations for the single-species harmonic-interaction model to a generic mixture. We have also obtained analytically the Gross-Pitaevskii solution for the ground state of the mixture. Thereafter, by taking the infinite-particle limit with respect to the two species, we were able to show that each of the species is indeed $100\%$ condensed, and that the many-body and Gross-Pitaevskii quantities for the energy per particle and reduced density matrices per particle coincide in this limit. When the infinite-particle limit is taken with respect to one of the species, only this species becomes $100\%$ condensed, whereas the other species remains correlated. Interestingly, in either of the infinite-particle-limit procedures the intra-species density matrix per particle becomes exactly the product of the intra-species reduced density matrices per particle.

It would be interesting to investigate other properties of the harmonic-interaction model for mixtures presented in the present work, e.g., quantities whose many-body and mean-field descriptions may not coincide in the infinite-particle limit Variance ; WAVE_OVLP . For finite systems, the model may prove deductive as well, for instance to investigate properties of an impurity made of a few interacting particles embedded inside a larger Bose-Einstein condensate, and to benchmark numerical tools.

Acknowledgements

I thank Shachar Klaiman, Alexej Streltsov, and Lorenz Cederbaum for fruitful discussions. This research was supported by the Israel Science Foundation (Grant No. 600/15).

Appendix A Further Details of Diagonalizing the Hamiltonain of the Mixture

Opening the braces of the particle-particle interaction terms in (II) and collecting the diagonal contributions together we have

[TABLE]

With the Jacobi-coordinate transformation (II) the harmonic trapping and kinetic energy remain diagonal, since from

[TABLE]

and

[TABLE]

we have

[TABLE]

Finally, using the quadratic relations

[TABLE]

the coupling terms in the Hamiltonian (52) are diagonalized too, and the ground-state wave-function transforms from the Jacobi-coordinate to the lab-frame representation, see (II) and (II).

Appendix B Further Details of the Solution of the Coupled Recurrence Relations

Within the horizontal reduction of $F_{N_{1},N_{2}}$, see (16)-(III.1), the recurrence relations between the corresponding constants are

[TABLE]

The recurrence relations (B) consist of one non-linear relation (for $C^{\prime}_{N_{1},j}$) and the linear relations for $D_{N_{1},j}$, $C_{N_{1},j}$, and $D^{\prime}_{N_{1},j}$. We see from (B) that first the relation for $C^{\prime}_{N_{1},j}$ needs to be solved, then for the $D_{N_{1},j}$, and then for the $C_{N_{1},j}$. The linear relation for $D^{\prime}_{N_{1},j}$ is trivial and does not depend on the other constants.

The recurrence relation for $C^{\prime}_{N_{1},j}$ has exactly the same structure as the recurrence relation emerging in the integration of the single-species harmonic-interaction model HIM_Cohen . Making use of this observation and substituting the result into $D_{N_{1},j}$ and then together into $C_{N_{1},j}$, relations (B) are solved. The final result reads

[TABLE]

where the initial conditions for $C^{\prime}_{N_{1},N_{2}}$, $D_{N_{1},N_{2}}$, $C_{N_{1},N_{2}}$, and $D^{\prime}_{N_{1},N_{2}}$ in (15) have been used. From (B) we obtain the constants $C^{\prime}_{N_{1},1}$, $D_{N_{1},1}$, $C_{N_{1},1}$, and $D^{\prime}_{N_{1},1}$ in (III.1) entering the auxiliary function $F_{N_{1},1}$, and $C_{N_{1},0}$ in (III.1) entering the auxiliary function $F_{N_{1},0}$.

Proceeding now to the vertical reduction of $F_{N_{1},1}$ to $F_{1,1}$, see (16), (III.1), and (III.1), we write

[TABLE]

The vertical recursion relations (B) consist again one non-linear relation (for $C_{j,1}$), and the linear relations for $C^{\prime}_{j,1}$, $D_{j,1}$, and $D^{\prime}_{j,1}$. The linear relation for $D^{\prime}_{j,1}$ is again trivial. Note the interchange of roles of the $C$ and $C^{\prime}$ constants when moving from horizontal to vertical reductions. We can hence interchange the order of horizontal and vertical reductions in order to prescribe $C^{\prime}_{1,1}$ based on the solution for $C_{1,1}$. Then, $C^{\prime}_{1,N_{2}}$ in (III.1) obtained from the vertical reduction is used as an initial condition for the recurrence relation. Combining all the above we find for the constants of the auxiliary function $F_{1,1}$

[TABLE]

Upon substitution we arrive at the final result (III.1).

Last are the constants of $F_{1,0}$ and $F_{0,1}$. From the recurrence relation

[TABLE]

we find

[TABLE]

where the initial conditions $C_{N_{1},0}$ and $C^{\prime}_{0,N_{2}}$ are given in (III.1), and where interchanging the order of the vertical and horizontal reductions allows one to obtain $C^{\prime}_{0,1}$ analogously to $C_{1,0}$. Substituting all quantities we arrive at the final expressions for the constants (III.1).

Appendix C Further Details of the Mean-Field Solution of the Mixture

Given the orbitals $\phi_{1}({\bf x})$ and $\phi_{2}({\bf y})$ (III.2), we can now evaluate the integrals $\int d{\bf x}^{\prime}|\phi_{1}({\bf x}^{\prime})|^{2}{\bf x}^{\prime 2}=\frac{3}{2\sqrt{\omega^{2}+\frac{2}{m_{1}}(\Lambda_{1}+\Lambda_{21})}}$ and $\int d{\bf y}^{\prime}|\phi_{2}({\bf y}^{\prime})|^{2}{\bf y}^{\prime 2}=\frac{3}{2\sqrt{\omega^{2}+\frac{2}{m_{1}}(\Lambda_{2}+\Lambda_{12})}}$ in (III.2) and determine the chemical potentials

[TABLE]

With this, the mean-field energy is given by

[TABLE]

where $\int d{\bf x}d{\bf x}^{\prime}|\phi_{1}({\bf x})|^{2}|\phi_{1}({\bf x}^{\prime})|^{2}({\bf x}-{\bf x}^{\prime})^{2}=\frac{3}{\sqrt{\omega^{2}+\frac{2}{m_{1}}(\Lambda_{1}+\Lambda_{21})}}$, $\int d{\bf y}d{\bf y}^{\prime}|\phi_{2}({\bf y})|^{2}|\phi_{2}({\bf y}^{\prime})|^{2}({\bf y}-{\bf y}^{\prime})^{2}=\frac{3}{\sqrt{\omega^{2}+\frac{2}{m_{2}}(\Lambda_{2}+\Lambda_{12})}}$, and $\int d{\bf x}d{\bf y}|\phi_{1}({\bf x})|^{2}|\phi_{2}({\bf y})|^{2}({\bf x}-{\bf y})^{2}=\frac{3}{2}\left[\frac{1}{\sqrt{\omega^{2}+\frac{2}{m_{1}}(\Lambda_{1}+\Lambda_{21})}}+\frac{1}{\sqrt{\omega^{2}+\frac{2}{m_{2}}(\Lambda_{2}+\Lambda_{12})}}\right]$ are used.