Quantum correlations and spatial localization in one-dimensional ultracold bosonic mixtures

The fascinating physics of interpenetrating superfluids has recently become a topic of great interest due to the experimental realization of multi-component, atomic Bose–Einstein condensates [1–5]. In the weakly interacting regime, these mixtures are well described by coupled mean-field Gross–Pitaevskii equations (GPEs), and within this framework, processes that lead to phase separation are well described [6–14].

While mean-field theories allow us to study weakly correlated systems, it is also important and interesting to examine quantum mixtures in strongly correlated regimes. In these regimes, analytic solutions can often be obtained only in limiting cases. Rather appealing results occur in strongly correlated regimes when the dimensionality is reduced. For quasi-one-dimensional (1D) gas mixtures, one finds that the Luttinger liquid theory predicts many interesting effects, which include de-mixing for repulsive interactions or spin-charge separation analogous to that found in 1D electronic quantum systems [15–18]. Other relevant effects include the presence of polarized ground states, which allow us to view the relative spatial oscillations as spin waves [19–23] and which have been experimentally observed [24–26].

Very strong correlations for single-component bosons are realized in the Tonks–Girardeau (TG) gas [27–29], which was recently observed experimentally [34, 35]. Bosonic mixtures in the strongly interacting limit have features common with the TG gas, and their ground-state wavefunction can be obtained analytically in certain interaction limits [36–38]. Experimental advances on Feshbach and confinement-induced resonances in recent years have made it possible to control both the intra-species interactions and the inter-species interactions over a wide range of parameters [39–41]. In the strongly interacting limit, a number of relevant phenomena have been described, including phase separation [15–17, 42], composite fermionization [43–45], a sharp crossover between both limits [46], and quantum magnetism [47].

In this work we focus on mixtures where the number of atoms is small. The recent successful experimental trapping of ensembles of a few atoms [48–52] has inspired an intense theoretical effort [53–67], and very recently even systems with SU(N) symmetry and $N\gt 2$ have been experimentally realized [68]. For mixtures of a few atoms, direct diagonalization methods [31, 42, 45] can be used together with other numerical methods efficient for larger numbers of atoms, like multiconfigurational Hartree–Fock methods (MCTDH) [70], density functional theory (DFT) [44], or quantum diffusion Monte Carlo (DMC) [69]. In the present work, we use direct numerical diagonalization to study the ground-state properties of a mixture of ultracold bosons confined in a 1D trap over a wide range of correlations regimes, determined by the scattering properties between the atoms. These are supplemented by DMC calculations to confirm trends for systems with larger particle numbers. While the extreme cases in which all correlations are either weak or strong are well known, here we calculate and discuss the full phase diagram and especially the transitions between the different regimes.

We study the ground-state wavefunction and pay particular attention to the one- and two-body correlations in the extreme limits and across the transitions between them. The quantum correlations between both components are characterized by means of the von Neumann entropy. This allows us to show that close to the crossover between the composite fermionization and phase separation, the ground state exhibits strong correlations between the two bosonic components.

Our paper is organized as follows. In section 2 we introduce the model Hamiltonian and a general analytical ansatz for the ground-state wavefunction. Focusing first on balanced mixtures, we discuss in section 3 the ground-state properties in terms of the densities, the coherence, the energies, the one- and two-body correlations, and the von Neumann entropy. In section 4, we then present results on how the ground-state properties change when one component is larger than the other, and finally summarize all our results in section 5.

Let us consider a mixture of two bosonic components, A and B, with a small, fixed number of atoms in each component, ${{N}_{{\rm A}}}$ and ${{N}_{{\rm B}}}$. We assume that the two components are two different hyperfine states of the same atomic species of mass m and that they are trapped in the same, 1D parabolic potential $V(x)=\frac{1}{2}m{{\omega }^{2}}{{x}^{2}}$. At low temperatures, all scattering processes between the atoms are assumed to be described by contact interactions $v_{{\rm int}}^{{\rm A}}={{g}_{{\rm A}}}\delta ({{x}_{j}}-{{x}_{{{j}^{\prime }}}})$, $v_{{\rm int}}^{{\rm B}}={{g}_{{\rm B}}}\delta ({{y}_{j}}-{{y}_{{{j}^{\prime }}}})$, and $v_{{\rm int}}^{{\rm AB}}={{g}_{{\rm AB}}}\delta ({{x}_{i}}-{{y}_{j}})$, where the positions of atoms of species A(B) are given by the coordinates ${{x}_{j}}({{y}_{j}})$. The 1D intra- and inter-species coupling constants ${{g}_{{\rm A}({\rm B})}}$ and ${{g}_{{\rm AB}}}$ are assumed to be tunable independently by means of confinement-induced resonances [39]. We will restrict our study to repulsive interactions. The many-body Hamiltonian is $\hat{\mathcal{H}}={{\hat{\mathcal{H}}}_{{\rm A}}}+{{\hat{\mathcal{H}}}_{{\rm B}}}+{{\hat{\mathcal{H}}}_{{\rm int}}}$, with:

$$\begin{eqnarray}\begin{array}{rcl} {{{\hat{\mathcal{H}}}}_{{\rm A}}} & = & \mathop{\sum }\limits_{j=1}^{{{N}_{{\rm A}}}}\left[ -\frac{{{\hbar }^{2}}}{2\;m}\frac{{{\partial }^{2}}}{\partial x_{j}^{2}}+V({{x}_{j}}) \right]+\mathop{\sum }\limits_{j\lt j^{\prime} }^{{{N}_{{\rm A}}}}v_{{\rm int}}^{{\rm A}}({{x}_{j}},{{x}_{{{j}^{\prime }}}}), \\ {{{\hat{\mathcal{H}}}}_{{\rm B}}} & = & \mathop{\sum }\limits_{j=1}^{{{N}_{{\rm B}}}}\left[ -\frac{{{\hbar }^{2}}}{2\;m}\frac{{{\partial }^{2}}}{\partial y_{j}^{2}}+V({{y}_{j}}) \right]+\mathop{\sum }\limits_{j\lt j^{\prime} }^{{{N}_{{\rm B}}}}v_{{\rm int}}^{{\rm B}}({{y}_{j}},{{y}_{{{j}^{\prime }}}}), \\ {{{\hat{\mathcal{H}}}}_{{\rm int}}} & = & \mathop{\sum }\limits_{j=1}^{{{N}_{{\rm A}}}}\mathop{\sum }\limits_{j^{\prime} =1}^{{{N}_{{\rm B}}}}v_{{\rm int}}^{{\rm AB}}({{x}_{j}}-{{y}_{{{j}^{\prime }}}}). \\ \end{array}\end{eqnarray} \tag{ 1 }$$

There are three coupling constants ${{g}_{{\rm A}}},{{g}_{{\rm B}}},{{g}_{{\rm AB}}}$, each of them ranging from g = 0 for the ideal Bose gas interaction to $g\to \infty$ for a strong Tonks–Girardeau interaction. This defines eight limits schematically shown in figure 1. The composite fermionization limit is reached when ${{g}_{{\rm AB}}}\to \infty$, with the other coupling constants vanishing [43–45]. We termed TG–BEC gas a system with one of the intra-species coupling constants large, while other coupling constants vanish [42]. If one of the intra-species coupling constants, together with the inter-species coupling constant, are large, the phase separation limit is reached [15–17, 42]. Finally, if all coupling constants tend to infinity, the wavefunction is known exactly and can be mapped to that of an ideal Fermi gas [36]. We call this limit full fermionization. Note that for a system of two A atoms and one B atom, the order in which the limits are taken (i.e., if one first takes ${{g}_{{\rm A}}}\to \infty$ and then ${{g}_{{\rm AB}}}\to \infty$ or v.v.) determines the particular ground state that is reached [38]. In the following we will calculate and discuss the complete phase diagram, which includes the transitions between these limits. To restrict the large number of free parameters, we note that the transition between TG and the phase separation limit is symmetric when switching the values of ${{g}_{{\rm A}}}$ and ${{g}_{{\rm B}}}$, and we can therefore circumscribe the discussion to the situation where ${{g}_{{\rm B}}}$ vanishes, and change ${{g}_{{\rm A}}}$. In the following, we will use harmonic oscillator units and scale all lengths in units of oscillator length ${{a}_{0}}=\sqrt{\hbar /(m\omega )}$ and all energies in units of level spacing $\hbar \omega$.

Zoom In Zoom Out Reset image size

To solve the Hamiltonian (1), we use two different numerical approaches: direct diagonalization [42] and DMC [69]. The former allows us to calculate the full density matrix of the system and therefore gives us access to all single- and multi-particle correlations. However, since it is limited to small particle numbers, the latter will be used to check for trends when the number of particles becomes larger. While DMC is well described in the literature, let us briefly explain our approach to direct diagonalization. For this we expand the second quantized field operators into eigenfunctions ${{\phi }_{n}}(x)$ of the single-particle (SP) Hamiltonian for the harmonic oscillator:

$$\begin{eqnarray}&&{{\hat{\psi }}_{{\rm A}}}(x)=\mathop{\sum }\limits_{n=1}^{{{n}_{{\rm A}}}}{{\hat{a}}_{n}}{{\phi }_{n}}(x),\qquad {\rm and}\qquad {{\hat{\psi }}_{{\rm B}}}(x)=\mathop{\sum }\limits_{n=1}^{{{n}_{{\rm B}}}}{{\hat{b}}_{n}}{{\phi }_{n}}(x),\end{eqnarray} \tag{ 2 }$$

where the creation and annihilation operators $\hat{a}_{k}^{\dagger }$ and ${{\hat{a}}_{k}}$ satisfy the bosonic commutation relations $[{{\hat{a}}_{k}},\hat{a}_{l}^{\dagger }]={{\delta }_{kl}}$, $[{{\hat{a}}_{k}},{{\hat{a}}_{l}}]=[\hat{a}_{k}^{\dagger },\hat{a}_{l}^{\dagger }]=0$ and similarly for $\hat{b}_{k}^{\dagger }$ and ${{\hat{b}}_{k}}$, while all commutators between operators belonging to different species vanish. Here, ${{n}_{{\rm A}({\rm B})}}$ is the number of modes used in the expansion. The Hamiltonian can then be written as [46]:

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm A}}}=\mathop{\sum }\limits_{k}\hat{a}_{k}^{\dagger }{{\hat{a}}_{k}}\hbar \omega \left( \frac{1}{2}+k \right)+\frac{1}{2}\mathop{\sum }\limits_{klmn}\hat{a}_{k}^{\dagger }\hat{a}_{l}^{\dagger }{{\hat{a}}_{m}}{{\hat{a}}_{n}}V_{klmn}^{{\rm A}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm B}}}=\mathop{\sum }\limits_{k}\hat{b}_{k}^{\dagger }{{\hat{b}}_{k}}\hbar \omega \left( \frac{1}{2}+k \right)+\frac{1}{2}\mathop{\sum }\limits_{klmn}\hat{b}_{k}^{\dagger }\hat{b}_{l}^{\dagger }{{\hat{b}}_{m}}{{\hat{b}}_{n}}V_{klmn}^{{\rm B}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm int}}}=\mathop{\sum }\limits_{klmn}\hat{a}_{k}^{\dagger }\hat{b}_{l}^{\dagger }{{\hat{b}}_{m}}{{\hat{a}}_{n}}V_{klmn}^{{\rm AB}}\;,\end{eqnarray} \tag{ 5 }$$

where:

$$\begin{eqnarray}&&V_{klmn}^{{\rm A}({\rm B})}={{g}_{{\rm A}({\rm B})}}\int {\rm d}x\ \phi _{k}^{*}(x)\phi _{l}^{*}(x){{\phi }_{m}}(x){{\phi }_{n}}(x),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&V_{klmn}^{{\rm AB}}={{g}_{{\rm AB}}}\int {\rm d}x\ \phi _{k}^{*}(x)\phi _{l}^{*}(x){{\phi }_{m}}(x){{\phi }_{n}}(x).\end{eqnarray} \tag{ 7 }$$

The ground state can be expressed in terms of Fock vectors ${{\Psi }_{0}}=\sum _{i=1}^{\Omega }{{c}_{i}}{{\Phi }_{i}}$ with:

$$\begin{eqnarray}&&{{\Phi }_{i}}=D_{i}^{{\rm A}}D_{i}^{{\rm B}}{{\left( \hat{a}_{1}^{\dagger } \right)}^{N_{1,i}^{{\rm A}}}}...{{\left( \hat{a}_{{{n}_{{\rm A}}}}^{\dagger } \right)}^{N_{{{n}_{{\rm A}}},i}^{{\rm A}}}}{{\left( \hat{b}_{1}^{\dagger } \right)}^{N_{1,i}^{{\rm B}}}}...{{\left( \hat{b}_{{{n}_{{\rm B}}}}^{\dagger } \right)}^{N_{{{n}_{{\rm B}}},i}^{{\rm B}}}}{{\Phi }_{0}},\end{eqnarray} \tag{ 8 }$$

where $D_{i}^{{\rm A}({\rm B})}={{\left( N_{1,i}^{{\rm A}({\rm B})}!...N_{{{n}_{{\rm A}({\rm B})}},i}^{{\rm A}({\rm B})}! \right)}^{-\frac{1}{2}}}$, and Φ₀ is the vacuum. The occupation numbers of the ${{n}_{{\rm A}}}$ (n_B) modes for each component are given by $N_{1,i}^{{\rm A}},...,N_{{{n}_{{\rm A}}},i}^{{\rm A}}$ ($N_{1,i}^{{\rm B}},...,N_{{{n}_{{\rm B}}},i}^{{\rm B}}$). The dimension of the Hilbert space is $\Omega ={{\Omega }_{{\rm A}}}{{\Omega }_{{\rm B}}}$ with ${{\Omega }_{{\rm A}({\rm B})}}=({{N}_{{\rm A}({\rm B})}}+{{n}_{{\rm A}({\rm B})}}-1)!/[{{N}_{{\rm A}({\rm B})}}!({{n}_{{\rm A}({\rm B})}}-1)!]$. Note that $\Omega$ increases exponentially with the number of particles and modes, which is the reason why the numerical solution using this approach is restricted to a small number of atoms.

A good ansatz for the unnormalized ground-state wavefunction of the mixture, when ${{g}_{{\rm B}}}=0$ and outside of the phase-separated regime, can be constructed using the solution for non-interacting atoms in the harmonic trap $\Phi (X)={\rm exp} [-\sum x_{i}^{2}/2]$, $X=\{{{x}_{i}}\}$ and $Y=\{{{y}_{i}}\}$, as in [46]:

$$\begin{eqnarray}&&\Psi (X,Y)=\Phi (X)\;\Phi (Y)\mathop{\prod }\limits_{j\lt k}^{{{N}_{{\rm A}}}}|{{x}_{k}}-{{x}_{j}}-{{a}_{{\rm A}}}|\quad \mathop{\prod }\limits_{k}^{{{N}_{{\rm A}}}}\mathop{\prod }\limits_{j}^{{{N}_{{\rm B}}}}|{{x}_{k}}-{{y}_{j}}-{{a}_{{\rm AB}}}|.\end{eqnarray} \tag{ 9 }$$

Here the 1D s-wave scattering length ${{a}_{{\rm A}}}$ for the intra-species interactions and ${{a}_{{\rm AB}}}$ for the inter-species interactions are related to the 1D coupling constants as ${{g}_{{\rm A}}}=-2{{\hbar }^{2}}/(m{{a}_{{\rm A}}})$ and ${{g}_{{\rm AB}}}=-2{{\hbar }^{2}}/(m{{a}_{{\rm AB}}})$, and we assume that both coupling constants are non-negative, corresponding to repulsive interactions. For practical purposes, we find that the coupling constant g = 20 is close enough to the infinite limit, and therefore we use this value in the direct diagonalization method for describing the large coupling constant limits.

In the following, we will first concentrate on systems in which both components have the same particle number. Unless otherwise stated, we will use ${{N}_{{\rm A}}}={{N}_{{\rm B}}}=2$.

3.1. Densities

The main feature of the density evolution in this system is the occurrence of phase separation for increasing inter-species interactions. However, this process takes two, fundamentally different, forms: in the composite fermionization limit, atoms of different species avoid each other, even though the species' densities still occupy the same space, whereas in the phase separation limit the overlap of the respective densities goes to zero.

The density along the transition from the BEC–BEC limit (all couplings small) to the composite fermionization limit (${{g}_{{\rm AB}}}$ large) is shown in figures 2(a) and (b). There are crucial differences in the evolution of the density along the transition from the TG–BEC to the phase separation limit (figures 2(c) and (d)). One immediately notices that the transition into the composite fermionization state happens at a finite value of ${{g}_{{\rm AB}}}\sim 2$, whereas the transition to the phase-separated regime happens already for very small values of ${{g}_{{\rm AB}}}$. Also, the final states reached in the composite fermionization or the phase separation limit are very different.

Zoom In Zoom Out Reset image size

This difference in the final states can be understood by looking at the one-body density matrix (OBDM) given by:

$$\begin{eqnarray}&&\rho _{1}^{{\rm A}}(x,x^{\prime} )={{N}_{{\rm A}}}\int {\rm d}{{x}_{2}}\;\cdots \;{\rm d}{{x}_{{{N}_{{\rm A}}}}}\;{\rm d}{{y}_{1}}\;\cdots \;{\rm d}{{y}_{{{N}_{{\rm B}}}}}|\Psi {{|}^{2}}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&={{\sum }_{k}}{{f}_{k}}(x){{f}_{k}}(x^{\prime} )\lambda _{{\rm A}}^{k}\end{eqnarray} \tag{ 11 }$$

with a similar expression for $\rho _{1}^{{\rm B}}(x,x^{\prime} )$. The decomposition in terms of natural orbitals f_k(x) of the OBDM and their corresponding occupations $\lambda _{{\rm A}}^{k}$ is given in equation (11). The densities shown in figure 2 are the diagonals of these matrices, calculated with direct diagonalization. As discussed in reference [46], the OBDMs of both components in the composite fermionization limit are identical and show two peaks. Conversely, in the phase separation limit, the OBDM of B shows a single peak located at the center of the trap, while the OBDM for A shows two peaks at the edges. The largest used value of the coupling constant g = 8 is large enough so that the density profiles shown in figure 2 are practically the same as for the infinite coupling constant.

Finally, the transition from the composite fermionization to the phase-separated regime is shown in figures 2(e) and (f). One can see that the spatial separation of the clouds happens for a finite value of ${{g}_{{\rm A}}}$. At the transition between both limits, the OBDMs of both species show a complicated structure, which we discuss in detail in subsection 3.4.

3.2. Coherence and entanglement

Since increasing the coupling constant will drive the system from the weakly to the strongly correlated regime, the coherence is a good quantity for identifying different regions in the phase diagram. It can be characterised by the largest eigenvalue of the OBDMs (11), $\lambda _{0}^{{\rm A}({\rm B})}$, which provides the largest occupation of a natural orbital. In our numerical calculations with direct diagonalization, we normalize the OBDM to 1 instead of the number of atoms. In figures 3(a) and (b) we show the largest occupation numbers for the A and the B species, respectively, over the whole range of interactions. Note that all eigenvalues of each component sum up to 1, in accordance with the chosen normalization.

Zoom In Zoom Out Reset image size

One can see from figure 3(a) that the coherence in the A species decreases monotonically along the transition from the BEC–BEC ($\lambda _{0}^{{\rm A}}=1$) to the TG–BEC ($\lambda _{0}^{{\rm A}}\sim 0.7$) limit, as well as to the composite fermionization limit ($\lambda _{0}^{{\rm A}}\sim 0.55$). However, the transition for increasing ${{g}_{{\rm A}}}$ at a finite ${{g}_{{\rm AB}}}$ shows that a maximum of coherence is reached for finite values of ${{g}_{{\rm A}}}\sim 5$, which corresponds roughly to the value where the cloud de-mixing happens (see figures 2(e) and (f)). This maximum in coherence within species A is very surprising, as usually the presence of interactions is thought of as detrimental to coherence. Here, however the presence of interactions within the A component to a certain degree 'counterbalance' the interactions between the species and therefore allows us to re-establish a higher degree of coherence again. Note that after the de-mixing transition, the coherence within species A goes down again, which is a clear indication that the enhancement is somehow mediated using the overlap with species B.

As expected, species B shows a large degree of coherence in all limits, except the composite fermionization one (see figure 3(b)). However, the re-establishment of coherence along the transition from composite fermionization to the phase-separated limit happens over a definite and narrow region, which corresponds to the area in which the coherence in species A shows a maximum.

One might at this point wonder how the transition to phase separation manifests itself during the transition from the TG–BEC to the phase-separated limit, as no obvious signature is visible in the coherence phase diagram. The answer is that phase separation happens already for small values of ${{g}_{{\rm AB}}}$, which can be seen in figure 2(c).

It is important to observe that there are no phase transitions in the whole phase diagram. The ground-state energy is always a continuous and smooth function of the parameters so that the transition between the different regimes is of a crossover type.

Closely related to the coherence in the sample is the entanglement between the two components. This can be quantified by calculating the von Neumann entropy, ${{S}_{{\rm A}}}=-{\rm Tr}\left( {{\rho }_{{\rm A}}}{\rm ln} {{\rho }_{{\rm A}}} \right)$, which is a function of the reduced density matrix for a single component:

$$\begin{eqnarray}&&{{\rho }_{{\rm A}}}={\rm T}{{{\rm r}}_{{\rm B}}}\rho =\mathop{\sum }\limits_{i}\langle \Phi _{i}^{{\rm B}}|{{\Psi }_{0}}\rangle \langle {{\Psi }_{0}}|\Phi _{i}^{{\rm B}}\rangle .\end{eqnarray} \tag{ 12 }$$

Here $\rho =|{{\Psi }_{0}}\rangle \langle {{\Psi }_{0}}|$ is the density matrix, Ψ₀ is the system ground state, and:

$$\begin{eqnarray}&&\Phi _{i}^{{\rm B}}=D_{i}^{{\rm B}}{{\left( \hat{b}_{1}^{\dagger } \right)}^{N_{1,i}^{{\rm B}}}}...{{\left( \hat{b}_{{{n}_{{\rm B}}}}^{\dagger } \right)}^{N_{{{n}_{{\rm B}}},i}^{{\rm B}}}}{{\Phi }_{0}}\end{eqnarray} \tag{ 13 }$$

is the Fock vector for species B only. This matrix is obtained by means of direct diagonalization. In figure 3(c) we show the von Neumann entropy S_A along the transition between BEC–BEC and composite fermionization. S_A can be seen to approach a constant value as ${{g}_{{\rm AB}}}$ is increased, corresponding to the large inter-species correlations present in the composite fermionization. The same plot also shows S_A along the transition between TG–BEC and the phase-separated limit. The two species are less correlated throughout this transition, but still S_A saturates to a constant value in the phase-separated limit. In figure 3(d) we plot S_A for different values of ${{g}_{{\rm AB}}}$ when ${{g}_{{\rm A}}}$ is tuned from zero to a large value. When ${{g}_{{\rm AB}}}=20$, this corresponds to the transition between composite fermionization and a phase-separated gas. We observe a peak which coincides with the crossover between both limits. This peak disappears as ${{g}_{{\rm AB}}}$ is reduced, as observed in the curves for ${{g}_{{\rm AB}}}=4,2$ in figure 3(d). For ${{g}_{{\rm AB}}}=0$, S_A is zero for every value of ${{g}_{{\rm A}}}$.

3.3. Interaction energies

An interesting question is how the interaction energy changes across the transitions between the different limits. The average interaction energy in species A is:

$$\begin{eqnarray}&&\langle {{U}_{{\rm A}}}\rangle =\left\langle \frac{1}{2}\mathop{\sum }\limits_{klmn}\hat{a}_{k}^{\dagger }\hat{a}_{l}^{\dagger }{{{\hat{a}}}_{m}}{{{\hat{a}}}_{n}}V_{klmn}^{{\rm A}} \right\rangle .\end{eqnarray} \tag{ 14 }$$

We display this energy in figure 4(a). For zero ${{g}_{{\rm A}}}$ there are no interactions between A atoms, and $\langle {{U}_{{\rm A}}}\rangle$ is equal to zero. By increasing ${{g}_{{\rm A}}}$, the energy $\langle {{U}_{{\rm A}}}\rangle$ first grows as correlations are being introduced. For larger repulsion, particles avoid each other, which leads to very strong correlations, and the interaction energy drops down to zero. Starting from the BEC–BEC region, this is a long process, however; for a finite value of ${{g}_{{\rm AB}}}$ this happens over a very well defined domain of the parameter ${{g}_{{\rm A}}}$, located at small values of ${{g}_{{\rm A}}}$. Note that for ${{g}_{{\rm A}}}=0$ and in the presence of interaction with species B, the particles in species A are much more localized than for ${{g}_{{\rm AB}}}=0$. Therefore, small increases in the interaction strength ${{g}_{{\rm A}}}$ lead to strong increases in the interaction energy $\langle {{U}_{{\rm A}}}\rangle$. This is also consistent with the maximum found in the correlation strength within component A.

Zoom In Zoom Out Reset image size

The interaction energy goes to zero in the TG–BEC limit, which is the behaviour expected for a single component gas [30–33], as the increased energy is now stored in the single particle harmonic oscillator energies. During the whole process the total energy is increased from:

$$\begin{eqnarray}&&{{E}_{{\rm BECBEC}}}=\frac{1}{2}({{N}_{{\rm A}}}+{{N}_{{\rm B}}})\end{eqnarray} \tag{ 15 }$$

to

$$\begin{eqnarray}&&{{E}_{{\rm TGBEC}}}=\frac{1}{2}({{N}_{{\rm B}}}+N_{{\rm A}}^{2}).\end{eqnarray} \tag{ 16 }$$

The energy for ${{N}_{{\rm A}}}={{N}_{{\rm B}}}=2$ is shown in figure 4(c). The energy obtained by the direct diagonalization and DMC methods coincides. For no interactions between different species, ${{g}_{{\rm AB}}}=0$, and the energy can be expressed as $E=\hbar \omega +{{E}_{2}}({{g}_{{\rm A}}})$, where ${{E}_{2}}({{g}_{{\rm A}}})$ is the energy of two trapped particles interacting with the coupling constant ${{g}_{{\rm A}}}$ [53]. In order to prove that the described limits exist in larger systems, we calculate the energy for ${{N}_{{\rm A}}}={{N}_{{\rm B}}}=10$ particles with the DMC method.

The energy per particle in the BEC–BEC limit (15) does not depend on the number of particles, ${{E}_{{\rm BECBEC}}}/N=1/2$. We show this in figures 4(c) and (e), with green circles for ${{N}_{A}}={{N}_{B}}=2$ and 10.

In figures 4(d) and (f), we depict the energy per particle as a function of ${{g}_{{\rm AB}}}$, starting from the BEC–BEC (solid line) and the TG–BEC (dashed line) limits. Here the green circles (red squares) indicate the energy per atom in the BEC–BEC (TG–BEC) limit. The energy in the TG–BEC limit given by equation (16) is ${{E}_{{\rm TGBEC}}}/N=3/4$ for ${{N}_{{\rm A}}}={{N}_{{\rm B}}}=2$ and ${{E}_{{\rm TGBEC}}}/N=11/4$ for ${{N}_{{\rm A}}}={{N}_{{\rm B}}}=10$. In the transition from the BEC–BEC limit to the composite fermionization one, the energy saturates to a certain value for which we do not have an analytical prediction. Also, a monotonic behavior is observed in the transition from the TG–BEC to the phase separation limit (figures 4(d) and (f)).

The average interaction energy between both species, given by:

$$\begin{eqnarray}&&\langle {{U}_{{\rm AB}}}\rangle =\left\langle \mathop{\sum }\limits_{klmn}\hat{a}_{k}^{\dagger }\hat{b}_{l}^{\dagger }{{{\hat{b}}}_{m}}{{{\hat{a}}}_{n}}V_{klmn}^{{\rm AB}} \right\rangle ,\end{eqnarray} \tag{ 17 }$$

is important to quantify the transition to the composite fermionization or the phase-separated regime. The interaction energy rapidly increases from zero to a maximum at ${{g}_{{\rm AB}}}\approx 2$ (see figure 4(b)) and decreases again toward zero for ${{g}_{{\rm AB}}}\to \infty$. For ${{g}_{{\rm A}}}=0$, this corresponds to building up strong correlations between the particles of different species in the composite fermionization limit, whereas in the limit of large ${{g}_{{\rm A}}}$, this reflects the transition to a macroscopic phase separation of the two components.

3.4. Correlation matrices

Since in the presence of strong interactions the system has non-trivial many-body correlations, it is interesting to look not only at single-particle densities, but also at pair-wise correlation functions. The single-particle densities are quantified by the OBDM (equation (10)). For particles of the same species, the two-particle correlations are quantified by the two-body distribution function (TBDF):

$$\begin{eqnarray}&&\rho _{2}^{{\rm A}}({{x}_{1}},{{x}_{2}})={{N}_{{\rm A}}}({{N}_{{\rm A}}}-1)\int {\rm d}{{x}_{3}}\cdots {\rm d}{{x}_{{{N}_{{\rm A}}}}}{\rm d}{{y}_{1}}\cdots {\rm d}{{y}_{{{N}_{{\rm B}}}}}\;|\Psi {{|}^{2}},\end{eqnarray} \tag{ 18 }$$

with an analogous expression for B. If the two atoms stem from different species, their pair-wise correlations are captured in the cross two-body distribution function (CTBDF) given by:

$$\begin{eqnarray}&&\rho _{2}^{{\rm AB}}({{x}_{1}},{{y}_{1}})={{N}_{{\rm A}}}\;{{N}_{{\rm B}}}\int {\rm d}{{x}_{2}}\cdots {\rm d}{{x}_{{{N}_{{\rm A}}}}}{\rm d}{{y}_{2}}\cdots {\rm d}{{y}_{{{N}_{{\rm B}}}}}\;|\Psi {{|}^{2}}.\end{eqnarray} \tag{ 19 }$$

Both functions are proportional to the joint probability for finding two atoms at two given positions.

It was shown in [43, 46] that the correlation functions are very useful for a description of the composite fermionization and the phase separation limits. In the following, we will carefully look at the transition between these two limits. The phase separation occurs for ${{g}_{{\rm A}}}$ and ${{g}_{{\rm AB}}}$ large and implies a density distribution with atoms of species B localized at the center of the trap, while the atoms of species A gather at the edges of the density of B. As discussed above, the de-mixing point can also be identified in the coherence, the interaction energies, and the entanglement.

In figure 5 we show the OBDMs, TBDFs and CTBDFs just before (${{g}_{{\rm A}}}=5$) and just after (${{g}_{{\rm A}}}=7$), the crossover. The upper and lower rows show numerical results, while the middle row represents the analytical results obtained from ansatz (9) with ${{a}_{{\rm A}}}={{a}_{{\rm AB}}}=0$. One can see that just before the crossover, the densities of both species, i.e., the diagonals of the OBDMs, significantly overlap (panels (a) and (b)), whereas the overlapping is greatly reduced after the crossover (panels (k) and (l)). The TBDFs and CTBDFs before and after the crossover (panels (c)–(e) and (m)–(o), respectively) demonstrate that the atoms of species A are anticorrelated with themselves and with the atoms of species B, as both functions vanish along the diagonal. Note that at the same time, atoms of species B are not strongly correlated. This is also captured by ansatz (9), where strong correlations are induced by zeros whenever A-A or A-B atoms overlap (see panels (f)–(j)). All densities and pair correlations computed with this ansatz qualitatively resemble the exact correlation functions just before demixing. However, the ansatz fails to describe the ground state of the system once it has phase separated.

Zoom In Zoom Out Reset image size

Let us note that the TBDF for the A species shown in figure 5(c) corresponding to the crossover for ${{N}_{{\rm A}}}={{N}_{{\rm B}}}=2$ looks similar to those obtained for ${{N}_{{\rm A}}}=4$ and a very heavy atom in component B (discussed in [56, 57]) or a large number of atoms in B (discussed in [42]). Those cases belong to the phase-separated limit, in which B formed a material barrier. Therefore, the two atoms of A stay at each side of B. Very differently in this case, there are only two atoms of A, and they can be localized in either side of B.

For ${{N}_{{\rm A}}}={{N}_{{\rm B}}}\gt 2$, the results discussed above remain qualitatively valid. We show in figure 6(a) the densities for the composite fermionization limit when ${{N}_{{\rm A}}}={{N}_{{\rm B}}}=2,4,6,8,10$ calculated with DMC. In this situation, the OBDMs are equal for both species. The two peaks present in the density tend to spatially separate as N is increased as a consequence of the large repulsion between both species, which increases with the number of atoms. In figures 6(b) and (c), we show the densities for B and A, respectively, in the phase-separated limit. As N is increased, the atoms of B have a greater tendency to localize in the center of the trap. The numerically calculated density for A shows that this component is localized at each side of B, forming two TG gases with $N/2$ atoms in each side. Note that this is reminiscent of effects in two-component fermionic systems [71], where it was recently shown that a single-particle minority component will be trapped in the center of the majority component due to correlations in the wave function.

Zoom In Zoom Out Reset image size

The difference in the energy between BEC–BEC and TG–BEC regimes is further increased in balanced systems of a larger size, ${{N}_{{\rm A}}}={{N}_{{\rm B}}}\gg 1$. Indeed, according to equation (15), the energy in the BEC–BEC scales linearly with the number of particles N, which is a typical behavior of weakly interacting bosons. Instead, in the TG–BEC limit, according to equation (16) the dependence on N is quadratic. This resembles the behavior of the energy of fermionic particles and is a manifestation of Girardeau mapping. Comparing the results for ${{N}_{{\rm A}}}={{N}_{{\rm B}}}=2$ with ${{N}_{{\rm A}}}={{N}_{{\rm B}}}=10$, we already observe how the difference in the energy between limits increases (see figure 4).

In the imbalanced case ${{N}_{{\rm B}}}\gt {{N}_{{\rm A}}}$, the wavefunction (9) can be equally used as an ansatz for the exact ground state of the systems. The four limits discussed above equally exist. Nevertheless, the weakly interacting species now has a greater tendency to localize in the center of the trap and condense, which modifies the boundaries between the different regimes associated to these limits. In figures 7(a)–(b) and (e)–(f) we report the largest eigenvalue of the OBDM for species A and B to quantify the coherence, covering the whole range of coupling constants, when ${{N}_{{\rm A}}}=2$ and ${{N}_{{\rm B}}}=3,4$, respectively. As ${{N}_{{\rm B}}}$ is increased we observe that the region in which B is not condensed is reduced (the light blue area in figures 7(b) and (f)). Moreover, the minimum value of $\lambda _{0}^{{\rm B}}$, which occurs in this non-condensed area, grows with ${{N}_{{\rm B}}}$ for fixed ${{N}_{{\rm A}}}$. Notice also that the area in which $\lambda _{0}^{{\rm A}}$ approaches the largest possible value ${{\lambda }_{0}}=1$, i.e., close to the ${{g}_{{\rm A}}}$ axis, is reduced as ${{N}_{{\rm B}}}$ is increased.

Zoom In Zoom Out Reset image size

In figures 8(a) and (b), we show the density profiles for A and B along the transition between the BEC–BEC limit and composite fermionization for ${{N}_{{\rm A}}}=2$ and ${{N}_{{\rm B}}}=4$. The atoms of species B are now more concentrated in the center than when both populations were equal, even though species B is still not fully condensed. The two peaks in species A appear at a smaller value of ${{g}_{{\rm AB}}}$ and are more spatially separated than in the case ${{N}_{{\rm B}}}={{N}_{{\rm A}}}$. We note that in the composite fermionization limit, the density of species A in the center for the balanced case is finite, while in the imbalanced case it vanishes (compare figures 2(a) and 8(a)). The density profiles along the transition between the TG–BEC and the phase-separated gas are presented in figures 8(c) and (d). Comparing with the balanced case plotted in figures 2(c) and (d), we notice that, in the phase separation limit, the two peaks in the density profile of A are now more separated, and the squeezing in the density of B is smaller. The average interaction energy $\langle {{U}_{{\rm AB}}}\rangle$ (figures 7(c) and (g)) tends to zero when phase separation occurs. Figures 8(e) and (f) report the density along the transition between composite fermionization and phase separation. We observe that the position of the two peaks in the density profile of A in the phase-separated, and the composite fermionization limit is closer than in the balanced case (compare with figures 2(e) and (f)). Also, the crossover occurs now at a smaller value of ${{g}_{{\rm A}}}$. The average interaction energy $\langle {{U}_{{\rm A}}}\rangle$ (figures 7(d) and (h)) decreases abruptly to zero after the crossover. We conclude that for larger imbalances, ${{N}_{{\rm B}}}\gg {{N}_{{\rm A}}}$, the composite fermionization region is highly suppressed, and, therefore, the surviving limits are those associated to BEC–BEC, TG–BEC, and the phase-separated mixtures.

Zoom In Zoom Out Reset image size

If the macroscopic limit is reached in such a way that the number of atoms in one of the species is fixed, the minority species plays the role of an impurity that perturbs the majority species. The relative contribution of the minority species to the energy becomes smaller, and a polaronic description might be applicable.

Current experimental advances in ultracold atomic physics allow one to scrutinize the onset and evolution of correlations in few-atom bosonic fluids. Small samples can be trapped, and their interactions can be largely tuned, thus providing a fantastic ground to understand how quantum many-body correlations build up in small samples. Binary mixtures are specially appealing, as they provide the first step toward understanding the effect of environments on quantum systems in a controlled way. To advance in that direction, we study the effect of embedding a quantum fluid (component A) within a second quantum fluid (B) with tunable intra- and inter-species interactions at zero temperature. We fix the coupling constant of B-B interactions to that of ideal bosons, ${{g}_{{\rm B}}}=0$, and vary A-A and A-B interactions in a wide range, $0\leqslant {{g}_{{\rm A}}}\lt \infty$, $0\leqslant {{g}_{{\rm AB}}}\lt \infty .$ This permits us to explore the phase diagram for a variety of regimes. The energy, one- and two-body correlation functions, density profiles, and von Neumann entropy are calculated exactly, using the diagonalization method. For larger system sizes, the results are complemented with the energy and density profiles obtained by the diffusion Monte Carlo method.

We have described the transition between the following four limits: (a) the BEC–BEC limit, where both components interact weakly and thus remain condensed, (b) the BEC–TG limit, where the two components interact weakly among each other and A has strong intra-species interaction, (c) the composite fermionization limit, where the interaction between both species is large, inducing strong correlations within both species, and (d) a phase separation limit, where both the intra-species interaction in A and the inter-species interactions are large. We show that the transition between the different limits involves sophisticated changes in the one- and two-body correlations. The energetic properties change in a smooth way, with the energy and its derivatives remaining continuous, which implies a transition of a crossover type rather than a true phase transition. At the same time, the entanglement between the two components has a much sharper dependence on the interactions. This is demonstrated by reporting the von Neumann entropy, which manifests a sharp peak along the transition between composite fermionization and phase separation. The evolution of the density profiles of A and B components is studied in detail both for the balanced and the imbalanced case. The effect of a large number of particles on the energy and the density profiles is discussed. We analyze the coherence properties by expanding the one-body density matrix in natural orbitals and obtaining the occupation numbers. We demonstrate that full condensation (largest occupation number equal to one) for A species is reached only in the BEC–BEC regime, while the weakly interacting B species also remains fully condensed in the TG–BEC regime, and the condensation is almost complete in the phase separation regime. We argue that the described picture of the transition between four mentioned regimes remains valid also in macroscopically large balanced mixtures, ${{N}_{{\rm A}}}={{N}_{{\rm B}}}\to \infty$. Conversely, when the macroscopic limit is reached by increasing the number of atoms of the weakly-interacting species ${{N}_{{\rm B}}}\to \infty$, the composite fermionization limit is suppressed. Therefore the phase diagram in this highly imbalanced case resembles that expected within a mean-field approach. The studied effects are relevant to ongoing and future experiments with small two-component systems.

This project was supported by Science Foundation Ireland under Project No. 10/IN.1/I2979. We acknowledge also partial financial support from the DGI (Spain) Grants No. FIS2011–25275, FIS2011–24154, and the Generalitat de Catalunya Grant No. 2009SGR-1003. GEA and BJD are supported by the Ramón y Cajal program, MEC (Spain).