Fragmentation of identical and distinguishable bosons' pairs and natural geminals of a trapped bosonic mixture

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 fragmentation of pairs in a trapped mixture of Bose-Einstein condensate 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 wavefunction 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.


I. INTRODUCTION
Condensation and fragmentation are basic and widely-studied concepts of Bose-Einstein condensate emanating from the properties of the reduced one-particle density matrix [1][2][3][4][5].
Consider now a mixture of two kinds of identical bosons, which are labeled species 1 and species 2. Mixtures of Bose-Einstein condensate is a highly investigated topic, see, e.g., . One may ask, just like for single-species bosons, about 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 oneparticle 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 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 of the intra-species reduced two-particle density matrices is needed. Summarizing, fragmentation of identical bosons and its manifestation in higherorder reduced density matrices stem from the properties of intra-species quantities.
But, a mixture of Bose-Einstein condensates offers a degree-of-freedom or many-particle construction which do not exist for single-species bosons, namely, inter-species reduced density matrices. Now, if fragmentation of identical bosons and pairs is defined as the macroscopic occupation of respective eigenvalues following the diagonalization of intra-species re-duced density matrices, we may analogously define fragmentation of distinguishable bosons' pairs as macroscopic occupation of the eigenvalues of the inter-species reduced two-particle density matrix. Obviously, the later 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 fragmentation of pairs of identical bosons and establish 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 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 wavefunction and discuss some of its properties at the limit of an infinite-number of particles where the mixture is 100% condensed; and (v) Achieving the first four goals analytically, using an exactly solvable model.
The structure of the paper is as follows. In Sec. II we construct and investigate fragmentation of intra-species and inter-species pair functions in the mixture. In Sec. III we extend the results and explore 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 Sec. III C a summary of the results and an outlook of some prospected research topics are provided. Finally, appendix A collects for comparison with the mixture the details of fragmentation of bosons and pairs in the single-species system.

A. 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]: There are N bosons of type 1, N bosons of type 2, and the mass of each boson is m. λ is the intra-species interaction strength, either between two bosons of type 1 or two bosons of type 2, and λ 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 dimensions. = 1 is used throughout. Employing Jacoby coordinates for the mixture and translating back to the laboratory frame, the 2N-boson wavefunction and corresponding many-particle density matrix are given by with , and the relative center-of-mass M 12 = m 2N and center-of-mass M = 2mN masses. The wavefunction and similarly the many-particle density of the mixture depend on two dressed frequencies, Ω and Ω 12 , and consist of three parts: One-body part with coefficient α, intraspecies two-body coupling with coefficient β, and inter-species two-body coupling with coefficient γ. Whereas α and β depend on the intra-species and inter-species interactions, γ depends on the inter-species interaction only.
Another issue worth mentioning is the stability of the mixture. For a stable, i.e., bound, mixture both dressed frequencies Ω and Ω 12 must be positive. This implies the conditions λ + λ 12 > − mω 2 2N and λ 12 > − mω 2 4N , respectively, on the interactions. In other words, the interspecies interaction λ 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 λ can take any value as long as the inter-species interaction is sufficiently attractive, i.e., λ > − mω 2 2N − λ 12 . We shall return to the dressed frequencies Ω and Ω 12 below.

B. 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, In a symmetric mixture, the corresponding reduced density matrices of species 2, ρ 2 (y, y ′ ) and ρ (2) 2 (y 1 , y 2 , y ′ 1 , y ′ 2 ), 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 lengthly and given in [74]. We start from the final expression for the intra-species reduced one-particle density matrix which is given by 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, ρ 1 (x, x ′ ) depends on the three parts of the many-boson wavefunction, i.e., on the α, β, and γ terms (2a). In the absence of interspecies interaction, i.e., for γ = 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 with s > 0 and, generally, for intra-species and inter-species reduced density matrices as well as later on for Schmidt decomposition of the wavefunction, 1 > ρ ≥ 0. H n are the Hermite polynomials.
Comparing the structure of the intra-species reduced one-particle density matrix ρ 1 (x, x ′ ) with that of Mehler's formula one readily has s (1) Here, 1 − ρ 1 is the condensate fraction of species 1 (and of species 2), i.e., the fraction of condensed bosons, and ρ (1) 1 is the depleted fraction, namely, the fraction of bosons residing outside the lowest, condensed mode. s (1) 1 is the scaling, or effective frequency, of the intraspecies 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 λ and inter-species λ 12 interaction strengths. A specific application of the general expressions (6) for the mixture appears below.
After the transformation (8), the first term of ρ is separable as a function of q 2 and q ′ 2 whereas, using Mehler's formula onto the variables q 1 and q ′ 1 , the second term can be diagonalized. Thus, comparing the second term in (9) and Eq. (5) we find With expressions (10), the decomposition of the intra-species reduced two-particle density matrix in terms of its natural geminals is explicitly given by Equation (11) is a general result on the inter-species natural geminals of the mixture. Together with (10) they imply that 1 − ρ is the fraction of condensed pairs of species 1 (and of species 2), ρ 1 is the fraction of depleted pairs, i.e., the fraction of pairs residing outside the lowest, condensed natural geminal, and s (2) 1 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 λ and inter-species λ 12 interactions. A specific application of the general decomposition (10), (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 follows the above pattern and will not be discussed further here.
Let us work out an explicit application where we shall find and analyze fragmentation of identical bosons' pairs. Consider the specific scenario where λ + λ 12 = 0, i.e., that the intraspecies 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, Ω = mω. Then, the coefficients of the three parts of the wavefunction simplify and one has α = mω + β = mω 1 + 1 2N Ω 12 ω − 1 and β = γ = m 2N (Ω 12 − ω). Consequently, 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 s (1) for the intra-species reduced one-particle density matrix, where α + C 1,0 = mω is used, and for the intra-species reduced two-particle density matrix, where α+β = mω 1 + 1 N Ω 12 are utilized. We see that fragmentation of identical pairs and bosons is governed by the ratio Ω 12 ω and its inverse ω Ω 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 Ω 12 ω and ω Ω 12 , see discussion below.
Let us analyze explicitly macroscopic fragmentation of geminals, i.e., when there is macroscopic occupation of more than a single intra-species natural pair function of . As a reference, we also refer to the corresponding and standardly defined macroscopic fragmentation of the intra-species natural orbitals of ρ (1) 1 (x, x ′ ). 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 ρ = 1 − ρ = 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 − ρ)ρ n , 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 Ω 12 ω , or, equivalently, for which inter-species interaction λ 12 = mω 2 4N Ω 12 ω 2 − 1 , 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 and working out (12b) for 50% natural-geminal fragmentation we get 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 Ω 12 ω is large, or near the border of stability for repulsions, that is when Ω 12 ω is close to zero. Also, 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 natural-geminal with 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 ρ 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 1, 1, 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, lim Ω 12 1,x −→ ∞ for 1, 1, where σ 1, 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 λ + λ 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 fragmentation of single-species bosons in the harmonic-interaction model, see appendix A.
C. 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 as For the harmonic-interaction model of the symmetric mixture it can be computed analytically and, starting from (2b), is given by [74] ρ (2) 12 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 ρ 12 (x, x ′ , y, y ′ ) one must couple and make linear combinations of coordinates associated with distinguishable bosons. Defining u = 1 (19): Consequently, we readily find the decomposition where the normalizations after and before diagonalization are, of course, equal. As might be expected, since the structure of ρ 12 (x, x ′ , y, y ′ ) is more involved than that of ρ 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 ′ and v, v ′ . We thus get where the "+" terms quantify the fragmentation in the u, u ′ part of the intra-species reduced two-particle density matrix and the "−" terms quantify the fragmentation in the v, v ′ part of the intra-species reduced two-particle density matrix, also see below. Equation (22)  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 12,n + ,n − (x, y) = All in all, (23) implies that the distinguishable-pair 'condensed fraction' is given by 12,− and the respective depleted fraction by 1− 1 − ρ 12,− . Each of the inter-species 'mixed' coordinates x±y √ 2 carry its own scaling, s 12,± . The distinguishable natural geminals Φ (2) 12,n + ,n − (x, y) are, needless to say, orthonormal to each other.
Furthermore, there are different numbers of pairs: 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 Ω 12 ω = 1 + 4N mω 2 λ 12 for which the inter-species reduced two-particle density matrix is 50% fragmented as in (13). Since ρ (2) 12,+ = 0 does not contribute, the only contribution to fragmentation comes from ρ =⇒ 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 λ + λ 12 = 0 in our example, the respective intra-species interaction is opposite in sign. Also, 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 ρ 12 (x, y) = ρ 12 (x, x ′ = x, y, y ′ = y), which is given by 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 .
To show the combined effect of the intra-species and inter-species interactions accompanying fragmentation of ρ 12 (x, x ′ , y, y ′ ), 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, lim Ω 12 where σ (2) 12, x+y √ 2 is independent of the interactions. We see that the size of the inter-species density saturates as well at about 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 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 Ω 12 ω 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 twoparticle densities (26) and (16).

III. PAIR OF DISTINGUISHABLE PAIRS AND SCHMIDT DECOMPOSITION OF THE WAVEFUNCTION
Following the results of the previous section on 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? In as much as single-species and intra-species fragmentations take place at the lowest-level reduced oneparticle 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 for 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 done by constructing the Schmidt decomposition of the mixture's wavefunction and discussing consequences of this decomposition at the limit of an infinite number of particles.
A. Inter-species fragmentation in higher-order reduced density matrices The inter-species reduced four-particle density matrix is defined as 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 ρ (4) where and α+C 1,1 ∓D 1,1 are given in (19b). The combinations of parameters α+β +2 (C 2,2 ∓ D 2,2 ) would appear below shortly.
To diagonalize ρ 12 (x 1 , x 2 , x ′ 1 , x ′ 2 , y 1 , y 2 , y ′ 1 , y ′ 2 ) 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, , we have for the different terms in (29a): Relations (30) imply that one could equally define inter-species linear combinations of the relative coordinates, since u 2 2 + v 2 2 = (x 1 −x 2 )+(y 1 −y 2 ) . We chose the former combinations.
Plugging (30) into (29) we readily find for the transformed inter-species reduced fourparticle density matrix ρ (4) where the normalization coefficients before and after diagonalization are, naturally, equal As can be seen in (31) and (21), the similarities and differences between the structures of ρ 12 (u, u ′ , v, v ′ ) clarify the issue of which coordinates are coupled and identify the coordinates that are not. In particular, just like for the twoparticle quantity, we can apply Mehler's formula twice, on the appropriately-constructed inter-species 'mixed coordinates' u 1 , u ′ 1 and v 1 , v ′ 1 , to diagonalize the inter-species reduced four-particle density matrix. When this is done one obtains s (4) where the "+" terms quantify the fragmentation in the u 1 , u ′ 1 part of the inter-species reduced four-particle density matrix and the "−" terms determine the fragmentation in the v 1 , v ′ 1 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 and v 2 , v ′ 2 . 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 will not be pursued further here.
We can now prescribe the decomposition of the intra-species reduced four-particle density matrix, in as much as the reduced two-particle density matrix was decomposed, into its natural four-particle functions made of distinguishable particles which is given by 12,n + ,n − (x 1 , x 2 , y 1 , y 2 ) = Equation (33) means that the pair-of-distinguishable-pairs 'condensed fraction' is given by the product 1 − ρ We proceed now to examine fragmentation in this higher-order inter-species reduced density matrix. We investigate as above the specific case of λ + λ 12 = 0. To compute we require the quantities [α + β + 2 (C 2,2 − D 2,2 )] = mω and α + β − 2D ′ 2,2 = mω for the "+" branch as well as [α + β + 2 (C 2,2 + D 2,2 )] = 12, 12,+ = 1, indicating that there is no contribution to fragmentation from the center-of-mass 'mixed coordinate' u 1 , u ′ 1 . This is additional to the no contribution to fragmentation coming from the relative coordinates u 2 , u ′ 2 and v 2 , v ′ 2 , see (31). For the relative center-of-mass 'mixed coordinate' v 1 , v ′ 1 , on the other end, one finds s (4) 12,− = mω namely, that the fragmentation of pairs of distinguishable pairs fully originates from the relative center-of-mass 'mixed coordinate' v 1 , v ′ 1 . We see that also the higher-order interspecies fragmentation is governed by the ratio Ω 12 ω and takes place both at the attractive and repulsive sectors of interactions. Concluding, higher-order inter-species fragmentation is proved. Now, one can compute the ratio Ω 12 ω = 1 + 4N mω 2 λ 12 for which the inter-species reduced four-particle density matrix is 50% fragmented as in (13). Since ρ (4) 12,+ = 0 does not contribute in this specific case, the only contribution to fragmentation comes from ρ (4) 12,− . Thus, solving (34b) for 50% distinguishable-four-particle-function fragmentation we obtain ρ (4) As for distinguishable pairs, there are two 'reciprocate' solutions, one for strong attractions and the second close to the border of stability for repulsions. Also, to achieve the same degree of 50% with a larger number 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 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 ρ 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 σ 12, , σ 12, , σ 12, (x 1 +x 2 )+(y 1 +y 2 ) 2 12, (x 1 +x 2 )−(y 1 +y 2 ) 2 To show the combined effect of the inter-species and intra-species interactions accompanying fragmentation of ρ 12 (x 1 , x 2 , x ′ 1 , x ′ 2 , y 1 , y 2 , y ′ 1 , y ′ 2 ), 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, 12, where σ 12, 12, , and σ 12, (x 1 +x 2 )+(y 1 +y 2 ) 2 do not depend on the interactions. We see that the size of the inter-species four-boson density saturates as well 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 hardy any shrinking of the density in comparison with the bare trap due to the condition λ + λ 12 = 0, i.e., that strong inter-species attractive interaction is accompanied by strong intra-species repulsion of equal magnitude. Summarizing, interspecies fragmentation in higher-order reduced density matrices is also governed by the ratio ω and takes place both at the attractive and repulsive sectors of interactions.

B. 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 obtained from the leading natural eigenvalues of the natural functions explicitly obtained above, see the general (6) The Jacoby coordinates of each species are given by 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 γ, namely, to attractive inter-species interactions for which Ω 12 > ω, and the minus sign is assigned to negative γ, i.e., to repulsive inter-species interactions where Ω 12 < ω.
For the symmetric mixture, given the above Jacobi coordinates of each species, Eq. (38), the wavefunction reads Mω π 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 intra-species center-ofmass Jacoby coordinates X N and Y N the Schmidt decomposition of (39a) is readily performed and given by We remind that the plus sign is for attraction and the minus for repulsion, which is what guarantees that ρ SD and consequently the Schmidt coefficients 1 − ρ 2 SD ρ n SD , n = 0, 1, 2, 3, . . . are always positive. s SD defines the inverse width of the individual species' center-of-mass Gaussians in the Schmidt basis Φ 1,n (X 1 , . . . , X N ) and Φ 2,n (Y 1 , . . . , Y N ).
Let us concisely discuss properties of the Schmidt decomposition of the mixture, Eq. (39b).
Clearly and interestingly, the Schmidt coefficients are independent of the intra-species dressed frequency Ω, 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 interspecies interactions, as one gets the same Schmidt coefficients for the inter-species frequency Ω 12 ω = 1 + 4N λ 12 mω 2 and inverse frequency ω 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λ 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 interspecies reduced density matrices per particle is 1, the mixture's wavefunction 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.

C. 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 manyparticle 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 in need. 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 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 more so 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 identi-cal or distinguishable bosons were explicitly contracted as a function of the inter-species and intra-species interactions. This was done by diagonalizing the corresponding intraspecies 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 magnitude, has been worked out explicitly. Fragmentation of bosons, pairs, and pairs of pairs in the mixture has been 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 wavefunction 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 present investigations suggest several directions for further developments. We have It also makes sense, in a generic mixture, to investigate fragmentation of aggregates with unequal numbers of bosons from each species, like, for instance, the analysis of inter-species reduced three-particle density matrices. Another extension foreseen is to mixtures with more species and, if feasible, to generic multi-species mixtures where, e.g., one species could serve as a bridge between two baths. 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 wavefunction. Recall that at the infiniteparticle-number limit any finite-order intra-species and inter-species reduced density matrix per particle is 100% condensed. Here, observables' variances and wavefunctions' overlaps have deepened our understanding of the differences between many-body and mean-field theories of Bose-Einstein condensates and mixtures at limit of an infinite number of particles, but are properties 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. The Hamiltonian of the single-species harmonic-interaction model is given by [78] H(x 1 , . . . ,

ACKNOWLEDGEMENTS
Employing single-species Jacoby coordinates and translating back to the laboratory frame, the N-boson wavefunction and corresponding density matrix are given by The stability of the system means that the interaction satisfies λ > − mω 2 2N . The reduced one-particle density matrix reads (A3) Comparing the structure of reduced single-particle density matrix with that of Mehler's formula [81,82] one readily has (A4) The reduced two-particle density matrix ρ (2) (x 1 , x 2 , x ′ 1 , x ′ 2 ) reads, after rotation of coordinates, where . Comparing the structure of the reduced two-particle density matrix with that of Mehler's formula (5) we readily find s (2) = (α + β) (α + β + 2C 2 ) = mΩ 1 + 2 where α + β = mΩ 1 + 2 N ω Ω − 1 . We see from (A4) and (A6) that fragmentation of bosons and pairs is governed, in the single-species harmonic-interaction model, by the ratio Ω ω and takes place both at the attractive and repulsive sectors of the interaction.
Similarly to the main text, we compute for which ratio Ω ω , or, equivalently, for which interaction λ = mω 2

2N
Ω ω 2 − 1 , the two-particle and one-particle reduced density matrices are macroscopically fragmented as in (13). Note the difference that, here, Ω ω is the singlespecies interaction and in the main text, for the mixtue, Ω 12 ω is the inter-species interaction. Thus, solving (A4) for 50% natural-orbital fragmentation we find and working out (A6) for 50% natural-geminal fragmentation we get There are two 'reciprocate' solutions for both the natural geminals and natural orbitals: Indeed, 50% fragmentation occurs for strong attractions, namely, when Ω ω is large, or in the vicinity of the border of stability for repulsions, i.e., when Ω ω is close to zero. Also, to achieve the same degree of 50% with a larger number N of bosons, a stronger attraction or repulsion is needed. Finally, comparing natural-geminal with natural-orbital fragmentation at the same 50% value, one sees from (A8) and (A7) that slightly weaker interactions, attractions or repulsions, are needed for the former, in a similar manner to intra-species fragmentation in the mixture discussed in the main text.
Finally, we prescribe the one-particle and two-particle densities, i.e., the diagonal parts ρ (1) (x) = ρ (1) (x, x ′ = x) and ρ (2)  (A9) Here, in the single-species case, λ governs the size of the densities which, as we now show, depend on the interaction stronger than for the mixture, also see the main text for comparison and further discussion.
We can deduce the size of pairs' and bosons' clouds using the widths of the corresponding Gaussian functions in the densities (A9). Hence, we have x 1 +x 2 √ 2 , σ x 1 −x 2 √ 2 = 1 2mΩ .
To describe further effects of the interaction λ accompanying fragmentation of the reduced density matrices (A4) and (A6), we compute the sizes (A10a) for strong attractions or repulsions at the border of stability. We obtain, respectively, as all widths (A10a) depend on the interaction strength. The size of the densities for strong attractions diminishes to much smaller values than the trap's size that do not depend on the strength of interaction. Thus, a high degree of fragmentation due to the strong attractive interaction is possible in the single-species system only together with shrinking of the density.
Alternatively, towards the edge of stability the density of the single-species system expands due to repulsive interaction unlimitedly as the degree of fragmentation increases. These should be compared to and contrasted with the results in the main text for the intra-species densities and the interplay between inter-species and intra-species interactions within intraspecies fragmentation in the mixture.