Time-of-flight expansion of binary Bose-Einstein condensates at finite temperature

Ultracold quantum gases provide a unique setting for studying and understanding the properties of interacting quantum systems. Here, we investigate a multi-component system of $^{87}$Rb--$^{39}$K Bose-Einstein condensates (BECs) with tunable interactions both theoretically and experimentally. Such multi-component systems can be characterized by their miscibility, where miscible components lead to a mixed ground state and immiscible components form a phase-separated state. Here we perform the first full simulation of the dynamical expansion of this system including both BECs and thermal clouds, which allows for a detailed comparison with experimental results. In particular we show that striking features emerge in time-of-flight for BECs with strong interspecies repulsion, even for systems which were separated in situ by a large gravitational sag. An analysis of the center of mass positions of the BECs after expansion yields qualitative agreement with the homogeneous criterion for phase-separation, but reveals no clear transition point between the mixed and the separated phases. Instead one can identify a transition region, for which the presence of a gravitational sag is found to be advantageous. Moreover we analyze the situation where only one component is condensed and show that the density distribution of the thermal component also show some distinct features. Our work sheds new light on the analysis of multi-component systems after time-of-flight and will guide future experiments on the detection of miscibility in these systems.

A key property of multi-component systems, is their miscibility . Miscible components overlap in space leading to a mixed ground state, while immiscible components form a phase-separated state. Experimentally, it is often difficult to detect the miscibleimmiscible transition by direct in situ imaging of two trapped components. Therefore most experiments have investigated the transition between mixed and separated phases by observing the time-of-flight (TOF) expansion of a mixture of Bose-Einstein condensates (BECs) [8,10,14,17,18]. These experiments have typically identified the transition by observing a dramatic increase in the separation between the centre of mass (COM) positions as the interaction strength was changed, with associated strong features at the interface between the expanded mixtures. The observed transition points were found to be in good quantitative agreement with those predicted by the criterion for the transition from miscible to immiscible in a homogeneous system. This criterion is given by the miscibility parameter ∆ = (g 11 g 22 /g 2 12 ) − 1 which depends on the intra-(g 11 , g 22 ) and interspecies (g 12 ) interaction strengths and determines whether two components mix (∆ > 0) or phase-separate (∆ < 0). However, we have previously shown that this criterion is not applicable in the trapped case and proposed an alternative criterion based on the in situ density distributions, which can be probed through the observation of dipole oscillations [44]. These simulations further indicated that the interactions significantly affect the expansion dynamics of these binary systems, and hence details of the relation between the in-situ density distributions and the COM positions after TOF expansion remained unclear.
In this work, we perform the first full simulation of the dynamical expansion of a multicomponent system including both BECs and thermal clouds and compare these results with experiments on a mixture of 87 Rb-39 K atoms. We show that striking features emerge in TOF for BECs with strong interspecies repulsion, even for two systems which were geometrically separated in situ by their large trap offsets. Moreover we show that a measurement of the centre of mass position of binary BECs after TOF expansion yields a result in qualitative agreement with the homogeneous criterion for phase-separation. In both cases we obtain good agreement between our simulation and the experimental results. This reconciles our previous theoretical analysis [44], with experimental observations [8,10,14,17,18].
In addition our analysis reveals a multitude of other interesting features. If only one component is condensed, the density distribution is shown to retain distinct features. Although the BEC in the partly-condensed component is found to be largely symmetric after expansion, the non-condensed species exhibits asymmetric features, which are enhanced in the immiscible case. Moreover, in the case of both components being partly-condensed, we show that a trap offset (e.g. due to a gravitational sag) is useful in yielding a pronounced shift of the COM position after expansion, thus facilitating an easier characterization of the properties of mixtures of different species. This was in fact implicitly used already in early binary mixture experiments (see e.g. Fig. 2 in Ref. [8]). Finally, we show that there is no abrupt transition which can be observed in TOF.
In the remaining sections, we first outline the experimental procedure used to investigate binary mixtures (Sec. 2.1) and give an overview of the theoretical tools used to model this scenario (Sec. 2.2). We then perform a detailed analysis of the spatial distributions and COM positions of binary mixtures after TOF expansion (Sec. 3). This includes the effect of temperature and of the gravitational sag on TOF expansion. We summarize our key findings in section 4.

Experimental and Theoretical Methodology
In the following a detailed description of the experimental production and detection of binary BECs is given. Moreover we provide an overview of the theoretical methods which enable us to perform a full numerical simulation of the TOF expansion of binary BECs in the presence of a thermal clouds. Further technical details are given in Appendix A.

Experimental realization of binary BECs
Our experiments with binary mixtures are performed with 87 Rb and 39 K atoms in an optical dipole trap. The optical dipole trap allows for free tuning of the external magnetic field and thus Feshbach resonances can be employed to tune the scattering length between the two species. Moreover the optical trap allows for particularly rapid switching of the external potential and thus a precise comparison of the TOF dynamics with theoretical results is possible. However, due to the different masses, the gravitational sags of the two species differ significantly. The effect from this relative trap offset is analysed further in Sec. 3.5. Figure 1 shows the relevant scattering properties. We perform the experiments with both species in the |F = 1, m F = −1 state in a range of magnetic fields where the scattering length of 39 K is positive, allowing for the production of binary BECs [17]. Moreover an interspecies Feshbach resonance is available in this region and thus the miscibility parameter can be tuned across the transition point. At a given magnetic field B, the s-wave scattering length between 39 K atoms is described by and the scattering length between 39 K and 87 Rb atoms is where a 0 is the Bohr radius [17,45]. The scattering length between 87 Rb atoms is constant at a Rb−Rb = 99a 0 . The interaction strengths are related to the scattering lengths as g i j = 2πh 2 a i j (m i + m j )/(m i m j ), where the indices refer to the two species with respective masses m i . This allows us to study mixtures with miscibility parameters ∆ = g 11 g 22 /g 2 12 − 1, from highly-miscible (∆ ∼ 1) to highly-immiscible (∆ ∼ −1), as shown in Fig. 1.
The details of the experimental apparatus and procedure are described in [17] and briefly summarized in the following. Initially, a dual-species magneto-optical trap simultaneously captures both species of atoms from a background vapour. This is followed by an optical molasses which cools them to sub-Doppler temperatures. Afterwards both species are optically pumped into the fully stretched |2, 2 state. The atoms are loaded into a magnetic quadrupole trap which transports them to a separate chamber where the remainder of the experimental procedure is performed. Here, they are loaded into a different quadrupole trap and evaporative cooling is performed. Microwave radiation is used to selectively cool 87 Rb atoms, which cool 39 K sympathetically. Next, an additional coil transforms the trap into a harmonic Ioffe-Pritchard configuration where further microwave cooling is performed. Before the atoms reach quantum degeneracy, they are loaded into an optical dipole trap consisting of two crossed laser beams at a wavelength of 1064 nm. A rapid adiabatic passage using radio frequency radiation transfers both species from the |2, 2 to the |2, −2 state. From here a radio frequency π-pulse and a rapid adiabatic passage transfer 87 Rb and 39 K atoms into the target state |1, −1 respectively.
The final evaporation is performed by lowering the power of the dipole trap laser, which preferably allows 87 Rb to leave the optical trap due to its larger gravitational sag. The rethermalisation of the two species during evaporation, is enhanced by setting the magnetic field to 118.7 G, which is in the vicinity of the interspecies Feshbach resonance. After evaporation, the field is adjusted to a desired target value, and the dipole trap power is increased to ensure further rethermalisation, resulting in trap frequencies of (ω ρ , ω z ) K /2π = (132, 186) Hz and (ω ρ , ω z ) Rb /2π = (93, 124) Hz for 39 K and 87 Rb, respectively. At these trap parameters the gravitational sag separates the centres of the potentials for the two species leading to an estimated trap offset of 7.1 ± 0.7 µm.
Finally, the TOF detection is realized by turning off the optical trap, and allowing the atoms to expand freely. The magnetic field is kept at the desired value for 8 ms allowing the atoms to interact during the initial expansion where interactions are significant. Subsequently the field is turned off and absorption images of 39 K and 87 Rb atoms are taken after 23 ms and 24.7 ms respectively. The time between these images is set by technical limitations of the camera system. Examples of absorption images of the two components are shown in Fig. 1 (c-f).
By tailoring the evaporation sequence in the dipole trap appropriately, it is possible to condense both species, only 87 Rb, or none of the two. When condensing both species, we typically obtain about 3 × 10 4 39 K and 1.2 × 10 5 87 Rb atoms at a temperature of 150-200 nK, resulting in BECs of respectively 8 × 10 3 39 K and 5 × 10 4 87 Rb atoms. The sum of the Thomas-Fermi radii of the condensates in the vertical direction is approximately 6 µm, which is on the order of the estimated trap offset.
To evaluate how the interspecies interaction affects the density distributions, we additionally record images of each species without the other one present. This allows us to evaluate the undisturbed COM position and spatial distribution of each species. These undisturbed COM positions are used as the origin of the coordinate systems for each species in all experimental data shown.

Theoretical model
The finite temperature dynamics of a partially-condensed system is well described in the context of a "two-gas" model ("Zaremba-Nikuni-Griffin", or "ZNG" model [46,47]), consisting of a BEC and a thermal cloud. Here we follow our earlier work [48][49][50] which appropriately generalized this to a binary mixture in such a way that both BECs and both thermal clouds are coupled within and between the mixture components through both meanfield and collisional interactions. In order to probe expansion dynamics, we limit our discussion here to the collisionless limit of the above theory: in this limit, each BEC is described by a generalized Gross-Pitaevskii equation (GPE), where j denote the component. The thermal clouds are described by a collisionless quantum Boltzmann equation where f j (p, r,t) = dr e ip · r /h δ † (r + r /2,t) δ (r − r /2,t) denotes the Wigner distribution of the thermal atoms (withδ j =Ψ j − φ j the fluctuation operator after the BEC mean field, φ j has been subtracted). The local thermal cloud density is then given byñ j (r) = dp/(2πh) 3 f j (p, r,t) while the density of the whole cloud n tot, j is obtained as the sum of the two, i.e. n tot, j = n c, j +ñ, where n c, j = |φ j | 2 . The total number of atoms (which is fixed in our model) is obtained as N j = dr n tot, j (r).
In contrast to the usual zero-temperature case, the BEC atoms now experience an effective potential U j c corrected by the presence of the thermal cloud, while the thermal atoms are modelled as classical particles moving in an effective potential U j n . These potentials, which are found to play a key role during the early stages of the expansion process, include both the external potential V j and the mean-field contributions, which are related to the BEC density n c, j and the thermal densityñ j , by The trapping potentials used here are defined by with trap centres z s, j , and radial and axial angular frequencies, ω ρ, j and ω z, j , respectively. The above expressions are consistent with standard multi-component Hartree-Fock theory [51], which includes an extra factor of 2 in the condensate mean-field potential (Eq. (5a)) associated with the direct and exchange interactions involving thermal and condensate atoms of the same species. A corresponding factor of 2 is included in the total density contributions appearing in the potential for the thermal atoms (Eq. (5b) within the same component. Such factors are however absent in the coupling between the components (final contributions in each of Eq. (5a-b)).
The experimental scenario is simulated as follows. The initial equilibrium states are obtained in the usual way [48][49][50]52]. At time t = 0, the atoms are released from the trap (i.e. V j (r) = 0) and expand freely. The expansion dynamics is simulated in a fixed twodimensional grid utilizing the cylindrical symmetry of the problem which allows us to study atomic clouds that have expanded up to approximately ten times their in-trap sizes without changing the grid spacing. Further expansion is computationally prohibitive due to the need to accurately capture the length scales associated with the non-negligible expansion speeds, as discussed in Appendix A, which also gives further details on the implemented numerical approach. Similar to the experimental procedure where the magnetic field is changed during expansion for imaging purposes, the scattering length in the simulation is abruptly changed after 6 ms, which is however found to have no noticeable effect on the COM position of the already diluted expanded clouds. Further details of our numerical implementation and the effect of scattering length change are given in Appendix A and Appendix B respectively.
The simulation parameters are based on our previous experimental work [17]. We thus consider a mixture of 1.2 × 10 5 87 Rb atoms and 4.2 × 10 4 39 K atoms in a harmonic trap with trapping frequencies of (ω ρ , ω z ) K /2π = (136, 189) Hz and (ω ρ , ω z ) Rb /2π = (96.9, 129) Hz. In this situation the trap centres are separated by an offset of 6.9 ± 0.5 µm due to the gravitational sag. Both species are assumed to be at a relatively high temperature T = 200 nK (unless otherwise stated) for which there is approximately 15% condensed fraction in 87 Rb atoms and 26-30% condensed fraction in 39 K atoms, with exact numbers slightly dependent on the precise values of interaction strengths. The scattering lengths are set by the magnetic field through Eqs. (1)-(2), and both finite-size and mean-field corrections are taken into account when calculating the critical temperatures for the two species [2].
In the following sections we compare our simulations to experimental data obtained according to section 2.1. The relevant parameters are similar, but the data was recorded using considerably longer expansion time, which allows spatial features to be resolved more precisely. Consequently, the experimental density distributions are approximately twice as large as the simulated distributions. However, this has no influence on main arguments and conclusions. Section 3.2 is based on a comparison with experiments performed at relatively short time of flight which were previously presented in [17].
For displaying our theoretical results, we use the centre between the two traps as zero point, unless otherwise stated.

Time-of-flight investigation of miscibility
The experimental realization of binary BECs of 87 Rb and 39 K in an optical dipole trap and the availability to tune the interaction strength allow for a detailed experimental investigation of miscibility in this system. Our simulations, which include the full dynamical evolution of the two density distributions and their mutual mean-field effects, can thus be compared to the experiment directly.

Density distributions of binary Bose-Einstein condensates after expansion
The simulated density distributions before and after TOF expansion for both a miscible and an immiscible mixture, are shown in Fig. 2, alongside the corresponding experimental images. A visual comparison of the distributions after TOF in the miscible and immiscible limits, immediately shows the good agreement between theory and experiment.   The in-situ distributions do not show any distinguishing feature between the miscible and immiscible mixtures for the chosen trap. The atomic clouds are generally elliptical with widths on the order of the single-component Thomas-Fermi radii, which have a larger value along the x-axis due to of the smaller radial trap frequencies. This lack of features arises from the fact that the BECs are separated by a distance, which is similar to the sum of the Thomas-Fermi radii.

Simulation Experiment
During the expansion of individual BECs the aspect ratio typically inverts [53,54], leading to a larger width along the axial direction. However, the repulsion between the two species strongly alters the dynamics and a flat interface between the two species is observed in our simulations and experiments as shown in Fig. 2 (b,c,e,f). Similar structures have previously been observed in other studies [8,10,14,17,18]. In particular in the case of an immiscible mixture, the 39 K BEC is strongly compressed and deviates dramatically from the predicted size due to self-similar expansion [53] as shown in Fig. 2 (b,c). Similarly the 87 Rb BEC is compressed with a flat interface towards the 39 K BEC. However, even in the miscible regime, the repulsive interactions are sufficient to affect the spatial distribution of the 39 K BEC after expansion as seen in Fig. 2 (e,f). Figure 3 shows the corresponding doubly-integrated densities. In the immiscible case (top row) the spatial features arising due to the coupled expansion are clearly visible, both in simulations and experiments. On the one hand, the density of 87 Rb exhibits a steep shoulder at the interface with the 39 K BEC. Moreover, distinct secondary peaks can be seen in the densities of both components in the directions facing away from the mixture interface. These are caused by the strong repulsion during the initial expansion, which accelerates the two components away from each other. For 39 K this effect can be rather pronounced, with the furthest fraction appearing to be almost spatially separated from the main 39 K BEC (see also Fig. 2 (c)). These features are reminiscent of dispersive shock waves, predicted to arise from the strong mutual expulsion of the two components [55,56].
To assess the effect of the thermal clouds we compare our full simulation with a simulation based on the GPE at T = 0. Figure 3 (dashed lines) shows that the presence of the thermal clouds partially suppresses the sharp features in the density distributions after coupled expansion. Nonetheless the same striking features are clearly visible in both simulated and the experimental density distributions.

Effects on the centre of mass positions
The COM position along the axial direction (z-axis) is a particularly suitable quantity to investigate the miscible-immiscible transition.
Following the numerical method outlined above we extract the relative COM positions of binary BECs and the associated thermal clouds after a full simulation of the TOF expansion. Since our computational grid coincides with a moving frame free-falling under gravity our results are given in terms of a relative COM position z . Figure 4 shows the time evolution of the relative COM positions for all components in the immiscible and the miscible case.
In general, the two species accelerate away from each other during the first 6 ms due to the interspecies repulsion, but travel at an approximately constant speed afterwards as the densities become negligible. We have checked that our simulated COM positions (including other sets of numerical data not shown here) obey the conservation of total momentum. It is interesting to note that the BEC atoms acquire a greater COM displacement from the interspecies repulsion than the thermal atoms. Since the BEC densities are larger and their distributions are spatially more compact, they experience a larger acceleration due to the interspecies repulsion than the thermal clouds. Figure 5 shows the extracted relative COM positions of the BECs compared to experimental data from [17]. In the simulation we extrapolate the COM positions to the experimental time of flight by performing a linear fit to the COM positions between 6 ms and 14 ms. We observe good agreement between simulation and experiment, confirming that the qualitative behaviour during TOF is captured well. However, a relatively small offset of 5.5 µm is required in the simulations in order to reproduce the experimental results. In Appendix C we additionally show the qualitative behaviour of the COM positions is independent of the offset.
We have also investigated whether it is physically meaningful to extract a transition point from this relative position. Specifically, we consider the derivative of the relative COM positions with respect to the miscibility, as discussed in more detail in Appendix C. We find no sharp increase of this derivative, which increases monotonically from the miscible to the immiscible regime.

Relevance of homogeneous miscibility criterion for trapped mixtures
To relate our findings to the homogeneous BEC phase-separation criterion, we restrict our discussion here to the case of two pure BECs (T = 0). Firstly, we recall that the criterion ∆ = 0 has been derived for a homogeneous T = 0 in-situ mixture at equilibrium [57], whereas the experimental data corresponds to non-equilibrium expanding density distributions of two BECs which were initially confined in a harmonic trap with different trapping frequencies for the two species. Upon release from the harmonic trap, the BECs expand, with each BEC affected by the mean field of the other one [8]. The shift in relative COM position arises primarily from the interaction within the interface region between the two species. We therefore argue that the transition region is predicted by the mechanical equilibrium of the BEC atoms in this interface region, i.e.

∂U Rb
which leads to the usual criterion ∆ = (g Rb−Rb g K−K /g 2 Rb−K ) − 1 = 0. This confirms that the agreement of the experimental observations [8,10,17,18] with the homogeneous criterion ∆ = 0 has its origin in the expansion dynamics, rather than in the in-situ density distribution.

Miscibility of thermal clouds and BECs
The two species in a mixture will generally have different critical temperatures, implying that the presence of a BEC in one component does not guarantee a BEC in the other. It is therefore interesting to investigate to which extent the presence of a BEC in one component affects the expansion dynamics of a BEC or thermal cloud in the other component coupled to it. This can also shed more light on the miscibility of a binary partially condensed mixture.
In our experiments, the typical range of temperatures and atom numbers is such that it enables us to probe all three cases below in the miscible and immiscible limits: Both (partly) condensed: T Rb < T c,Rb and T K < T c,K Rb (partly) condensed; K thermal: T Rb < T c,Rb and T K > T c,K Both thermal: T Rb > T c,Rb and T K > T c,K The density distributions of the two species after expansion for different miscibilities and temperatures are shown in Fig. 6. A comparison of our numerical simulations (left columns) with the corresponding experimental measurements (right columns) shows excellent agreement for both immiscible (top rows) and miscible (bottom rows) mixtures.
Comparing the three numerical density distributions in the immiscible regime Fig. 6 (top left panels), there is a considerable difference in the axial width and shape of the two expanded systems for different temperature regimes. In particular, if the temperature of 39 K is above the critical temperature, the flat repulsive interface disappears in the 87 Rb distribution but remains visible in the central part of the thermal 39 K distribution. If the temperature of 87 Rb is also above the critical temperature, the density distributions are circular in shape, reflecting the isotropic nature of the velocity distribution of the thermal atoms. These features are also present in our experimental TOF images (top right panels).
We also investigate the interplay of miscibility and temperature, by comparing expanded distributions in different temperature combinations for systems expected to be immiscible ( Fig. 6 top panels) and miscible (bottom panels). When both species are partially condensed, strong features between miscible and immiscible mixtures are observed as discussed in Fig. 2. Moreover, we find a slightly-dented 39 K thermal density for the immiscible mixture when only the 87 Rb atoms are partially condensed and no difference when both systems are thermal. This shows that the TOF expansion has little influence on the low density thermal atoms.
In both immiscible and miscible cases, experimental observations (Fig 6 right columns), reveal broadly the same features as in our simulations discussed above. This includes a dented 39 K thermal distribution in the presence of a strongly repulsive 87 Rb BEC.
We explain these observations by the difference in the mean-field forces, which are proportional to the density gradients (∇U j c and ∇U j n ). For example, assuming that the BEC density can be approximated by a Thomas-Fermi distribution ∝ (z − z s,k ) 2 and the thermal density is approximated by a Gaussian distribution ∝ exp[−m k ω 2 z,k (z − z s,k ) 2 /(2k B T k )] , we expect the contribution to the mean-field forces to scale like z − z s,k from a BEC, and (z − z s,k ) exp[−m k ω 2 z,k (z − z s,k ) 2 /(2k B T k )] from a thermal cloud, where the indices j and k refers to the two species. Since the two species will meet and interact at a distance z far from z s,k , we expect the contributions from the BEC to be large while those from the thermal cloud will be minimal due to the extra exponential factor. Hence the presence of the BECs can be inferred from the density distribution of the other species after TOF.
Our analysis confirms that the repelling features seen in Fig. 2 are a result of the mutual BEC mean field repulsion during expansion (since the original BECs were spatially separated), and that the vertical width of the 39 K BEC is significantly compressed by the repulsion. In case only one species is condensed, the difference between miscible and immiscible interactions is less distinguishable.
Although such behaviour is broadly expected, the excellent verification of our numerical model against experiments enables us to make further numerical predictions about the importance of the trap offset (e.g. due to gravitational sag) in experiments.

Effect of the gravitational sag
The comparison of our experimental data with the full simulation in the preceding sections has shown that it is vital to appropriately include the initial gravitational sag between the two species. The role of trap sag has previously been theoretically investigated in Refs. [30,31], while the effect of temperature has been considered in Refs. [27,42,[48][49][50]. Figure 7 provides a more complete analysis of the effect of the sag on the COM position of the expanded partially condensed 39 K cloud. In both the miscible and the immiscible case, the repulsion of the clouds leads to a shift of the COM position. Not surprisingly this shift is larger for small sag, corresponding to a larger initial overlap of the repulsive clouds, and tends slowly to zero with increasing distance between the two trap centres. This suggests that a small nonzero sag, e.g. one that is considerably smaller than each of the individual Thomas-Fermi radii (see Refs. [8,10]) is ideal for achieving a maximal relative COM position in experiments. Comparing the COM shift between miscible and immiscible regimes in Fig. 7 for fixed sag and expansion time, we see a clear enhancement for immiscible clouds, due to the stronger intraspecies scattering length g 12 which generates a stronger mutually-repulsive potential during expansion. Our simulations also enable us to address the effect of an initial sag on the density distributions of the expanded clouds. Typical examples of such distributions are shown in Figure 7 for two different sag values in both immiscible (top) and miscible (bottom) regimes. In each case both the in situ and the expanded distributions are shown. In the miscible case neither the in situ, nor the expanded distributions are significantly affected by the sag. Contrary to this, a change in the sag significantly affects the distributions in the immiscible case. Denoting the BEC radii (finite temperature Thomas-Fermi radii) by R Rb and R K , and the trap sag by d, we distinguish two cases: (i) For small trap offset d < R Rb , R K the in situ distributions already reveal evidence of immiscibility; upon expansion, there is maximum relative COM position, with the curvature of the Rb distribution in the interface region becoming smoother, but not quite planar, while the repelled expanding 39 K BEC remains largely spherical. (ii) Once the trap offset becomes comparable to the sum of the effective BEC radii d (R Rb + R K ), we see a very different picture: the in situ distributions are largely decoupled due to their large geometrical separation, having elliptical shapes (ellipticity set by the trap aspect ratio). However, upon expansion the interface between the two BECs becomes largely planar, and the widths of both BECs in the vertical direction are highly compressed by the mean-field potentials during expansion. It is precisely those distributions that were analysed in Fig. 2, which confirmed the experimental validity of our findings.
We thus conclude that although a small trap sag may be ideal for maximizing the COM displacement between two immiscible BECs, a cleaner signature of the immiscibility of a given mixture is obtained for a larger sag, with the dynamical mean-field pressure upon expansion generating a planar interface between the two initially-disjoint BECs. Moreover the strong acceleration leads to highly non-trivial density distributions with clearly distinguishable side peaks.

Conclusion
To summarize, we have performed the first ab initio theoretical analysis of two components which fully includes the mean-field dynamics during TOF expansion, showing very good agreement with experimental results. Our analysis explains the striking features which emerge in time-of-flight for BECs with strong interspecies repulsion, showing that the expansion can lead to the formation of distinct spatial features, which closely resemble dispersive shock waves [55]. Moreover our analysis of the centre of mass positions of the BECs after expansion demonstrates that the homogeneous phase-separation criterion also emerges in mixtures during expansion as a result of mechanical equilibrium across the BEC interface region. This is an important observation since it reconciles the previously-reported lack of applicability of this criterion for in situ density distributions of inhomogeneous mixtures with the experimental mixture observations reported to date. However, we find that no clear transition point between the mixed and the separated phases can be detected in the centre of mass positions after expansion, but that this is a more gradual effect, where g 11 g 22 ≈ g 2 12 marks an approximate transition between miscible and immiscible. The presence of a gravitational sag is found to be advantageous for identifying this transition region. Moreover, we find that this sag should be a fraction of the combined Thomas-Fermi radii of the two BECs to obtain a large shift of the COM position. Finally we analysed the situation where only one component is condensed. In this case distinct features prevail, with the BEC largely symmetric after TOF expansion, but a slightly asymmetric non-condensed species, which is more pronounced in the immiscible case.
This work sheds new light on the evaluation of multi-component systems after time-offlight. In particular it shows that the expansion has a major effect on the density distributions and COM position observed after TOF. Thus our analysis will guide future experiments on the detection of miscibility in these systems. It can be discretized on a two-dimensional spatial grid,

Acknowledgments
with the radial grid spacing dρ, the axial grid spacing dz and a total of N ρ × N z grid points. The time-propagation is then solved with alternating direction implicit method. In particular, we use the Crank-Nicolson method [58] for the radial direction and the Fourier spectral splitstep method [59] for the axial direction.
In order to apply the Crank-Nicolson method, we approximate the radial differential operators with a finite-difference approximation. Using the notation φ j (l, m) ≡ φ j (ρ l , z m ), we arrive at 1 We further choose φ j (0, m) = φ j (1, m) to enforce the boundary condition ∂ φ j /∂ ρ = 0 at ρ = 0 such as to preserve the cylindrical symmetry. Equations (A.6) and (A.7) then become exact at ρ = dρ/2 if φ j = c 0 + c 2 ρ 2 for any constant c 0 and c 2 .
On the other hand, the quantum Boltzmann equation (4) is solved with the direct simulation Monte Carlo method [60,61]. A large number of test particles (typically of the order of millions) are sampled with the acceptance-rejection method according to the Bose-Einstein distribution [e (ε j p − µ j )/(k B T j ) − 1] −1 , where ε j p = p 2 2m j + U j n is the local energy of a thermal atom. These test particles are then evolved in time according to Newton's equation of motion using the symplectic leapfrog method.
In contrast to the BEC propagation, simulating the test particle dynamics with cylindrical polar coordinates requires a computational grid with irregular radial grid spacings, where the radial grid pointsρ l are located at a constant multiple of the zeros ζ l of the Bessel function J 0 . The main reason is the need for an efficient numerical method to compute the discrete Hankel transformation (see details in the next paragraph). Arranging the zeros ζ l in increasing order, with ζ 1 being the smallest zero (ζ 1 = 2.4048), we havẽ ρ l = ζ l ζ N ρ +1 ρ max , for l = 1, 2, . . . , N ρ (A. 8) where ρ max gives the radial boundary of the computational grid. Three technical points are worth mentioning here: (1) binning of the test particles into spatial cells, (2) smoothing of the resulting distributions so as to approximateñ j and (3) computing the mean-field force acting on the test particles from the spatial derivative of the effective potential U j n .
(i) Binning of test particles Suppose that we are simulating the dynamics ofÑ j thermal atoms using N tp, j test particles. Each test particle therefore carries a weight w j =Ñ j /N tp, j . For a test particle with a position vector (x, y, z) ≡ (ρ, z), we assign its contribution to D j (ρ, z) ≡ 2πρñ j (ρ, z) using the cloud-in-cell method: forρ l−1 < ρ <ρ l and z m−1 < z < z m , we add the contribution

Appendix C. Characterization of the miscible-immiscible crossover
To investigate whether it is physically meaningful to extract a transition point from the relative COM position of expanded clouds, we interpolate our simulated data points by cubic splines (Fig. C1 top) and extract the first-order derivative of the relative COM positions with respect to the miscibility (Fig. C1 bottom). A well defined transition point should manifest itself as a sharp increase of the derivative. While the relative shift is continuously increasing from the miscible (∆ > 0) to the immiscible (∆ < 0) regime, the derivative does not show a sharp transition point. Fig. C1 also shows that the expansion dynamics depends on the chosen offset. The best agreement between theory and experiment is obtained for an offset of 5.5 µm in our simulation, despite the fact that the estimated experimental offset is 6.9 µm. Several sources contribute to this discrepancy, including inaccurate knowledge of the precise scattering lengths, a systematic error in the evaluation of the centre of mass condensate positions and incomplete knowledge about the precise trap offset.