Fragmentation of Identical and Distinguishable Bosons’ Pairs and Natural Geminals of a Trapped Bosonic Mixture

Abstract

In a mixture of two kinds of identical bosons, there are two types of pairs: identical bosons’ pairs, of either species, and pairs of distinguishable bosons. In the present work, the fragmentation of pairs in a trapped mixture of Bose–Einstein condensates is investigated using a solvable model, the symmetric harmonic-interaction model for mixtures. The natural geminals for pairs made of identical or distinguishable bosons are explicitly contracted by diagonalizing the intra-species and inter-species reduced two-particle density matrices, respectively. Properties of pairs’ fragmentation in the mixture are discussed, the role of the mixture’s center-of-mass and relative center-of-mass coordinates is elucidated, and a generalization to higher-order reduced density matrices is made. As a complementary result, the exact Schmidt decomposition of the wave function of the bosonic mixture is constructed. The entanglement between the two species is governed by the coupling of their individual center-of-mass coordinates, and it does not vanish at the limit of an infinite number of particles where any finite-order intra-species and inter-species reduced density matrix per particle is 100% condensed. Implications are briefly discussed.

1. Introduction

Condensation and fragmentation are basic and widely-studied concepts of Bose–Einstein condensates emanating from the properties of the reduced one-particle density matrix [1,2,3,4,5]. The bosons are said to be condensed if there is a single macroscopic eigenvalue of the reduced one-particle density matrix [6] and fragmented if there are two or more such macroscopic eigenvalues [7]. These eigenvalues are commonly called natural occupation numbers, and the respective eigenfunctions of the reduced one-particle density matrix are referred to as natural orbitals. The fragmentation of Bose–Einstein condensates has been investigated, e.g., in [8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27].

The condensation and especially fragmentation of the reduced two-particle density matrix of interacting identical bosons is less studied; see, e.g., [28]. Here, the analysis of the reduced two-particle density matrix would determine whether pairs of bosons are condensed or fragmented. The respective eigenfunctions of the reduced two-particle density matrix are often called natural geminals. We note that natural geminals in electronic systems have long been explored; see, e.g., [29,30,31,32,33,34,35,36,37,38,39,40].

Consider now a mixture of two kinds of identical bosons, which are labeled species 1 and species 2. Mixtures of Bose–Einstein condensates is a highly investigated topic; see, e.g., [41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70]. One may ask, just like for the single-species bosons, about the condensation or fragmentation of each of the species and how, for instance, one species is affected by the presence of the other species and vice versa. To answer this question, the intra-species reduced one-particle density matrices of species 1 and 2 are required, i.e., analyzing the intra-species occupation numbers and natural orbitals. Following the above line, one could also investigate the fragmentation of higher-order intra-species reduced density matrices in the mixture. For instance, to investigate whether pairs of identical bosons, of either species 1 or species 2, are fragmented, diagonalizing the intra-species reduced two-particle density matrices is needed. In summary, the fragmentation of identical bosons and its manifestation in higher-order reduced density matrices stem from the properties of intra-species quantities.

However, a mixture of Bose–Einstein condensates offers a degree-of-freedom or many-particle construction which does not exist for single-species bosons, namely, inter-species reduced density matrices. Now, if the fragmentation of identical bosons and pairs is defined as the macroscopic occupation of respective eigenvalues following the diagonalization of intra-species reduced density matrices, we may analogously define the fragmentation of distinguishable bosons’ pairs as a macroscopic occupation of the eigenvalues of the inter-species reduced two-particle density matrix. Obviously, the latter is the lowest-order inter-species quantity, since at least one particle of each species is needed to build an inter-species entity.

The above discussion defines the goals of the present work, which are: (i) to investigate the fragmentation of pairs of identical bosons and establish the fragmentation of pairs of distinguishable bosons in a mixture of Bose–Einstein condensates; (ii) to construct the respective natural geminals of the mixture, for identical pairs and for distinguishable pairs; (iii) to show that the fragmentation of distinguishable bosons’ pairs in the mixture persists with higher-order inter-species reduced density matrices; (iv) to construct the Schmidt decomposition of the mixture’s wave function and discuss some of its properties at the limit of an infinite-number of particles where the mixture is $100\%$ condensed; and (v) to achieve the first four goals analytically, using an exactly solvable model.

To this end, we recruited the harmonic-interaction model for mixtures [71,72,73,74,75,76], or, more precisely here, a symmetric version of which [77]. The harmonic-interaction model for single-species bosons (and fermions) has been used extensively in the literature, including for investigating the properties of Bose–Einstein condensates [78,79,80,81,82,83,84,85,86,87,88,89,90,91,92]. In our work, we built on results obtained and techniques used for the reduced density matrices of single-species bosons within the harmonic-interaction model [78,79,80,81,82], and, among others, generalized and extended them for the intra-species and particularly the inter-species reduced density matrices of mixtures [74].

The structure of the paper is as follows: In Section 2, we construct and investigate the fragmentation of intra-species and inter-species pair functions in the mixture. In Section 3, we extend the results and explore the fragmentation of pairs of distinguishable pairs. Furthermore, a complementary result for the Schmidt decomposition of the mixture’s wavefucntion at the limit of an infinite number of particles is offered. In Section 4, a summary of the results and an outlook of some prospected research topics are provided. Finally, Appendix A collects the details of fragmentation of bosons and pairs in the single-species system for a comparison with the mixture.

2. Intra-Species and Inter-Species Natural Pair Functions

2.1. The Symmetric Two-Species Harmonic-Interaction Model

We consider a mixture of two Bose–Einstein condensates described by the Hamiltonian of the symmetric two-species harmonic-interaction model [74,77]:

$$\begin{array}{cc}\hfill & \widehat{H}(x_1,\dots ,x_N,y_1,\dots ,x_N)={\sum }_{j=1}^N\left(-\frac{1}{2m}\frac{{\partial }^2}{\partial x_j^2}+\frac{1}{2}m{\omega }^2{x}_j^2\right)+\lambda {\sum }_{1\le j<k}^N{\left(x_j-x_k\right)}^2+\hfill \\ \hfill & +{\sum }_{j=1}^N\left(-\frac{1}{2m}\frac{{\partial }^2}{\partial y_j^2}+\frac{1}{2}m{\omega }^2{y}_j^2\right)+\lambda {\sum }_{1\le j<k}^N{\left(y_j-y_k\right)}^2+{\lambda }_{12}{\sum }_{j=1}^N{\sum }_{k=1}^N{\left(x_j-y_k\right)}^2.\hfill \end{array}$$

(1)

There are N bosons of type 1, N bosons of type 2, and the mass of each boson is m. $\lambda$ is the intra-species interaction strength, either between two bosons of type 1 or two bosons of type 2, and ${\lambda }_{12}$ is the inter-species interaction strength between type 1 and type 2 bosons. Dimensionality plays no role in the present work; hence, we work in one spatial dimension. $\hslash =1$ is used throughout. Employing Jacoby coordinates for the mixture and translating back to the laboratory frame, the $2N$-boson wave function and corresponding many-particle density matrix are given by

$$\begin{array}{cc}\hfill & \Psi (x_1,\dots ,x_N,y_1,\dots ,y_N)={\left(\frac{m\Omega }{\pi }\right)}^{\frac{N-1}{2}}{\left(\frac{M_{12}{\Omega }_{12}}{\pi }\right)}^{\frac{1}{4}}{\left(\frac{M\omega }{\pi }\right)}^{\frac{1}{4}}\hfill \\ \hfill & \times e^{-\frac{\alpha }{2}{\sum }_{j=1}^N\left(x_j^2+y_j^2\right)-\beta {\sum }_{1\le j<k}^N\left(x_j{x}_k+y_j{y}_k\right)+\gamma {\sum }_{j=1}^N{\sum }_{k=1}^N{x}_j{y}_k},\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill & \Psi (x_1,\dots ,x_N,y_1,\dots ,y_N){\Psi }^{*}(x_1^{\prime },\dots ,x_N^{\prime },y_1^{\prime },\dots ,y_N^{\prime })={\left(\frac{m\Omega }{\pi }\right)}^{N-1}{\left(\frac{M_{12}{\Omega }_{12}}{\pi }\right)}^{\frac{1}{2}}{\left(\frac{M\omega }{\pi }\right)}^{\frac{1}{2}}\hfill \\ \hfill & \times e^{-\frac{\alpha }{2}{\sum }_{j=1}^N\left(x_j^2+{x_j^{\prime }}^2+y_j^2+{y_j^{\prime }}^2\right)-\beta {\sum }_{1\le j<k}^N\left(x_j{x}_k+x_j^{\prime }{x}_k^{\prime }+y_j{y}_k+y_j^{\prime }{y}_k^{\prime }\right)+\gamma {\sum }_{j=1}^N{\sum }_{k=1}^N\left(x_j{y}_k+x_j^{\prime }{y}_k^{\prime }\right)},\hfill \end{array}$$

(2b)

with

$$\begin{array}{cc}\hfill & \Omega =\sqrt{{\omega }^2+\frac{2N}{m}\left(\lambda +{\lambda }_{12}\right)},{\Omega }_{12}=\sqrt{{\omega }^2+\frac{4N}{m}{\lambda }_{12}},\hfill \\ \hfill & \alpha =m\Omega +\beta ,\beta =\frac{m}{2N}\left({\Omega }_{12}+\omega -2\Omega \right),\gamma =\frac{m}{2N}\left({\Omega }_{12}-\omega \right),\hfill \end{array}$$

(2c)

and the relative center-of-mass $M_{12}=\frac{m}{2N}$ and center-of-mass $M=2mN$ masses. The wave function and similarly the many-particle density of the mixture depend on two dressed frequencies, $\Omega$ and ${\Omega }_{12}$, and consist of three parts: One-body part with coefficient $\alpha$, intra-species two-body coupling with coefficient $\beta$, and inter-species two-body coupling with coefficient $\gamma$, whereas $\alpha$ and $\beta$ depend on the intra-species and inter-species interactions, $\gamma$ depends on the inter-species interaction only.

Another issue worth mentioning is the stability region of the mixture, namely when the mixture is bound. The condition that a bound solution of the mixture exists requires that all $2N$ degrees-of-freedom (oscillators after diagonalization with the Jacoby coordinates) are bound. Namely, the $2(N-1)$ relative-coordinate oscillators, which depend on both the intra-species and the inter-species interactions $\lambda$ and ${\lambda }_{12}$, should be bound; the relative center-of-mass oscillator, which depends on the inter-species interaction ${\lambda }_{12}$ only, should be bound; and the center-of-mass oscillator which does not depend on any of the two interactions and is always bound. Thus, for a bound mixture, both dressed frequencies $\Omega$ and ${\Omega }_{12}$ must be positive. This implies the conditions $\lambda +{\lambda }_{12}>-\frac{m{\omega }^2}{2N}$ and ${\lambda }_{12}>-\frac{m{\omega }^2}{4N}$, respectively, on the interactions. In other words, the inter-species interaction ${\lambda }_{12}$ is bound from below, implying that the mutual repulsion between the two species cannot be too strong but is not bound from above, meaning that the mutual attraction between the two species can be unlimitedly strong. Furthermore, the intra-species interaction $\lambda$ can take any value as long as the inter-species interaction is sufficiently attractive, i.e., $\lambda >-\frac{m{\omega }^2}{2N}-{\lambda }_{12}$. We shall return to the dressed frequencies $\Omega$ and ${\Omega }_{12}$ below. Note that in our work, we use the notions bound and stable interchangeably. For example, the border of stability of the mixture is the above-described border of parameters within which the ground state exists.

2.2. Intra-Species Natural Pair Functions

The intra-species reduced density matrices are defined when all bosons of the other type are integrated out. We concentrate in what follows on the reduced one-particle and in particular the two-particle density matrices of species 1,

$$\begin{array}{cc}\hfill & {\rho }_1^{\left(1\right)}(x,x^{\prime })=N\int d{x}_2\cdots d{x}_{N}d{y}_1\cdots d{y}_N\Psi (x,x_2,\dots ,x_N,y_1,\dots ,y_N)\hfill \\ \hfill & \times {\Psi }^{*}(x^{\prime },x_2,\dots ,x_N,y_1,\dots ,y_N),\hfill \\ \hfill & {\rho }_1^{\left(2\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime })=N(N-1)\int d{x}_3\cdots d{x}_{N}d{y}_1\cdots d{y}_N\Psi (x_1,x_2,x_3,\dots ,x_N,y_1,\dots ,y_N)\hfill \\ \hfill & \times {\Psi }^{*}(x_1^{\prime },x_2^{\prime },x_3,\dots ,x_N,y_1,\dots ,y_N).\hfill \end{array}$$

(3)

In a symmetric mixture, the corresponding reduced density matrices of species 2, ${\rho }_2^{\left(1\right)}(y,y^{\prime })$ and ${\rho }_2^{\left(2\right)}(y_1,y_2,y_1^{\prime },y_2^{\prime })$, are the same and need not be repeated.

The reduction of the many-particle density (2b) to its finite-order reduced density matrices is somewhat lengthy and given in [74]. We start from the final expression for the intra-species reduced one-particle density matrix, which is given by

$$\begin{array}{cc}\hfill & {\rho }_1^{\left(1\right)}(x,x^{\prime })=N{\left(\frac{\alpha +C_{1,0}}{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{\alpha }{2}\left(x^2+{x^{\prime }}^2\right)}{e}^{-\frac{1}{4}{C}_{1,0}{\left(x+x^{\prime }\right)}^2}\hfill \\ \hfill & =N{\left(\frac{\alpha +C_{1,0}}{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{\alpha +\frac{C_{1,0}}{2}}{2}\left(x^2+{x^{\prime }}^2\right)}{e}^{-\frac{1}{2}{C}_{1,0}x{x}^{\prime }},\hfill \\ \hfill & \alpha +C_{1,0}=(\alpha -\beta )\frac{{\left[(\alpha -\beta )+N\beta \right]}^2-N^2{\gamma }^2}{\left[(\alpha -\beta )+N\beta \right]\left[(\alpha -\beta )+(N-1)\beta \right]-N(N-1){\gamma }^2}.\hfill \end{array}$$

(4)

The coefficient $C_{1,0}$ governs the properties of the intra-species reduced one-particle density matrix and reminds one that all bosons of type 2 and all but a single boson of type 1 are integrated out. As might be expected, ${\rho }_1^{\left(1\right)}(x,x^{\prime })$ depends on the three parts of the many-boson wave function, i.e., on the $\alpha$, $\beta$, and $\gamma$ terms (2a). In the absence of inter-species interaction, i.e., for $\gamma =0$, the coefficient $C_{1,0}$ boils down to that of the single-species harmonic-interaction model; see Appendix A for further discussion.

Just as for the case of single-species bosons [81,82], the intra-species reduced one-particle density matrix (4) can be diagonalized using Mehler’s formula. Mehler’s formula can be written as follows:

$$\begin{array}{cc}\hfill & {\left[\frac{(1-\rho )s}{(1+\rho )\pi }\right]}^{\frac{1}{2}}{e}^{-\frac{1}{2}\frac{(1+{\rho }^2)s}{1-{\rho }^2}\left(x^2+{x^{\prime }}^2\right)}{e}^{+\frac{2\rho s}{1-{\rho }^2}x{x}^{\prime }}\hfill \\ \hfill & ={\sum }_{n=0}^{\infty }(1-\rho ){\rho }^n\frac{1}{\sqrt{2^{n}n!}}{\left(\frac{s}{\pi }\right)}^{\frac{1}{4}}{H}_n\left(\sqrt{s}x\right)e^{-\frac{1}{2}s{x}^2}\frac{1}{\sqrt{2^{n}n!}}{\left(\frac{s}{\pi }\right)}^{\frac{1}{4}}{H}_n\left(\sqrt{s}{x}^{\prime }\right)e^{-\frac{1}{2}s{x^{\prime }}^2},\hfill \end{array}$$

(5)

with $s>0$ and, generally, for intra-species and inter-species reduced density matrices as well as later on for Schmidt decomposition of the wave function, $1>\rho \ge 0$. $H_n$ are the Hermite polynomials.

Comparing the structure of the intra-species reduced one-particle density matrix ${\rho }_1^{\left(1\right)}(x,x^{\prime })$ with that of Mehler’s formula one readily has

$$\begin{array}{cc}\hfill & s_1^{\left(1\right)}=\sqrt{\alpha \left(\alpha +C_{1,0}\right)}=\sqrt{\frac{\alpha (\alpha -\beta )\left\{{\left[(\alpha -\beta )+N\beta \right]}^2-N^2{\gamma }^2\right\}}{\left[(\alpha -\beta )+N\beta \right]\left[(\alpha -\beta )+(N-1)\beta \right]-N(N-1){\gamma }^2}},\hfill \\ \hfill & {\rho }_1^{\left(1\right)}=\frac{\alpha -s_1^{\left(1\right)}}{\alpha +s_1^{\left(1\right)}}=\frac{\sqrt{\frac{\alpha \left\{\left[(\alpha -\beta )+N\beta \right]\left[(\alpha -\beta )+(N-1)\beta \right]-N(N-1){\gamma }^2\right\}}{(\alpha -\beta )\left\{{\left[(\alpha -\beta )+N\beta \right]}^2-N^2{\gamma }^2\right\}}}-1}{\sqrt{\frac{\alpha \left\{\left[(\alpha -\beta )+N\beta \right]\left[(\alpha -\beta )+(N-1)\beta \right]-N(N-1){\gamma }^2\right\}}{(\alpha -\beta )\left\{{\left[(\alpha -\beta )+N\beta \right]}^2-N^2{\gamma }^2\right\}}}+1},\hfill \\ \hfill & 1-{\rho }_1^{\left(1\right)}=\frac{2s_1^{\left(1\right)}}{\alpha +s_1^{\left(1\right)}},\hfill \end{array}$$

(6)

where $1-{\rho }_1^{\left(1\right)}$ is the condensate fraction of species 1 (and of species 2), i.e., the fraction of condensed bosons, and ${\rho }_1^{\left(1\right)}$ is the depleted fraction, namely, the fraction of bosons residing outside the lowest, condensed mode. $s_1^{\left(1\right)}$ is the scaling, or effective frequency, of the intra-species natural orbitals. The condensate fraction, depleted fraction, and scaling of the natural orbitals are all given in closed form as a function of the number of bosons N, and the intra-species $\lambda$ and inter-species ${\lambda }_{12}$ interaction strengths. A specific application of the general expressions (6) for the mixture appears below.

For the intra-species two-particle reduced density matrix, we have:

$$\begin{array}{cc}\hfill & {\rho }_1^{\left(2\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime })=N(N-1){\left(\frac{\alpha +C_{1,0}}{\pi }\right)}^{\frac{1}{2}}{\left(\frac{\alpha +C_{2,0}}{\pi }\right)}^{\frac{1}{2}}\hfill \\ \hfill & \times e^{-\frac{\alpha }{2}\left(x_1^2+x_2^2+{x_1^{\prime }}^2+{x_2^{\prime }}^2\right)}{e}^{-\beta \left(x_1{x}_2+x_1^{\prime }{x}_2^{\prime }\right)}{e}^{-\frac{1}{4}{C}_{2,0}{\left(x_1+x_2+x_1^{\prime }+x_2^{\prime }\right)}^2},\hfill \\ \hfill & \alpha +\beta +2C_{2,0}=(\alpha -\beta )\frac{{\left[(\alpha -\beta )+N\beta \right]}^2-N^2{\gamma }^2}{\left[(\alpha -\beta )+N\beta \right]\left[(\alpha -\beta )+(N-2)\beta \right]-N(N-2){\gamma }^2},\hfill \end{array}$$

(7)

where $C_{1,0}$ is the coefficient of the intra-species reduced one-particle density matrix (4), and the combination $\left(\alpha +\beta +2C_{2,0}\right)$ would appear shortly after.

To obtain the natural geminals of ${\rho }_1^{\left(2\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime })$, we define the variables $q_1=\frac{1}{\sqrt{2}}\left(x_1+x_2\right)$, $q_2=\frac{1}{\sqrt{2}}\left(x_1-x_2\right)$ and $q_1^{\prime }=\frac{1}{\sqrt{2}}\left(x_1^{\prime }+x_2^{\prime }\right)$, $q_2^{\prime }=\frac{1}{\sqrt{2}}\left(x_1^{\prime }-x_2^{\prime }\right)$, i.e., the center-of-mass and relative coordinate of two identical bosons. With this rotation of coordinates, we have for the different terms in (7):

$$\begin{array}{cc}\hfill & x_1^2+x_2^2+{x_1^{\prime }}^2+{x_2^{\prime }}^2=q_1^2+{q_1^{\prime }}^2+q_2^2+{q_2^{\prime }}^2,\hfill \\ \hfill & x_1{x}_2+x_1^{\prime }{x}_2^{\prime }=\frac{1}{2}\left(q_1^2+{q_1^{\prime }}^2-q_2^2-{q_2^{\prime }}^2\right),\hfill \\ \hfill & {\left(x_1+x_2+x_1^{\prime }+x_2^{\prime }\right)}^2=2\left(q_1^2+{q_1^{\prime }}^2+2q_1{q}_1^{\prime }\right).\hfill \end{array}$$

(8)

Consequently, one readily finds the diagonal form

$$\begin{array}{cc}\hfill & {\rho }_1^{\left(2\right)}(q_1,q_1^{\prime },q_2,q_2^{\prime })=N(N-1){\left(\frac{\alpha -\beta }{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{\alpha -\beta }{2}\left(q_2^2+{q_2^{\prime }}^2\right)}\hfill \\ \hfill & \times {\left(\frac{\alpha +\beta +2C_{2,0}}{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{\alpha +\beta +C_{2,0}}{2}\left(q_1^2+{q_1^{\prime }}^2\right)}{e}^{-C_{2,0}{q}_1{q}_1^{\prime }},\hfill \end{array}$$

(9)

where the normalization coefficients before and after diagonalization are, of course, equal and satisfy $\left(\alpha +C_{1,0}\right)\left(\alpha +C_{2,0}\right)=\left(\alpha -\beta \right)\left(\alpha +\beta +2C_{2,0}\right)$.

After the transformation (8), the first term of ${\rho }_1^{\left(2\right)}(q_1,q_1^{\prime },q_2,q_2^{\prime })$ is separable as a function of $q_2$ and $q_2^{\prime }$, whereas by using Mehler’s formula onto the variables $q_1$ and $q_1^{\prime }$, the second term can be diagonalized. Thus, comparing the second term in (9) and Equation (5), we find

$$\begin{array}{cc}\hfill & s_1^{\left(2\right)}=\sqrt{\left(\alpha +\beta \right)\left(\alpha +\beta +2C_{2,0}\right)}=\sqrt{\frac{({\alpha }^2-{\beta }^2)\left\{{\left[(\alpha -\beta )+N\beta \right]}^2-N^2{\gamma }^2\right\}}{\left[(\alpha -\beta )+N\beta \right]\left[(\alpha -\beta )+(N-2)\beta \right]-N(N-2){\gamma }^2}},\hfill \\ \hfill & {\rho }_1^{\left(2\right)}=\frac{\left(\alpha +\beta \right)-s_1^{\left(2\right)}}{\left(\alpha +\beta \right)+s_1^{\left(2\right)}}=\frac{\sqrt{\frac{\left(\alpha +\beta \right)\left\{\left[(\alpha -\beta )+N\beta \right]\left[(\alpha -\beta )+(N-2)\beta \right]-N(N-2){\gamma }^2\right\}}{(\alpha -\beta )\left\{{\left[(\alpha -\beta )+N\beta \right]}^2-N^2{\gamma }^2\right\}}}-1}{\sqrt{\frac{\left(\alpha +\beta \right)\left\{\left[(\alpha -\beta )+N\beta \right]\left[(\alpha -\beta )+(N-2)\beta \right]-N(N-2){\gamma }^2\right\}}{(\alpha -\beta )\left\{{\left[(\alpha -\beta )+N\beta \right]}^2-N^2{\gamma }^2\right\}}}+1},\hfill \\ \hfill & 1-{\rho }_1^{\left(2\right)}=\frac{2s_1^{\left(2\right)}}{\alpha +s_1^{\left(2\right)}}.\hfill \end{array}$$

(10)

With expressions (10), the decomposition of the intra-species reduced two-particle density matrix in terms of its natural geminals is explicitly given by

$$\begin{array}{cc}\hfill & {\rho }_1^{\left(2\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime })=N(N-1){\sum }_{n=0}^{\infty }\left(1-{\rho }_1^{\left(2\right)}\right){\left({\rho }_1^{\left(2\right)}\right)}^n{\Phi }_{1,n}^{\left(2\right)}(x_1,x_2){\Phi }_{1,n}^{\left(2\right),*}(x_1^{\prime },x_2^{\prime }),\hfill \\ \hfill & {\Phi }_{1,n}^{\left(2\right)}(x_1,x_2)=\frac{1}{\sqrt{2^{n}n!}}{\left(\frac{s_1^{\left(2\right)}}{\pi }\right)}^{\frac{1}{4}}{H}_n\left(\sqrt{\frac{s_1^{\left(2\right)}}{2}}\left(x_1+x_2\right)\right)e^{-\frac{1}{4}{s}_1^{\left(2\right)}{\left(x_1+x_2\right)}^2}\hfill \\ \hfill & \times {\left(\frac{\alpha -\beta }{\pi }\right)}^{\frac{1}{4}}{e}^{-\frac{1}{4}\left(\alpha -\beta \right){\left(x_1-x_2\right)}^2}.\hfill \end{array}$$

(11)

Equation (11) is a general result on the intra-species natural geminals of the mixture. Together with (10), they imply that $1-{\rho }_2^{\left(1\right)}$ is the fraction of condensed pairs of species 1 (and of species 2), ${\rho }_1^{\left(2\right)}$ is the fraction of depleted pairs, i.e., the fraction of pairs residing outside the lowest, condensed natural geminal, and $s_1^{\left(2\right)}$ is the scaling, or effective frequency, of the intra-species natural pair functions. The intra-species natural geminals along with their condensate and depleted fractions are prescribed as explicit functions of the number of bosons N, and the intra-species $\lambda$ and inter-species ${\lambda }_{12}$ interactions. A specific application of the general decomposition (10) and (11) to natural geminals of the mixture is provided below. Finally, we point out that the generalization to higher-order intra-species reduced density matrices and corresponding natural functions follow the above pattern and are not discussed further here.

Let us work out an explicit application where we shall find and analyze fragmentation of identical boson’s pairs. Consider the specific scenario where $\lambda +{\lambda }_{12}=0$, i.e., that the intra-species interaction is inverse to and ‘compensates’ the effect of the inter-species interaction on each of the species in the manner that the intra-species frequency is that of non-interacting particles, $\Omega =\omega$. This implies that the frequency of $2N-1$ oscillators, the $2(N-1)$ relative-coordinate oscillators, and the center-of-mass oscillator is $\omega$, and that only the relative center-of-mass oscillator carries an interaction-dressed frequency which is ${\Omega }_{12}$. We would see the consequences below. Then, the coefficients of the three parts of the wave function simplify, and one has $\alpha =m\omega +\beta =m\omega \left[1+\frac{1}{2N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]$ and $\beta =\gamma =\frac{m}{2N}\left({\Omega }_{12}-\omega \right)$. Consequently, the expressions (6) and (10) simplify, and the intra-species reduced one-particle and two-particle density matrices can be evaluated further. Thus, we readily find

$$\begin{array}{cc}\hfill & s_1^{\left(1\right)}=m\omega \sqrt{\frac{1+\frac{1}{2N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)}{1+\frac{1}{2N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}},\hfill \\ \hfill & {\rho }_1^{\left(1\right)}=\frac{\sqrt{\left[1+\frac{1}{2N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{1}{2N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}-1}{\sqrt{\left[1+\frac{1}{2N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{1}{2N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}+1},\hfill \\ \hfill & 1-{\rho }_1^{\left(1\right)}=\frac{2}{\sqrt{\left[1+\frac{1}{2N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{1}{2N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}+1}\hfill \end{array}$$

(12a)

for the intra-species reduced one-particle density matrix, where $\alpha +C_{1,0}=m\omega \frac{1}{1+\frac{1}{2N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}$ is used, and

$$\begin{array}{cc}\hfill & s_1^{\left(2\right)}=m\omega \sqrt{\frac{1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)}{1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}},\hfill \\ \hfill & {\rho }_1^{\left(2\right)}=\frac{\sqrt{\left[1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}-1}{\sqrt{\left[1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}+1},\hfill \\ \hfill & 1-{\rho }_1^{\left(2\right)}=\frac{2}{\sqrt{\left[1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}+1}\hfill \end{array}$$

(12b)

for the intra-species reduced two-particle density matrix, where $\alpha +\beta =m\omega \left[1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]$ and $\alpha +\beta +2C_{2,0}=m\omega \frac{1}{1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}$ are utilized. We see that the fragmentation of identical pairs and bosons are governed by the ratio $\frac{{\Omega }_{12}}{\omega }$ and its inverse $\frac{\omega }{{\Omega }_{12}}$, meaning that it takes place both at the attractive and repulsive sectors of interactions. Moreover, the condensed and depleted fractions of the pairs and bosons are symmetric to interchanging $\frac{{\Omega }_{12}}{\omega }$ and $\frac{\omega }{{\Omega }_{12}}$; see discussion below.

Let us analyze explicitly macroscopic fragmentation of geminals, i.e., when there is a macroscopic occupation of more than a single intra-species natural pair function of ${\rho }_1^{\left(2\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime })$. As a reference, we also refer to the corresponding and standardly defined macroscopic fragmentation of the intra-species natural orbitals of ${\rho }_1^{\left(1\right)}(x,x^{\prime })$. The structure of the eigenvalues, emanating from Mehler’s formula and its applicability to the various reduced density matrices, suggests that, say, the ‘middle’ value $\rho =1-\rho =\frac{1}{2}$, i.e., when the condensed and depleted fractions are equal, is a convenient manifestation of macroscopic fragmentation. Indeed, for this value the first few natural occupation fractions $(1-\rho ){\rho }^n$, $n=0,1,2,3,4,\dots$ are

$$\frac{1}{2},\frac{1}{4},\frac{1}{8},\frac{1}{16},\frac{1}{32},\dots ,$$

(13)

namely, there is $50\%$ occupation of the first natural geminal, $25\%$ occupation of the second, $12.5\%$ of the third, $6.25\%$ of the fourth, $3.125\%$ of the fifth, and so on. For brevity, we refer to the fragmentation values in (13) as $50\%$ fragmentation.

Now, one can compute for which ratio $\frac{{\Omega }_{12}}{\omega }$, or, equivalently, for which inter-species interaction ${\lambda }_{12}=\frac{m{\omega }^2}{4N}\left[{\left(\frac{{\Omega }_{12}}{\omega }\right)}^2-1\right]$, the intra-species reduced two-particle and one-particle density matrices are macroscopically fragmented as in (13). Thus, solving (12a) for $50\%$ natural-orbital fragmentation, we find

$$\begin{array}{c}\hfill {\rho }_1^{\left(1\right)}=\frac{1}{2}⟹\frac{{\Omega }_{12}}{\omega }=\left(1+\frac{8N^2}{N-\frac{1}{2}}\right)\pm \sqrt{{\left(1+\frac{8N^2}{N-\frac{1}{2}}\right)}^2-1},\end{array}$$

(14)

and working out (12b) for $50\%$ natural-geminal fragmentation, we obtain

$$\begin{array}{c}\hfill {\rho }_1^{\left(2\right)}=\frac{1}{2}⟹\frac{{\Omega }_{12}}{\omega }=\sqrt{1+\frac{4N}{m{\omega }^2}{\lambda }_{12}}=\left(1+\frac{4N^2}{N-1}\right)\pm \sqrt{{\left(1+\frac{4N^2}{N-1}\right)}^2-1}.\end{array}$$

(15)

There are two ‘reciprocate’ solutions for both the natural geminals and natural orbitals: We see that $50\%$ fragmentation occurs for strong attractions, i.e., when $\frac{{\Omega }_{12}}{\omega }$ is large, or near the border of stability for repulsions, which is when $\frac{{\Omega }_{12}}{\omega }$ is close to zero. In addition, to achieve the same degree of $50\%$ with a larger number N of species 1 (and species 2) bosons, a stronger attraction or repulsion is needed. Finally, comparing the natural-geminal with the natural-orbital fragmentation at the same $50\%$ value, one sees from (15) and (14) that slightly weaker interactions, attractions or repulsions, are needed for the former.

It is also useful to register the one-particle and two-particle densities, i.e., the diagonal parts ${\rho }_1^{\left(1\right)}\left(x\right)={\rho }_1^{\left(1\right)}(x,x^{\prime }=x)$ and ${\rho }_1^{\left(2\right)}(x_1,x_2)={\rho }_1^{\left(2\right)}(x_1,x_2,x_1^{\prime }=x_1,x_2^{\prime }=x_2)$, which read

$$\begin{array}{cc}\hfill & {\rho }_1^{\left(1\right)}\left(x\right)=N{\left(\frac{\alpha +C_{1,0}}{\pi }\right)}^{\frac{1}{2}}{e}^{-\left(\alpha +C_{1,0}\right)x^2}=N{\left(\frac{m\omega }{\pi \left[1+\frac{1}{2N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}\right)}^{\frac{1}{2}}{e}^{-\frac{m\omega }{1+\frac{1}{2N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}{x}^2},\hfill \\ \hfill & {\rho }_1^{\left(2\right)}(x_1,x_2)=N(N-1){\left(\frac{\alpha -\beta }{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{\alpha -\beta }{2}{\left(x_1-x_2\right)}^2}{\left(\frac{\alpha +\beta +2C_{2,0}}{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{\alpha +\beta +2C_{2,0}}{2}{\left(x_1+x_2\right)}^2}\hfill \\ \hfill & =N(N-1){\left(\frac{m\omega }{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{m\omega }{2}{\left(x_1-x_2\right)}^2}{\left(\frac{m\omega }{\pi \left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}\right)}^{\frac{1}{2}}{e}^{-\frac{m\omega }{2\left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}{\left(x_1+x_2\right)}^2}.\hfill \end{array}$$

(16)

From the densities (16), we can infer a measure for the size of identical pairs’ and bosons’ clouds using the widths of the respective Gaussian functions therein. Thus, we have

$$\begin{array}{cc}\hfill & {\sigma }_{1,x}^{\left(1\right)}=\sqrt{\frac{1+\frac{1}{2N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}{2m\omega }},\hfill \\ \hfill & {\sigma }_{1,\frac{x_1+x_2}{\sqrt{2}}}^{\left(2\right)}=\sqrt{\frac{1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}{2m\omega }},{\sigma }_{1,\frac{x_1-x_2}{\sqrt{2}}}^{\left(2\right)}=\sqrt{\frac{1}{2m\omega }}.\hfill \end{array}$$

(17a)

We can understand intuitively why the width ${\sigma }_{1,\frac{x_1-x_2}{\sqrt{2}}}^{\left(2\right)}$ is the typical length of the harmonic trap. For this, recall that ${\sigma }_{1,\frac{x_1-x_2}{\sqrt{2}}}^{\left(2\right)}$ is associated with the relative two-boson coordinate $q_2,q_2^{\prime }$. The general intra-species two-particle reduced density matrix (9) is diagonal in this relative coordinate. The respective exponent is $\alpha -\beta =m\Omega$. Then, since for the specific case studied $\lambda +{\lambda }_{12}=0$ the exponent boils down to that of the bare harmonic potential, i.e., $\Omega =\omega$, we readily obtain the result.

To assess the combined impact of the intra-species and inter-species interactions atop the fragmentation of the reduced density matrices, it is useful to compute the sizes (17a) for large inter-species attractions or inter-species repulsions at the border of stability. One finds, respectively,

$$\begin{array}{cc}\hfill & \underset{\frac{{\Omega }_{12}}{\omega }\to \infty }{lim}{\sigma }_{1,x}^{\left(1\right)}=\sqrt{\frac{1-\frac{1}{2N}}{2m\omega }},{\sigma }_{1,x}^{\left(1\right)}⟶\infty \mathrm{for}\frac{{\Omega }_{12}}{\omega }\to 0^{+},\hfill \\ \hfill & \underset{\frac{{\Omega }_{12}}{\omega }\to \infty }{lim}{\sigma }_{1,\frac{x_1+x_2}{\sqrt{2}}}^{\left(2\right)}=\sqrt{\frac{1-\frac{1}{N}}{2m\omega }},{\sigma }_{1,\frac{x_1+x_2}{\sqrt{2}}}^{\left(2\right)}⟶\infty \mathrm{for}\frac{{\Omega }_{12}}{\omega }\to 0^{+},\hfill \end{array}$$

(17b)

where ${\sigma }_{1,\frac{x_1-x_2}{\sqrt{2}}}^{\left(2\right)}$ as discussed above is independent of the interactions. Interestingly, the size of the densities for strong inter-species attractions, which is accompanied by strong intra-species repulsions because $\lambda +{\lambda }_{12}=0$, saturates at about the trap’s size and does not depend on the strengths of interactions. In other words, a high degree of fragmentation is possible in the mixture without shrinking of the density due to strong inter-species attractive interaction or expansion of the intra-species densities due to strong intra-species repulsive interaction. For the sake of comparative analysis, it is instructive to make contact with the fragmentation of single-species bosons in the harmonic-interaction model; see Appendix A.

2.3. Inter-Species Natural Pair Functions

As mentioned above, in a mixture of two types of identical bosons, there are other kinds of pairs, namely, pairs of distinguishable particles. If we are to examine the lowest-order inter-species reduced density matrix, we can ask regarding distinguishable pairs questions analogous to those asked concerning identical pairs. The purpose of this subsection is to derive the relevant tools and answer such questions.

The inter-species reduced two-particle density matrix, i.e., the lowest-oder inter-species quantity, is defined from the all-particle density matrix as

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(2\right)}(x,x^{\prime },y,y^{\prime })=N^2\int d{x}_2\cdots d{x}_{N}d{y}_2\cdots d{y}_N\Psi (x,x_2,\dots ,x_N,y,y_2,\dots ,y_N)\hfill \\ \hfill & \times {\Psi }^{*}(x^{\prime },x_2,\dots ,x_N,y^{\prime },y_2,\dots ,y_N).\hfill \end{array}$$

(18)

For the harmonic-interaction model of the symmetric mixture, it can be computed analytically and, starting from (2b), is given by [74]

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(2\right)}(x,x^{\prime },y,y^{\prime })=N^2{\left[\frac{{({\alpha }_1+C_{1,1})}^2-D_{1,1}^2}{{\pi }^2}\right]}^{\frac{1}{2}}\hfill \\ \hfill & \times e^{-\frac{{\alpha }_1}{2}\left(x^2+{x^{\prime }}^2+y^2+{y^{\prime }}^2\right)}{e}^{-\frac{1}{4}{C}_{1,1}\left[{\left(x+x^{\prime }\right)}^2+{\left(y+y^{\prime }\right)}^2\right]}{e}^{+\frac{1}{2}{D}_{1,1}\left(x+x^{\prime }\right)\left(y+y^{\prime }\right)}{e}^{+\frac{1}{2}{D}_{1,1}^{\prime }\left(x-x^{\prime }\right)\left(y-y^{\prime }\right)},\hfill \end{array}$$

(19a)

where

$$\begin{array}{cc}\hfill & \alpha +C_{1,1}\mp D_{1,1}=(\alpha -\beta )\frac{(\alpha -\beta )+N(\beta \mp \gamma )}{(\alpha -\beta )+(N-1)(\beta \mp \gamma )},\hfill \\ \hfill & D_{1,1}^{\prime }=\gamma .\hfill \end{array}$$

(19b)

We see that the structure of the inter-species reduced two-particle density matrix is more involved than that of the intra-species reduced two-particle density matrix, as well as that of the product of the two, species 1 and species 2 intra-species reduced one-particle density matrices. Nonetheless, it can be diagonalized.

To diagonalize ${\rho }_{12}^{\left(2\right)}(x,x^{\prime },y,y^{\prime })$, one must couple and make linear combinations of coordinates associated with distinguishable bosons. Defining $u=\frac{1}{\sqrt{2}}\left(x+y\right)$, $v=\frac{1}{\sqrt{2}}\left(x-y\right)$ and $u^{\prime }=\frac{1}{\sqrt{2}}\left(x^{\prime }+y^{\prime }\right)$, $v^{\prime }=\frac{1}{\sqrt{2}}\left(x^{\prime }-y^{\prime }\right)$, we have for the different terms in (19):

$$\begin{array}{cc}\hfill & x^2+y^2+{x^{\prime }}^2+{y^{\prime }}^2=u^2+{u^{\prime }}^2+v^2+{v^{\prime }}^2,\hfill \\ \hfill & {\left(x+x^{\prime }\right)}^2+{\left(y+y^{\prime }\right)}^2={\left(u+u^{\prime }\right)}^2+{\left(v+v^{\prime }\right)}^2=u^2+{u^{\prime }}^2+v^2+{v^{\prime }}^2+2\left(u{u}^{\prime }+v{v}^{\prime }\right),\hfill \\ \hfill & \left(x\pm x^{\prime }\right)\left(y\pm y^{\prime }\right)=\frac{1}{2}\left[{\left(u\pm u^{\prime }\right)}^2-{\left(v\pm v^{\prime }\right)}^2\right]=\frac{1}{2}\left(u^2+{u^{\prime }}^2-v^2-{v^{\prime }}^2\right)\pm \left(u{u}^{\prime }-v{v}^{\prime }\right).\hfill \end{array}$$

(20)

Consequently, we readily find the decomposition

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(2\right)}(u,u^{\prime },v,v^{\prime })=N^2{\left(\frac{{\alpha }_1+C_{1,1}-D_{1,1}}{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{{\alpha }_1+\frac{C_{1,1}}{2}-\frac{D_{1,1}+D_{1,1}^{\prime }}{2}}{2}\left(u^2+{u^{\prime }}^2\right)}{e}^{-\frac{1}{2}\left[C_{1,1}-\left(D_{1,1}-D_{1,1}^{\prime }\right)\right]u{u}^{\prime }}\hfill \\ \hfill & \times {\left(\frac{{\alpha }_1+C_{1,1}+D_{1,1}}{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{{\alpha }_1+\frac{C_{1,1}}{2}+\frac{D_{1,1}+D_{1,1}^{\prime }}{2}}{2}\left(v^2+{v^{\prime }}^2\right)}{e}^{-\frac{1}{2}\left[C_{1,1}+\left(D_{1,1}-D_{1,1}^{\prime }\right)\right]v{v}^{\prime }},\hfill \end{array}$$

(21)

where the normalizations after and before diagonalization are, of course, equal. As might be expected, since the structure of ${\rho }_{12}^{\left(2\right)}(x,x^{\prime },y,y^{\prime })$ is more involved than that of ${\rho }_1^{\left(2\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime })$, the diagonalization of the former is more intricate. Fortunately, we can do that using the application of Mehler’s formula twice, on the appropriately-constructed inter-species ‘mixed’ coordinates $u,u^{\prime }$ and $v,v^{\prime }$. We thus obtain

$$\begin{array}{cc}\hfill & s_{12,\pm }^{\left(2\right)}=\sqrt{\left(\alpha \mp D_{1,1}^{\prime }\right)\left(\alpha +C_{1,1}\mp D_{1,1}\right)}=\sqrt{(\alpha \mp \gamma )(\alpha -\beta )\frac{(\alpha -\beta )+N(\beta \mp \gamma )}{(\alpha -\beta )+(N-1)(\beta \mp \gamma )}},\hfill \\ \hfill & {\rho }_{12,\pm }^{\left(2\right)}=\frac{\left(\alpha \mp D_{1,1}^{\prime }\right)-s_{12,\pm }^{\left(2\right)}}{\left(\alpha \mp D_{1,1}^{\prime }\right)+s_{12,\pm }^{\left(2\right)}}=\frac{\frac{(\alpha \mp \gamma )\left[(\alpha -\beta )+(N-1)(\beta \mp \gamma )\right]}{(\alpha -\beta )\left[(\alpha -\beta )+N(\beta \mp \gamma )\right]}-1}{\frac{(\alpha \mp \gamma )\left[(\alpha -\beta )+(N-1)(\beta \mp \gamma )\right]}{(\alpha -\beta )\left[(\alpha -\beta )+N(\beta \mp \gamma )\right]}+1},\hfill \\ \hfill & 1-{\rho }_{12,\pm }^{\left(2\right)}=\frac{2s_{12,\pm }^{\left(2\right)}}{\left(\alpha \mp D_{1,1}^{\prime }\right)+s_{12,\pm }^{\left(2\right)}},\hfill \end{array}$$

(22)

where the “+” terms quantify the fragmentation in the $u,u^{\prime }$ part of the inter-species reduced two-particle density matrix and the “−” terms quantify the fragmentation in the $v,v^{\prime }$ part of the inter-species reduced two-particle density matrix; also see below. Equation (22) is one of the main results of the present work and bears a clear and appealing physical meaning: that pairs made of distinguishable bosons can be fragmented and that this fragmentation is governed by the center-of-mass and by a relative coordinate of distinguishable bosons. We shall return to this point in what follows.

We can now prescribe the decomposition of the inter-species reduced two-particle density matrix to its distinguishable natural pair functions, which is given by

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(2\right)}(x,x^{\prime },y,y^{\prime })\hfill \\ \hfill & =N^2{\sum }_{n_{+}=0}^{\infty }{\sum }_{n_{-}=0}^{\infty }\left(1-{\rho }_{12,+}^{\left(2\right)}\right)\left(1-{\rho }_{12,-}^{\left(2\right)}\right){\left({\rho }_{12,+}^{\left(2\right)}\right)}^{n_{+}}{\left({\rho }_{12,-}^{\left(2\right)}\right)}^{n_{-}}{\Phi }_{12,n_{+},n_{-}}^{\left(2\right)}(x,y){\Phi }_{12,n_{+},n_{-}}^{\left(2\right),*}(x^{\prime },y^{\prime }),\hfill \\ \hfill & {\Phi }_{12,n_{+},n_{-}}^{\left(2\right)}(x,y)=\frac{1}{\sqrt{2^{n_{+}}{n}_{+}!}}{\left(\frac{s_{12,+}^{\left(2\right)}}{\pi }\right)}^{\frac{1}{4}}{H}_{n_{+}}\left(\sqrt{\frac{s_{12,+}^{\left(2\right)}}{2}}\left(x+y\right)\right)e^{-\frac{1}{4}{s}_{12,+}^{\left(2\right)}{\left(x+y\right)}^2}\hfill \\ \hfill & \times \frac{1}{\sqrt{2^{n_{-}}{n}_{-}!}}{\left(\frac{s_{12,-}^{\left(2\right)}}{\pi }\right)}^{\frac{1}{4}}{H}_{n_{-}}\left(\sqrt{\frac{s_{12,-}^{\left(2\right)}}{2}}\left(x-y\right)\right)e^{-\frac{1}{4}{s}_{12,-}^{\left(2\right)}{\left(x-y\right)}^2}.\hfill \end{array}$$

(23)

All in all, (23) implies that the distinguishable-pair ‘condensed fraction’ is given by $\left(1-{\rho }_{12,+}^{\left(2\right)}\right)\left(1-{\rho }_{12,-}^{\left(2\right)}\right)$ and the respective depleted fraction by $1-\left(1-{\rho }_{12,+}^{\left(2\right)}\right)\left(1-{\rho }_{12,-}^{\left(2\right)}\right)={\rho }_{12,+}^{\left(2\right)}+{\rho }_{12,-}^{\left(2\right)}-{\rho }_{12,+}^{\left(2\right)}{\rho }_{12,-}^{\left(2\right)}$. Each of the inter-species ‘mixed’ coordinates $\frac{x\pm y}{\sqrt{2}}$ carries its own scaling, $s_{12,\pm }^{\left(2\right)}$. The distinguishable natural geminals ${\Phi }_{12,n_{+},n_{-}}^{\left(2\right)}(x,y)$ are, needless to say, orthonormal to each other.

We proceed now for an application. We considered above the specific case of $\lambda +{\lambda }_{12}=0$, which leads to $\Omega =\omega$, $\alpha =m\omega \left[1+\frac{1}{2N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]$, and $\beta =\gamma =\frac{m}{2N}\left({\Omega }_{12}-\omega \right)$. Recall that this implies that the frequency of $2N-1$ coordinates, the $2(N-1)$ relative coordinates, and the center-of-mass coordinate is $\omega$ and that only the relative center-of-mass coordinate has an interaction-dressed frequency, ${\Omega }_{12}$. To evaluate ${\rho }_{12}^{\left(2\right)}(x,x^{\prime },y,y^{\prime })$, we also need the combinations $\left(\alpha +C_{1,1}-D_{1,1}\right)=m\omega$ and $\left(\alpha -D_{1,1}^{\prime }\right)=m\omega$ for the “+” branch, as well as $\left(\alpha +C_{1,1}+D_{1,1}\right)=m\omega \frac{1}{1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}$ and $\left(\alpha +D_{1,1}^{\prime }\right)=m\omega \left[1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]$ for the “−” branch.

Thus, expressions (22) can readily be evaluated, and the following picture of inter-species fragmentation is found:

$$\begin{array}{ccc}\hfill & & s_{12,+}^{\left(2\right)}=m\omega ,{\rho }_{12,+}^{\left(2\right)}=0,1-{\rho }_{12,+}^{\left(2\right)}=1,\hfill \end{array}$$

(24a)

indicating that there is no contribution to fragmentation from the symmetric ‘mixed’ coordinate $u,u^{\prime }$. On the other end,

$$\begin{array}{cc}\hfill & s_{12,-}^{\left(2\right)}=m\omega \sqrt{\frac{1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)}{1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}},\hfill \\ \hfill & {\rho }_{12,-}^{\left(2\right)}=\frac{\sqrt{\left[1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}-1}{\sqrt{\left[1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}+1},\hfill \\ \hfill & 1-{\rho }_{12,-}^{\left(2\right)}=\frac{2}{\sqrt{\left[1+\frac{1}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}+1},\hfill \end{array}$$

(24b)

namely, that the fragmentation fully originates from the asymmetric ‘mixed’ coordinate $v,v^{\prime }$. We conclude that, whereas the fragmentation of identical pairs is associated with their center-of-mass coordinate, the fragmentation of distinguishable pairs is linked, in this explicit case, only with a relative coordinate between two distinguishable bosons. Interestingly, the degree of intra-species and inter-species pair fragmentation is the same in the specific case considered, despite pertaining to different parts of the mixtures’ many-boson wave function. Furthermore, there are different numbers of pairs: $\frac{N}{2}$ intra-species identical pairs (for each of the species) and N inter-species pairs of distinguishable bosons.

Now, one can compute the ratio $\frac{{\Omega }_{12}}{\omega }=\sqrt{1+\frac{4N}{m{\omega }^2}{\lambda }_{12}}$ for which the inter-species reduced two-particle density matrix is $50\%$ fragmented as in (13). Since ${\rho }_{12,+}^{\left(2\right)}=0$ does not contribute, the only contribution to fragmentation comes from ${\rho }_{12,-}^{\left(2\right)}$. Thus, solving (24b) for $50\%$ distinguishable-pair-function fragmentation, we obtain

$$\begin{array}{c}\hfill {\rho }_{12,-}^{\left(2\right)}=\frac{1}{2}⟹\frac{{\Omega }_{12}}{\omega }=\left(1+\frac{4N^2}{N-1}\right)\pm \sqrt{{\left(1+\frac{4N^2}{N-1}\right)}^2-1}.\end{array}$$

(25)

As above, there are two ‘reciprocate’ solutions: one for strong inter-species attraction and the second close to the border of stability for intermediate-strength inter-species repulsion. We remind that since $\lambda +{\lambda }_{12}=0$ in our example, the respective intra-species interaction is opposite in sign. In addition, to achieve the same degree of $50\%$ fragmentation with a larger number N of distinguishable pairs, a stronger inter-species attraction or repulsion is needed. Furthermore, as discussed above, comparing distinguishable-pair and identical-pair fragmentation at the same $50\%$ value in this example, one sees from (25) and (15) that the same interaction is needed.

Finally, we prescribe the inter-species two-particle density, namely, the diagonal part ${\rho }_{12}^{\left(2\right)}(x,y)={\rho }_{12}^{\left(2\right)}(x,x^{\prime }=x,y,y^{\prime }=y)$, which is given by

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(2\right)}(x,y)=N^2{\left(\frac{{\alpha }_1+C_{1,1}-D_{1,1}}{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{{\alpha }_1+C_{1,1}-D_{1,1}}{2}{\left(x+y\right)}^2}\hfill \\ \hfill & \times {\left(\frac{{\alpha }_1+C_{1,1}+D_{1,1}}{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{{\alpha }_1+C_{1,1}+D_{1,1}}{2}{\left(x-y\right)}^2}\hfill \\ \hfill & =N^2{\left(\frac{m\omega }{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{m\omega }{2}{\left(x+y\right)}^2}{\left(\frac{m\omega }{\pi \left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}\right)}^{\frac{1}{2}}{e}^{-\frac{m\omega }{2\left[1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}{\left(x-y\right)}^2}.\hfill \end{array}$$

(26)

Next, the size of the distinguishable pairs’ cloud can be assessed from the density (26) using the widths of the respective Gaussian functions. Accordingly, we find

$$\begin{array}{ccc}\hfill & & {\sigma }_{12,\frac{x+y}{\sqrt{2}}}^{\left(2\right)}=\sqrt{\frac{1}{2m\omega }},{\sigma }_{12,\frac{x-y}{\sqrt{2}}}^{\left(2\right)}=\sqrt{\frac{1+\frac{1}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}{2m\omega }}.\hfill \end{array}$$

(27a)

The explanation why the width ${\sigma }_{12,\frac{x+y}{\sqrt{2}}}^{\left(2\right)}$ is also the typical length of the harmonic potential is somewhat less intuitive. The distinguishable pair center-of-mass coordinate $u,u^{\prime }$ appears in the general inter-species reduced two-particle density matrix (21) in a coupled form which requires diagonalization via Mehler’s formula. We find in the specific case $\lambda +{\lambda }_{12}=0$ that $\gamma =\beta$. This leads to a neat cancellation of terms, and the retaining of a single term in the Mehler’s expansion with the typical length of the harmonic potential.

To show the combined effect of the intra-species and inter-species interactions accompanying the fragmentation of ${\rho }_{12}^{\left(2\right)}(x,x^{\prime },y,y^{\prime })$, it is useful to compute the sizes (27a) for large inter-species attractions or inter-species repulsions at the border of stability. We obtain, respectively,

$$\begin{array}{ccc}\hfill & & \underset{\frac{{\Omega }_{12}}{\omega }\to \infty }{lim}{\sigma }_{12,\frac{x-y}{\sqrt{2}}}^{\left(2\right)}=\sqrt{\frac{1-\frac{1}{N}}{2m\omega }},{\sigma }_{12,\frac{x-y}{\sqrt{2}}}^{\left(2\right)}⟶\infty \mathrm{for}\frac{{\Omega }_{12}}{\omega }\to 0^{+},\hfill \end{array}$$

(27b)

where ${\sigma }_{12,\frac{x+y}{\sqrt{2}}}^{\left(2\right)}$ as discussed above is independent of the interactions. We see that the size of the inter-species density saturates as does the trap’s size and does not depend on the strengths of interactions in the limit of strong inter-species attractions. Analogously to identical pairs, a strong fragmentation of distinguishable pairs is possible in the mixture without the shrinking of the inter-species density due to strong inter-species attractive interaction. At the other end, when the inter-species repulsion is close to the border of stability, the inter-species density expands boundlessly. Summarizing, inter-species fragmentation is governed by the ratio $\frac{{\Omega }_{12}}{\omega }$ and takes place both at the attractive and repulsive sectors of interactions. For the sake of analysis, we compared the results for inter-species pair fragmentation with intra-species pair fragmentation and discussed the similarity and differences between the respective two-particle densities (16) and (26).

3. Pair of Distinguishable Pairs and Schmidt Decomposition of the wave function

Following the results of the previous section on the fragmentation of distinguishable pairs, there are two questions that warrant answers. The first is whether inter-species fragmentation persists beyond distinguishable pairs, say, to pairs of distinguishable pairs? Inasmuch as single-species and intra-species fragmentations take place at the lowest-level reduced one-particle density matrix and persist at higher-level single-species reduced density matrices, we wish to establish the result of inter-species fragmentation at the level of higher-order reduced density matrices. After all, the reduced two-particle density matrix is the lowest-order inter-species one. The second question deals with the nature of the inter-species coordinates governing fragmentation. At the level of distinguishable-pair fragmentation, i.e., within the inter-species reduced two-particle density matrix, one cannot unambiguously tell whether the relative center-of-mass coordinate of the two species is involved or whether other relative inter-species coordinates govern fragmentation. This is because in a pair of distinguishable particles, one cannot distinguish between the two types of coordinates.

As seen in the previous section, the inter-species reduced two-particle density matrix is more intricate than the intra-species ones, and consequently, its diagonalization is more involved. We derive now the inter-species reduced four-particle density matrix and examine which ‘normal coordinates’ govern its diagonalization. Then, the natural four-particle functions are obtained explicitly and investigated.

Finally, and as a complementary result of the techniques used for inter-species fragmentation, we carry the connection between inter-species and intra-species center-of-mass coordinates, in conjunction with the usage of Mehler’s formula within a mixture, further. This is performed by constructing the Schmidt decomposition of the mixture’s wave function and discussing the consequences of this decomposition at the limit of an infinite number of particles.

3.1. Inter-Species Fragmentation in Higher-Order Reduced Density Matrices

The inter-species reduced four-particle density matrix is defined as

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(4\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime },y_1,y_2,y_1^{\prime },y_2^{\prime })=N^2{\left(N-1\right)}^2\int d{x}_3\cdots d{x}_{N}d{y}_3\cdots d{y}_N\hfill \\ \hfill & \times \Psi (x_1,x_2,x_3,\dots ,x_N,y_1,y_2,y_3,\dots ,y_N){\Psi }^{*}(x_1^{\prime },x_2^{\prime },x_3,\dots ,x_N,y_1^{\prime },y_2^{\prime },y_3,\dots ,y_N).\hfill \end{array}$$

(28)

Note that here we only treat the four-particle quantity with two identical bosons per each species. Integrating the harmonic interaction-model for symmetric mixtures, we find the final expression explicitly

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(4\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime },y_1,y_2,y_1^{\prime },y_2^{\prime })=N^2{\left(N-1\right)}^2{\left[\frac{{(\alpha +C_{1,1})}^2-D_{1,1}^2}{{\pi }^2}\right]}^{\frac{1}{2}}{\left[\frac{{(\alpha +C_{2,2})}^2-D_{2,2}^2}{{\pi }^2}\right]}^{\frac{1}{2}}\hfill \\ \hfill & \times e^{-\frac{\alpha }{2}\left(x_1^2+x_2^2+{x_1^{\prime }}^2+{x_2^{\prime }}^2+y_1^2+y_2^2+{y_1^{\prime }}^2+{y_2^{\prime }}^2\right)}{e}^{-\beta \left(x_1{x}_2+x_1^{\prime }{x}_2^{\prime }+y_1{y}_2+y_1^{\prime }{y}_2^{\prime }\right)}\hfill \\ \hfill & \times e^{-\frac{1}{4}{C}_{2,2}\left[{\left(x_1+x_2+x_1^{\prime }+x_2^{\prime }\right)}^2+{\left(y_1+y_2+y_1^{\prime }+y_2^{\prime }\right)}^2\right]}\hfill \\ \hfill & \times e^{+\frac{1}{2}{D}_{2,2}\left(x_1+x_2+x_1^{\prime }+x_2^{\prime }\right)\left(y_1+y_2+y_1^{\prime }+y_2^{\prime }\right)}{e}^{+\frac{1}{2}{D}_{2,2}^{\prime }\left(x_1+x_2-x_1^{\prime }-x_2^{\prime }\right)\left(y_1+y_2-y_1^{\prime }-y_2^{\prime }\right)},\hfill \end{array}$$

(29a)

where

$$\begin{array}{cc}\hfill & \alpha +\beta +2\left(C_{2,2}\mp D_{2,2}\right)=(\alpha -\beta )\frac{(\alpha -\beta )+N(\beta \mp \gamma )}{(\alpha -\beta )+(N-2)(\beta \mp \gamma )},\hfill \\ \hfill & D_{2,2}^{\prime }=\gamma ,\hfill \end{array}$$

(29b)

and $\alpha +C_{1,1}\mp D_{1,1}$ are given in (19b). The combinations of parameters $\alpha +\beta +2\left(C_{2,2}\mp D_{2,2}\right)$ would appear below shortly.

To diagonalize ${\rho }_{12}^{\left(4\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime },y_1,y_2,y_1^{\prime },y_2^{\prime })$, we need to mix and rotate the coordinates of the two species into new coordinates appropriately. Thus, defining the new coordinates as the center-of-mass, relative center-of-mass, and relative coordinates of two identical pairs, one for each of the species, $u_1=\frac{1}{2}\left[\left(x_1+x_2\right)+\left(y_1+y_2\right)\right]$, $v_1=\frac{1}{2}\left[\left(x_1+x_2\right)-\left(y_1+y_2\right)\right]$, $u_2=\frac{1}{\sqrt{2}}\left(x_1-x_2\right)$, $v_2=\frac{1}{\sqrt{2}}\left(y_1-y_2\right)$ and $u_1^{\prime }=\frac{1}{2}\left[\left(x_1^{\prime }+x_2^{\prime }\right)\right.+\left.\left(y_1^{\prime }+y_2^{\prime }\right)\right]$, $v_1^{\prime }=\frac{1}{2}[\left(x_1^{\prime }+x_2^{\prime }\right)$$-\left(y_1^{\prime }+y_2^{\prime }\right)]$, $u_2^{\prime }=\frac{1}{\sqrt{2}}\left(x_1^{\prime }-x_2^{\prime }\right)$, $v_2^{\prime }=\frac{1}{\sqrt{2}}\left(y_1^{\prime }-y_2^{\prime }\right)$, we have for the different terms in (29a):

$$\begin{array}{cc}\hfill & x_1^2+x_2^2+y_1^2+y_2^2+{x_1^{\prime }}^2+{x_2^{\prime }}^2+{y_1^{\prime }}^2+{y_2^{\prime }}^2=u_1^2+{u_1^{\prime }}^2+v_1^2+{v_1^{\prime }}^2+u_2^2+{u_2^{\prime }}^2+v_2^2+{v_2^{\prime }}^2,\hfill \\ \hfill & x_1{x}_2+x_1^{\prime }{x}_2^{\prime }+y_1{y}_2+y_1^{\prime }{y}_2^{\prime }=\frac{1}{2}\left(u_1^2+{u_1^{\prime }}^2+v_1^2+{v_1^{\prime }}^2-u_2^2-{u_2^{\prime }}^2-v_2^2-{v_2^{\prime }}^2\right),\hfill \\ \hfill & {\left(x_1+x_2+x_1^{\prime }+x_2^{\prime }\right)}^2+{\left(y_1+y_2+y_1^{\prime }+y_2^{\prime }\right)}^2=2\left[{\left(u_1+u_1^{\prime }\right)}^2+{\left(v_1+v_1^{\prime }\right)}^2\right],\hfill \\ \hfill & \left[\left(x_1+x_2\right)\pm \left(x_1^{\prime }+x_2^{\prime }\right)\right]\left[\left(y_1+y_2\right)\pm \left(y_1^{\prime }+y_2^{\prime }\right)\right]={\left(u_1\pm u_1^{\prime }\right)}^2-{\left(v_1\pm v_1^{\prime }\right)}^2.\hfill \end{array}$$

(30)

Relations (30) imply that one could equally define inter-species linear combinations of the relative coordinates, since $u_2^2+v_2^2={\left[\frac{\left(x_1-x_2\right)+\left(y_1-y_2\right)}{2}\right]}^2+{\left[\frac{\left(x_1-x_2\right)-\left(y_1-y_2\right)}{2}\right]}^2$ and ${u_2^{\prime }}^2+{v_2^{\prime }}^2={\left[\frac{\left(x_1^{\prime }-x_2^{\prime }\right)+\left(y_1^{\prime }-y_2^{\prime }\right)}{2}\right]}^2+{\left[\frac{\left(x_1^{\prime }-x_2^{\prime }\right)-\left(y_1^{\prime }-y_2^{\prime }\right)}{2}\right]}^2$. We chose the former combinations.

Plugging (30) into (29), we readily find for the transformed inter-species reduced four-particle density matrix

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(4\right)}(u_1,u_1^{\prime },v_1,v_1^{\prime },u_2,u_2^{\prime },v_2,v_2^{\prime })=N^2{\left(N-1\right)}^2\hfill \\ \hfill & \times {\left(\frac{\alpha -\beta }{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{\alpha -\beta }{2}\left(u_2^2+{u_2^{\prime }}^2\right)}{\left(\frac{\alpha -\beta }{\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{\alpha -\beta }{2}\left(v_2^2+{v_2^{\prime }}^2\right)}\hfill \\ \hfill & \times {\left[\frac{\alpha +\beta +2\left(C_{2,2}-D_{2,2}\right)}{\pi }\right]}^{\frac{1}{2}}{e}^{-\frac{\alpha +\beta +C_{2,2}-\left(D_{2,2}+D_{2,2}^{\prime }\right)}{2}\left(u_1^2+{u_1^{\prime }}^2\right)}{e}^{-\left[C_{2,2}-\left(D_{2,2}-D_{2,2}^{\prime }\right)\right]u_1{u}_1^{\prime }}\hfill \\ \hfill & \times {\left[\frac{\alpha +\beta +2\left(C_{2,2}+D_{2,2}\right)}{\pi }\right]}^{\frac{1}{2}}{e}^{-\frac{\alpha +\beta +C_{2,2}+\left(D_{2,2}+D_{2,2}^{\prime }\right)}{2}\left(v_1^2+{v_1^{\prime }}^2\right)}{e}^{-\left[C_{2,2}+\left(D_{2,2}-D_{2,2}^{\prime }\right)\right]v_1{v}_1^{\prime }},\hfill \end{array}$$

(31)

where the normalization coefficients before and after diagonalization are, naturally, equal and fulfill $\left[\alpha +\left(C_{1,1}\mp D_{1,1}\right)\right]\left[\alpha +\left(C_{2,2}\mp D_{2,2}\right)\right]=\left(\alpha -\beta \right)\left[\alpha +\beta +2\left(C_{2,2}\mp D_{2,2}\right)\right]$.

As can be seen in (31) and (21), the similarities and differences between the structures of ${\rho }_{12}^{\left(4\right)}(u_1,u_1^{\prime },v_1,v_1^{\prime },u_2,u_2^{\prime },v_2,v_2^{\prime })$ and ${\rho }_{12}^{\left(2\right)}(u,u^{\prime },v,v^{\prime })$ clarify the issue of which coordinates are coupled and identify the coordinates that are not. In particular, just like for the two-particle quantity, we can apply Mehler’s formula twice, on the appropriately constructed inter-species ‘mixed coordinates’ $u_1,u_1^{\prime }$ and $v_1,v_1^{\prime }$, to diagonalize the inter-species reduced four-particle density matrix. When this is performed, one obtains

$$\begin{array}{cc}\hfill & s_{12,\pm }^{\left(4\right)}=\sqrt{\left(\alpha +\beta \mp 2D_{2,2}^{\prime }\right)\left[\alpha +\beta +2\left(C_{2,2}\mp D_{2,2}\right)\right]}\hfill \\ \hfill & =\sqrt{(\alpha +\beta \mp 2\gamma )(\alpha -\beta )\frac{(\alpha -\beta )+N(\beta \mp \gamma )}{(\alpha -\beta )+(N-2)(\beta \mp \gamma )}},\hfill \\ \hfill & {\rho }_{12,\pm }^{\left(4\right)}=\frac{\left(\alpha +\beta \mp 2D_{2,2}^{\prime }\right)-s_{12,\pm }^{\left(4\right)}}{\left(\alpha +\beta \mp 2D_{2,2}^{\prime }\right)+s_{12,\pm }^{\left(4\right)}}=\frac{\frac{(\alpha +\beta \mp 2\gamma )\left[(\alpha -\beta )+(N-2)(\beta \mp \gamma )\right]}{(\alpha -\beta )\left[(\alpha -\beta )+N(\beta \mp \gamma )\right]}-1}{\frac{(\alpha +\beta \mp 2\gamma )\left[(\alpha -\beta )+(N-2)(\beta \mp \gamma )\right]}{(\alpha -\beta )\left[(\alpha -\beta )+N(\beta \mp \gamma )\right]}+1},\hfill \\ \hfill & 1-{\rho }_{12,\pm }^{\left(4\right)}=\frac{2s_{12,\pm }^{\left(4\right)}}{\left(\alpha +\beta \mp 2D_{2,2}^{\prime }\right)+s_{12,\pm }^{\left(4\right)}},\hfill \end{array}$$

(32)

where the “+” terms quantify, the fragmentation in the $u_1,u_1^{\prime }$ part of the inter-species reduced four-particle density matrix and the “−” terms determine the fragmentation in the $v_1,v_1^{\prime }$ part of the inter-species reduced four-particle density matrix. As found and shown in (31), there is no fragmentation due to the relative-coordinate parts $u_2,u_2^{\prime }$ and $v_2,v_2^{\prime }$. Equation (32) adds to the main results of the present work and bears a transparent and appealing physical meaning: in the mixture, inter-species fragmentation is quantified by the eigenvalues obtained from Mehler’s formula when the latter is applied to the mixture’s center-of-mass and relative center-of-mass coordinates of distinguishable pairs of pairs. Extensions to larger distinguishable aggregates of species 1 and species 2 identical bosons in the mixture is possible along the above lines and are not pursued further here.

We can now prescribe the decomposition of the inter-species reduced four-particle density matrix, inasmuch as the reduced two-particle density matrix was decomposed, into its natural four-particle functions made of distinguishable particles. The final result is given by

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(4\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime },y_1,y_2,y_1^{\prime },y_2^{\prime })=N^2{(N-1)}^2{\sum }_{n_{+}=0}^{\infty }{\sum }_{n_{-}=0}^{\infty }\left(1-{\rho }_{12,+}^{\left(4\right)}\right)\left(1-{\rho }_{12,-}^{\left(4\right)}\right)\hfill \\ \hfill & \times {\left({\rho }_{12,+}^{\left(4\right)}\right)}^{n_{+}}{\left({\rho }_{12,-}^{\left(4\right)}\right)}^{n_{-}}{\Phi }_{12,n_{+},n_{-}}^{\left(4\right)}(x_1,x_2,y_1,y_2){\Phi }_{12,n_{+},n_{-}}^{\left(4\right),*}(x_1^{\prime },x_2^{\prime },y_1^{\prime },y_2^{\prime }),\hfill \\ \hfill & {\Phi }_{12,n_{+},n_{-}}^{\left(4\right)}(x_1,x_2,y_1,y_2)\hfill \\ \hfill & =\frac{1}{\sqrt{2^{n_{+}}{n}_{+}!}}{\left(\frac{s_{12,+}^{\left(4\right)}}{\pi }\right)}^{\frac{1}{4}}{H}_{n_{+}}\left(\frac{\sqrt{s_{12,+}^{\left(4\right)}}}{2}\left[\left(x_1+x_2\right)+\left(y_1+y_2\right)\right]\right)e^{-\frac{1}{8}{s}_{12,+}^{\left(4\right)}{\left[\left(x_1+x_2\right)+\left(y_1+y_2\right)\right]}^2}\hfill \\ \hfill & \times \frac{1}{\sqrt{2^{n_{-}}{n}_{-}!}}{\left(\frac{s_{12,-}^{\left(4\right)}}{\pi }\right)}^{\frac{1}{4}}{H}_{n_{-}}\left(\frac{\sqrt{s_{12,-}^{\left(4\right)}}}{2}\left[\left(x_1+x_2\right)-\left(y_1+y_2\right)\right]\right)e^{-\frac{1}{8}{s}_{12,-}^{\left(4\right)}{\left[\left(x_1+x_2\right)-\left(y_1+y_2\right)\right]}^2}\hfill \\ \hfill & \times {\left(\frac{\alpha -\beta }{\pi }\right)}^{\frac{1}{4}}{e}^{-\frac{\alpha -\beta }{4}{\left(x_1-x_2\right)}^2}{\left(\frac{\alpha -\beta }{\pi }\right)}^{\frac{1}{4}}{e}^{-\frac{\alpha -\beta }{4}{\left(y_1-y_2\right)}^2}.\hfill \end{array}$$

(33)

Equation (33) means that the pair-of-distinguishable-pairs ‘condensed fraction’ is given by the product $\left(1-{\rho }_{12,+}^{\left(4\right)}\right)\left(1-{\rho }_{12,-}^{\left(4\right)}\right)$ and the respective depleted fraction is $1-\left(1-{\rho }_{12,+}^{\left(4\right)}\right)\left(1-{\rho }_{12,-}^{\left(4\right)}\right)={\rho }_{12,+}^{\left(4\right)}+{\rho }_{12,-}^{\left(4\right)}-{\rho }_{12,+}^{\left(4\right)}{\rho }_{12,-}^{\left(4\right)}$. The center-of-mass and relative center-of-mass coordinates of the two pairs, $\frac{\left(x_1+x_2\right)\pm \left(y_1+y_2\right)}{2}$, carry the respective scalings $s_{12,\pm }^{\left(4\right)}$. The natural four-particle functions ${\Phi }_{12,n_{+},n_{-}}^{\left(4\right)}(x_1,x_2,y_1,y_2)$ are enumerated by the two quantum numbers $n_{+}$, $n_{-}$ and are obviously orthonormal to each other.

We proceed now to examine the fragmentation in this higher-order inter-species reduced density matrix. We investigate, as mentioned above, the specific case of $\lambda +{\lambda }_{12}=0$. To compute ${\rho }_{12}^{\left(4\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime },y_1,y_2,y_1^{\prime },y_2^{\prime })$, we require the quantities $\left[\alpha +\beta +2\left(C_{2,2}-D_{2,2}\right)\right]$$=m\omega$ and $\left(\alpha +\beta -2D_{2,2}^{\prime }\right)=m\omega$ for the “+” branch, as well as $\left[\alpha +\beta +2\left(C_{2,2}+D_{2,2}\right)\right]=m\omega \frac{1}{1+\frac{2}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}$ and $\left(\alpha +\beta +2D_{2,2}^{\prime }\right)=m\omega \left[1+\frac{2}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]$ for the “−” branch.

Now, expressions (32) can readily be evaluated, and the following picture of higher-order inter-species fragmentation is found:

$$\begin{array}{ccc}\hfill & & s_{12,+}^{\left(4\right)}=m\omega ,{\rho }_{12,+}^{\left(4\right)}=0,1-{\rho }_{12,+}^{\left(4\right)}=1,\hfill \end{array}$$

(34a)

indicating that there is no contribution to fragmentation from the center-of-mass ‘mixed coordinate’ $u_1,u_1^{\prime }$. This is additional to the no contribution to fragmentation coming from the relative coordinates $u_2,u_2^{\prime }$ and $v_2,v_2^{\prime }$; see (31). For the relative center-of-mass ‘mixed coordinate’ $v_1,v_1^{\prime }$, on the other end, one finds

$$\begin{array}{cc}\hfill & s_{12,-}^{\left(4\right)}=m\omega \sqrt{\frac{1+\frac{2}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)}{1+\frac{2}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}},\hfill \\ \hfill & {\rho }_{12,-}^{\left(4\right)}=\frac{\sqrt{\left[1+\frac{2}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{2}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}-1}{\sqrt{\left[1+\frac{2}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{2}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}+1},\hfill \\ \hfill & 1-{\rho }_{12,-}^{\left(4\right)}=\frac{2}{\sqrt{\left[1+\frac{2}{N}\left(\frac{{\Omega }_{12}}{\omega }-1\right)\right]\left[1+\frac{2}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}+1},\hfill \end{array}$$

(34b)

namely, that the fragmentation of pairs of distinguishable pairs fully originates from the relative center-of-mass ‘mixed coordinate’ $v_1,v_1^{\prime }$. We see that also the higher-order inter-species fragmentation is governed by the ratio $\frac{{\Omega }_{12}}{\omega }$ and takes place both at the attractive and repulsive sectors of interactions. In conclusion, the higher-order inter-species fragmentation is proved.

Now, one can compute the ratio $\frac{{\Omega }_{12}}{\omega }=\sqrt{1+\frac{4N}{m{\omega }^2}{\lambda }_{12}}$ for which the inter-species reduced four-particle density matrix is $50\%$ fragmented, as in (13). Since ${\rho }_{12,+}^{\left(4\right)}=0$ does not contribute in this specific case, the only contribution to fragmentation comes from ${\rho }_{12,-}^{\left(4\right)}$. Thus, solving (34b) for $50\%$ distinguishable-four-particle-function fragmentation, we obtain

$$\begin{array}{c}\hfill {\rho }_{12,-}^{\left(4\right)}=\frac{1}{2}⟹\frac{{\Omega }_{12}}{\omega }=\left(1+\frac{2N^2}{N-2}\right)\pm \sqrt{{\left(1+\frac{2N^2}{N-2}\right)}^2-1}.\end{array}$$

(35)

As for distinguishable pairs, there are two ‘reciprocate’ solutions: one for strong attractions and the second close to the border of stability for repulsions. In addition, to achieve the same degree of $50\%$ fragmentation with a larger number $\frac{N}{2}$ of distinguishable four-boson aggregates, a stronger attraction or repulsion is needed. Furthermore, comparing distinguishable-four-boson and distinguishable-two-boson fragmentation at the same $50\%$ value, one sees from (35) and (25) that slightly weaker interactions–attractions or repulsions–are needed for the former. This behavior of the fragmentation of increasing orders of inter-species reduced density matrices is analogous to and generalizes that of intra-species and single-species reduced density matrices; see the previous section and the appendix, respectively.

Finally, we present for completeness the inter-species four-particle density, i.e., the diagonal part ${\rho }_{12}^{\left(4\right)}(x_1,x_2,y_1,y_2)={\rho }_{12}^{\left(4\right)}(x_1,x_2,x_1^{\prime }=x_1,x_2^{\prime }=x_2,y_1,y_2,y_1^{\prime }=y_1,y_2^{\prime }=y_2)$, which is given by

$$\begin{array}{cc}\hfill & {\rho }_{12}^{\left(4\right)}(x_1,x_2,y_1,y_2)=N^2{\left(N-1\right)}^2\left(\frac{\alpha -\beta }{\pi }\right)e^{-\frac{\alpha -\beta }{2}{\left(x_1-x_2\right)}^2}{e}^{-\frac{\alpha -\beta }{2}{\left(y_1-y_2\right)}^2}\hfill \\ \hfill & \times {\left[\frac{\alpha +\beta +2\left(C_{2,2}-D_{2,2}\right)}{\pi }\right]}^{\frac{1}{2}}{e}^{-\frac{\alpha +\beta +2\left(C_{2,2}-D_{2,2}\right)}{4}{\left[\left(x_1+x_2\right)+\left(y_1+y_2\right)\right]}^2}\hfill \\ \hfill & \times {\left[\frac{\alpha +\beta +2\left(C_{2,2}+D_{2,2}\right)}{\pi }\right]}^{\frac{1}{2}}{e}^{-\frac{\alpha +\beta +2\left(C_{2,2}+D_{2,2}\right)}{4}{\left[\left(x_1+x_2\right)-\left(y_1+y_2\right)\right]}^2}\hfill \\ \hfill & =N^2{\left(N-1\right)}^2{\left(\frac{m\omega }{\pi }\right)}^{\frac{3}{2}}{e}^{-\frac{m\omega }{2}{\left(x_1-x_2\right)}^2}{e}^{-\frac{m\omega }{2}{\left(y_1-y_2\right)}^2}{e}^{-\frac{m\omega }{4}{\left[\left(x_1+x_2\right)+\left(y_1+y_2\right)\right]}^2}\hfill \\ \hfill & \times {\left(\frac{m\omega }{\pi \left[1+\frac{2}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}\right)}^{\frac{1}{2}}{e}^{-\frac{m\omega }{4\left[1+\frac{2}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)\right]}{\left[\left(x_1+x_2\right)-\left(y_1+y_2\right)\right]}^2}.\hfill \end{array}$$

(36)

To proceed, the size of the distinguishable four-boson cloud can be estimated from the widths of the respective Gaussian functions in the density (36). Thus, we obtain

$$\begin{array}{cc}\hfill & {\sigma }_{12,\frac{x_1-x_2}{\sqrt{2}}}^{\left(4\right)}=\sqrt{\frac{1}{2m\omega }},{\sigma }_{12,\frac{y_1-y_2}{\sqrt{2}}}^{\left(4\right)}=\sqrt{\frac{1}{2m\omega }},{\sigma }_{12,\frac{\left(x_1+x_2\right)+\left(y_1+y_2\right)}{2}}^{\left(4\right)}=\sqrt{\frac{1}{2m\omega }},\hfill \\ \hfill & {\sigma }_{12,\frac{\left(x_1+x_2\right)-\left(y_1+y_2\right)}{2}}^{\left(4\right)}=\sqrt{\frac{1+\frac{2}{N}\left(\frac{\omega }{{\Omega }_{12}}-1\right)}{2m\omega }}.\hfill \end{array}$$

(37a)

The explanations why the widths ${\sigma }_{12,\frac{x_1-x_2}{\sqrt{2}}}^{\left(4\right)}$, ${\sigma }_{12,\frac{y_1-y_2}{\sqrt{2}}}^{\left(4\right)}$ with intra-species relative coordinates and, separately, the width ${\sigma }_{12,\frac{\left(x_1+x_2\right)-\left(y_1+y_2\right)}{2}}^{\left(4\right)}$ with inter-species center-of-mass coordinate are all the typical length of the harmonic potential follow those given above, alongside Equations (17a) and (27a), respectively, and are not reproduced here.

To show the combined effect of the inter-species and intra-species interactions accompanying the fragmentation of ${\rho }_{12}^{\left(4\right)}(x_1,x_2,x_1^{\prime },x_2^{\prime },y_1,y_2,y_1^{\prime },y_2^{\prime })$, it is instrumental to compute the sizes (37a) for large inter-species attractions or inter-species repulsions at the border of stability. We obtain, respectively,

$$\begin{array}{ccc}\hfill & & \underset{\frac{{\Omega }_{12}}{\omega }\to \infty }{lim}{\sigma }_{12,\frac{\left(x_1+x_2\right)-\left(y_1+y_2\right)}{2}}^{\left(4\right)}=\sqrt{\frac{1-\frac{2}{N}}{2m\omega }},{\sigma }_{12,\frac{\left(x_1+x_2\right)-\left(y_1+y_2\right)}{2}}^{\left(4\right)}⟶\infty \mathrm{for}\frac{{\Omega }_{12}}{\omega }\to 0^{+},\hfill \end{array}$$

(37b)

where ${\sigma }_{12,\frac{x_1-x_2}{\sqrt{2}}}^{\left(4\right)}$, ${\sigma }_{12,\frac{y_1-y_2}{\sqrt{2}}}^{\left(4\right)}$, and ${\sigma }_{12,\frac{\left(x_1+x_2\right)+\left(y_1+y_2\right)}{2}}^{\left(4\right)}$, as mentioned above, do not depend on the interactions. We see that the size of the inter-species four-boson density also saturates at about the trap’s size, and does not depend on the strengths of interactions in the limit of strong inter-species attractions. As for the pair of distinguishable bosons, a strong fragmentation is possible in the mixture with hardly any shrinking of the density in comparison with that of the bare trap due to the condition $\lambda +{\lambda }_{12}=0$, namely when strong inter-species attractive interaction is accompanied by strong intra-species repulsion of equal magnitude. In summary, inter-species fragmentation in higher-order reduced density matrices is also governed by the ratio $\frac{{\Omega }_{12}}{\omega }$ and takes place both at the attractive and repulsive sectors of interactions.

3.2. Inter-Species Entanglement and the Limit of an Infinite Number of Particles

In the previous sections, the reduced density matrices for identical and distinguishable pairs of bosons were diagonalized, and the intra-species and inter-species fragmentations explored. Both kinds of fragmentations are critical phenomena in the sense that, going to the limit of an infinite number of particles while keeping the interaction parameters (products of the number of particles times the interaction strengths) constant, the respective reduced density matrix per particle becomes $100\%$ condensed [74]. This can be easily found from the leading natural eigenvalues of the natural functions explicitly obtained above; see the general (6), (10), (22), (32) and specific (12), (24), (34) expressions, which are all equal to 1 in this limit.

In the present, concluding subsection, we touch upon a property of the mixture which does not diminish at the limit of an infinite number of particles. Classifying properties of Bose–Einstein condensates and their mixtures at the limit of an infinite number of particles, and especially when many-body and mean-field theories do not coincide, is an active field of research, where variances of observables and the overlap between the many-body and mean-field wave functions are discussed elsewhere; see [75,76,77,93,94,95,96,97,98,99,100,101]. Here, combining the techniques used in the previous sections, we apply Mehler’s formula to perform the Schmidt decomposition of the wave function.

Let us examine the mixture’s wave function, for which the coordinates of the two species are coupled to each other owing to the inter-species interaction; see the last term in (2a). As a reminder, the wave function is obtained by representing the Hamiltonian (1) with the mixture’s Jacoby coordinates, for which it is fully diagonalized, and translating it back to the laboratory frame. To decouple the coordinates of each species, in the sense of prescribing the Schmidt decomposition of the wave function, it is useful to go ‘half a step’ backward, and express (2a), using the individual species’ Jacoby coordinates.

The Jacoby coordinates of each species are given by

$$\begin{array}{cc}\hfill & X_k=\frac{1}{\sqrt{k(k+1)}}{\sum }_{j=1}^k\left(x_{k+1}-x_j\right),1\le k\le N-1,X_N=\frac{1}{\sqrt{N}}{\sum }_{j=1}^N{x}_j,\hfill \\ \hfill & Y_k=\frac{1}{\sqrt{k(k+1)}}{\sum }_{j=1}^k\left(y_{k+1}-y_j\right),1\le k\le N-1,Y_N=\pm \frac{1}{\sqrt{N}}{\sum }_{j=1}^N{y}_j,\hfill \end{array}$$

(38)

where, for the derivation given below, it is useful to distinguish between the two cases for the definition of, say, $Y_N$: The plus sign is assigned to positive $\gamma$, namely, to attractive inter-species interactions for which ${\Omega }_{12}>\omega$, and the minus sign is assigned to negative $\gamma$, i.e., to repulsive inter-species interactions where ${\Omega }_{12}<\omega$.

For the symmetric mixture, given the above Jacobi coordinates of each species, Equation (38), the wave function reads

$$\begin{array}{cc}\hfill & \Psi (X_1,\dots ,X_N,Y_1,\dots ,Y_N)={\left(\frac{m\Omega }{\pi }\right)}^{\frac{N-1}{2}}{\left(\frac{M_{12}{\Omega }_{12}}{\pi }\right)}^{\frac{1}{4}}{\left(\frac{M\omega }{\pi }\right)}^{\frac{1}{4}}\hfill \\ \hfill & \times e^{-\frac{1}{2}m\Omega {\sum }_{k=1}^{N-1}\left(X_k^2+Y_k^2\right)}{e}^{-\frac{\frac{1}{2}m\left({\Omega }_{12}+\omega \right)}{2}\left(X_N^2+Y_N^2\right)}{e}^{\pm \frac{1}{2}m\left({\Omega }_{12}-\omega \right)X_N{Y}_N}.\hfill \end{array}$$

(39a)

Indeed, all relative coordinates are decoupled, and the only coupling due to the inter-species interaction is between the center-of-mass $X_N$ of species 1 bosons and the center-of-mass $Y_N$ of species 2 bosons. Consequently, applying Mehler’s formula to the terms with the intra-species center-of-mass Jacoby coordinates $X_N$ and $Y_N$ the Schmidt decomposition of (39a) is readily performed and given by

$$\begin{array}{cc}\hfill & \Psi (X_1,\dots ,X_N,Y_1,\dots ,Y_N)={\sum }_{n=0}^{\infty }\sqrt{1-{\rho }_{SD}^2}{\rho }_{SD}^n{\Phi }_{1,n}(X_1,\dots ,X_N){\Phi }_{2,n}(Y_1,\dots ,Y_N),\hfill \\ \hfill & {\Phi }_{1,n}(X_1,\dots ,X_N)={\left(\frac{m\Omega }{\pi }\right)}^{\frac{N-1}{4}}{e}^{-\frac{1}{2}m\Omega {\sum }_{k=1}^{N-1}{X}_k^2}\frac{1}{\sqrt{2^{n}n!}}{\left(\frac{s_{SD}}{\pi }\right)}^{\frac{1}{4}}{H}_n\left(\sqrt{s_{SD}}{X}_N\right)e^{-\frac{1}{2}{s}_{SD}{X}_N^2},\hfill \\ \hfill & {\Phi }_{2,n}(Y_1,\dots ,Y_N)={\left(\frac{m\Omega }{\pi }\right)}^{\frac{N-1}{4}}{e}^{-\frac{1}{2}m\Omega {\sum }_{k=1}^{N-1}{Y}_k^2}\frac{1}{\sqrt{2^{n}n!}}{\left(\frac{s_{SD}}{\pi }\right)}^{\frac{1}{4}}{H}_n\left(\sqrt{s_{SD}}{Y}_N\right)e^{-\frac{1}{2}{s}_{SD}{Y}_N^2},\hfill \\ \hfill & \sqrt{1-{\rho }_{SD}^2}=\frac{2\sqrt{\frac{{\Omega }_{12}}{\omega }}}{1+\frac{{\Omega }_{12}}{\omega }},{\rho }_{SD}=\frac{{\left(\frac{{\Omega }_{12}}{\omega }\right)}^{\pm 1}-1}{{\left(\frac{{\Omega }_{12}}{\omega }\right)}^{\pm 1}+1},s_{SD}=m\sqrt{\omega {\Omega }_{12}}.\hfill \end{array}$$

(39b)

We repeat that the plus sign is for attraction and the minus for repulsion, which is what guarantees that ${\rho }_{SD}$ and consequently the Schmidt coefficients $\sqrt{1-{\rho }_{SD}^2}{\rho }_{SD}^n$, $n=0,1,2,3,\dots$ are always positive. $s_{SD}$ defines the inverse width of the individual species’ center-of-mass Gaussians in the Schmidt basis ${\Phi }_{1,n}(X_1,\dots ,X_N)$ and ${\Phi }_{2,n}(Y_1,\dots ,Y_N)$.

Let us concisely discuss the properties of the Schmidt decomposition of the mixture, Equation (39b). Clearly and interestingly, the Schmidt coefficients are independent of the intra-species dressed frequency $\Omega$, which only appears in conjunction with intra-species relative coordinates, i.e., the Schmidt coefficients depend solely on the inter-species interaction. Furthermore, there is a kind of symmetry between respective attractive and repulsive inter-species interactions, as one obtains the same Schmidt coefficients for the inter-species frequency $\frac{{\Omega }_{12}}{\omega }=\sqrt{1+\frac{4N{\lambda }_{12}}{m{\omega }^2}}$ and inverse frequency $\frac{\omega }{{\Omega }_{12}}=\frac{1}{\sqrt{1+\frac{4N{\lambda }_{12}}{m{\omega }^2}}}$.

Last but not least, the same Schmidt coefficients are obtained when the product of the number of bosons in each species times the inter-species interaction strength, $N{\lambda }_{12}$, is held fixed, and N is increased to infinity. In other words, whereas identical and distinguishable bosons, pairs, four-particle aggregates, etc. are $100\%$ condensed at the limit of an infinite number of particles, i.e., the leading eigenvalue of all finite-order intra-species and inter-species reduced density matrices per particle is 1, the mixture’s wave function exhibits a fixed amount of entanglement at the infinite-particle-number limit. This is a good place to bring the present study to an end.

4. Summary and Outlook

The present work aims at developing and combining concepts from quantum theory of many-particle systems with novel results on the physics of trapped mixtures of Bose–Einstein condensates. The notions of natural orbitals and natural geminals are fundamental to many-particle systems made of identical particles. These natural functions entail the diagonalization of the reduced one-particle and two-particle density matrices, respectively. In a mixture of two kinds of identical particles—here explicitly two types of bosons—there are, naturally, identical bosons and pairs made of indistinguishable bosons of either species. To find their natural orbitals and natural geminals, the construction and subsequent diagonalization of respective intra-species reduced density matrices is necessary. In the mixture, there are, additionally, pairs made of distinguishable bosons. Analogously, their theoretical description would require assembling, diagonalizing, and analyzing the inter-species reduced two-particle density matrix. In the present work, we have investigated pairs made of identical or distinguishable bosons in a mixture of Bose–Einstein condensates, covering both the structure of the respective natural pair functions, on the more formal theoretical side, and the exploration of pairs’ fragmentation. Like identical bosons, which can, depending on whether the reduced one-particle density matrix has one or more macroscopic eigenvalues, be condensed or fragmented, so do the pairs of bosons. We showed in the present work that, in the mixture, both pairs made of identical bosons and pairs consisting of distinguishable bosons can be condensed and moreover fragmented.

To tackle the above and other questions, we employed a solvable model, the symmetric harmonic-interaction model for mixtures. The natural geminals for pairs made of identical or distinguishable bosons were explicitly contracted as a function of the inter-species and intra-species interactions. This was performed by diagonalizing the corresponding intra-species and inter-species reduced two-particle density matrices using applications of Mehler’s formula on appropriately constructed linear combinations of intra-species and inter-species coordinates. Here, the role of the mixture’s center-of-mass and relative center-of-mass coordinates was identified and explained. The structure of identical and distinguishable pairs in the mixture was discussed, and a generalization to pairs of distinguishable pairs using the inter-species reduced four-body density matrix was made. A particular case, where attractive and repulsive inter-species and intra-species interactions are (opposite in sign and) equal in magnitude, was worked out explicitly. The fragmentation of bosons, pairs, and pairs of pairs in the mixture was proven, and the size of the respective densities analyzed. Last but not least, as a complementary investigation, the exact Schmidt decomposition of the mixture’s wave function was performed. The entanglement between the two species was shown to be governed by the coupling of their individual center-of-mass coordinates and, consequently, not to vanish at the limit of an infinite number of particles where any finite-order intra-species and inter-species reduced density matrix per particle is 100% condensed.

The results reported in this work were obtained analytically, because the one-body trapping potential and the two-body inter-particle interactions are all harmonic. It is relevant to inquire whether the results would be robust in general. We know that fragmentation of single-species bosons takes place in traps of various shapes; see, e.g., [13], for the inter-particle interactions of various ranges [24], and at the level of the one-particle and two-particle reduced density matrices [28]. Furthermore, as seen in Appendix A, fragmentation of bosons and pairs occurs, respectively, at the level of the one-particle and two-particle reduced density matrices within the single-species harmonic-interaction model. We can therefore quite confidently foresee that the results obtained here analytically, of intra-species and inter-species pairs’ fragmentation using the (symmetric) harmonic-interaction model for mixtures, would persist for trapped interacting bosonic mixtures in general. Demonstrating that explicitly we would have to resort to many-body numerical tools, see, e.g., in [50,59,62,67]. Similarly, our analytical and numerical experience for studying single-species trapped bosons at the limit of an infinite-number of particles [95,96,99,101], along with the corresponding analytical results for variances and overlaps in mixtures [75,77], strongly suggest that a similar generality would hold true for the Schmidt decomposition of trapped interacting mixtures at the infinite-particle-number limit.

The present investigations suggest several directions for further developments. We have treated the symmetric mixture, and an anticipated extension to generic trapped mixtures, with different numbers of bosons, masses, and interaction strengths for each species, would be in place. In what capacity can the fragmentation of identical pairs in the different species be made to differ, and to what extent would the fragmentation of distinguishable pairs become more complex? Do the center-of-mass and relative center-of-mass coordinates keep their role in the diagonalization of inter-species reduced density matrices for a generic mixture? It also makes sense, in a generic mixture, to investigate the fragmentation of aggregates with unequal numbers of bosons from each species, such as, for instance, the analysis of inter-species reduced three-particle density matrices. Another foreseen extension is about mixtures with more species and, if feasible, about generic multi-species mixtures where, e.g., one species could serve as a bridge between two baths. Take for instance the (symmetric) three-species harmonic-interaction model. The diagonalization of the Hamiltonian would require relative coordinates for each of the species, the center-of-mass coordinate of all bosons, and two relative center-of-mass coordinates to be determined [71]. The major challenge would be the reduction of the all-particle density matrix to the various intra-species and inter-species reduced density matrices. Recall that, for single-species bosons, the reduction of the all-particle density matrix leads to ‘vector’ recursive relations (the coefficients of the reduced density matrices depend on one index) [78] and that for a two-species mixture, the reduction of the all-particle density matrix results in ‘matrix’ recursive relations (the coefficients of the various intra-species and inter-species reduced density matrices depend on two indices) [74]. It is thus reasonable to expect that in the three-species mixture, the reduction of the all-particle density matrix would require ‘tensor’ recursive relations. Finally, for a three-species mixture, the question whether an extended or general Schmidt decomposition [102] is at all feasible would become relevant, and if the answer is positive, whether in conjunction, the three individual species’ center-of-mass coordinates would become disentangled. Surely, exciting generalizations are awaiting for further investigations.

Finally, one could forecast that the topic of Bose–Einstein condensates and mixtures in the limit of an infinite number of particles would be enriched by exploring the Schmidt decomposition of the wave function. Recall that at the infinite-particle-number limit, any finite-order intra-species and inter-species reduced density matrix per particle is $100\%$ condensed. Here, studying the variances of observables and overlaps of wave functions has deepened our understanding of the differences between many-body and mean-field theories of Bose–Einstein condensates and mixtures at the limit of an infinite number of particles. However, these properties are already defined for single-species bosons. The Schmidt decomposition, on the other hand, is a property that enters the topic of the infinite-particle-number limit starting, obviously, only from a two-species mixture—all of which paves the way for further intriguing investigations to come.