Correlation Effects in the Quench-Induced Phase Separation Dynamics of a Two-Species Ultracold Quantum Gas

We explore the quench dynamics of a binary Bose-Einstein condensate crossing the miscibility-immiscibility threshold and vice versa, both within and in particular beyond the mean-field approximation. Increasing the interspecies repulsion leads to the filamentation of the density of each species, involving shorter wavenumbers and longer spatial scales in the many-body approach. These filaments appear to be strongly correlated and exhibit domain-wall structures. Following the reverse quench process multiple dark-antidark solitary waves are spontaneously generated and subsequently found to decay in the many-body scenario. We simulate single-shot images to connect our findings to possible experimental realizations. Finally, the growth rate of the variance of a sample of single-shots probes the degree of entanglement inherent in the system.


I. INTRODUCTION
The realm of atomic Bose-Einstein condensates (BECs) has offered over the past two decades a fertile testbed for the examination of phenomena involving the role of nonlinearity in wave dynamics and phase transitions [1][2][3][4][5][6]. Phase separation dynamics in the case of multi-component BECs has held a prominent role among the relevant studies and is a topic that by now has been summarized in various reviews [1,2,4,7]. Nevertheless, the majority of the relevant studies has focused on a mean-field (MF) description, while the role of many-body (MB) effects in such transitions is much less understood.
Since the early days of the experimental realization of BECs, experimental achievements include binary mixtures of e.g. two hyperfine states of 23 Na [8] and of 87 Rb [9]. Progress of the experimental control over the relevant multi-component settings enabled detailed observations of phase separation phenomena and related dynamical manifestations [10][11][12][13][14][15][16][17][18]. In recent years, external coupling fields have been utilized to control and modify the thresholds for mixing-demixing dynamics in pseudo-spinor (two-component) [19,20] and even in spinor systems [21]. Moreover, the quench dynamics across the phase separation transition has been a focal point of studies examining the scaling properties of suitable correlation functions and associated universality properties [22][23][24].
More recently the inclusion of correlations in multicomponent few boson systems enabled a microscopic characterization of their static properties. A variety of novel phases have been realized in these settings such as altered phase separation processes [25][26][27][28][29], composite fermionization [30][31][32], or even the crossover between the two [33,34]. Also the dynamical properties of such MB ultracold mixtures have been studied including, among others, the dependence of the tunnelling dynamics on the mass ratio [35,36] or the intra-and interspecies interactions [37], as well as the emergence of Anderson's orthogonality catastrophe upon quenching the interspecies repulsion [38]. On the other hand, far less emphasis has been placed on the MB character of the quench-induced phase separation phenomenology. It is the latter apparent gap in the literature that the present work aims at addressing for both few and larger bosonic ensembles.
To incorporate the quantum fluctuations due to correlations [39][40][41][42] emerging when quenching the binary BEC system, we bring to bear the Multi-Layer Multi-Configuration Time-Dependent Hartree Method for bosons (ML-MCTDHB) [43,44] designed for simulating the quantum dynamics of bosonic mixtures. We explore different scenarios, emphasizing the case where the inter-species interaction is quenched from the miscible to the immiscible regime (positive quench) or vice versa (negative quench). We find significant variations in the MB scenario in comparison to the MF one. In the positive quench scenario the unstable dynamics leads to the filamentation of the density of each species and the dominant wavenumber associated with the emerging phase separated state appears to generically be higher in the MF case. The one-and the two-body correlation functions indicate the presence of correlations between the filaments of the same or different species signalling the presence of fragmentation and entanglement respectively. In particular, strong one-body correlations appear between non-parity symmetric (with respect to the trap center) filaments formed indicating their tendency of localization. These filaments are found to be strongly anticorrelated at the two-body level indicating a negligible probability of finding two bosons of the same species one residing in an outer and one in an inner filament. More importantly, combining the behavior of one-and twobody correlations supports the formation of domain-walls i.e. interfaces that separate these distinct filaments [45][46][47].
In sharp contrast to the above dynamical manifestation of the phase separation, in the negative quench scenario multiple dark-antidark (DAD), i.e. density humps on top of the BEC background, solitary waves [48,49] are spontaneously generated both within and beyond the MF approximation. At the MB level many decay events, at the early stages of the dynamics, increase the production of DAD solitary waves with the product of each decay being a slow and a fast DAD structure [50]. The latter increase results in multiple collisions and interference events between these matter waves, and most of them are lost during evolution. Furthermore, in both the positive and the negative quench scenarios, single-shot simulations, utilized here for the first time for binary mixtures, offer a link to potential experimental realizations of the above-observed dynamics. In particular, the growth rate of the variance of single-shots resembles the growth rate of the entanglement inherent in the system. Additionally, deviations between the variances of the two species reveal the fragmented nature of the binary system. Last, but not least the case of quenches within the immiscible regime, are explored showcasing the one-dimensional (1D) analogue of the so-called "ball" and "shell" structure appearing in higher-dimensional binary BECs [12].
Our presentation is structured as follows. In section II, we provide the details of the binary setup and the corresponding MB ansatz, briefly addressing the ML-MCTDHB approach. In section III we examine the different quench scenarios focusing on the miscible to immiscible quench as well as the reverse quench dynamics. Section IV provides a summary of our findings and a number of proposed directions for future study. In Appendix A we present the details of the single-shot procedure, and in Appendix B we show how the quench induced phase separation dynamics is altered for small particle numbers. Finally, in Appendix C we address the convergence of the ML-MCTDHB results.

II. SETUP AND MANY-BODY ANSATZ
To explore the correlated out-of-equilibrium quantum dynamics in a relevant experimental setting, we consider a binary bosonic gas trapped in a 1D harmonic oscillator potential. The MB Hamiltonian consisting of N A , N B bosons with masses m A , m B for the species A, B respectively, reads In the s-wave scattering limit [1] both the intra and interspecies interactions are modeled by a contact poten-tial, where the effective coupling constants are denoted by g AA , g BB , and g AB respectively. Experimentally g σσ can be tuned either via the three-dimensional scattering length with the aid of Feshbach resonances [51,52] or via the corresponding transversal confinement frequency and the resulting confinement-induced resonances [53,54]. Moreover, here we assume that both species possess the same mass, i.e. m A =m B =m, and are confined in the same external potential, i.e. ω A =ω B =Ω. Throughout this work the trapping frequency is fixed to Ω = 0.1 ≈ 2π ×20Hz assuming a transversal confinement ω ⊥ = 2π × 200Hz. Furthermore, we fix the intraspecies interactions to g AA = 1.004 and g BB = 0.9544, which are the values for a binary BEC of 87 Rb atoms prepared in the internal states |F = 1, m F = −1 and |F = 2, m F = 1 [15], while g AB is left to arbitrarily vary upon a quench taking values within the interval g AB = [0, 2]. We remark that in the following the Hamiltonian of Eq. (1) is rescaled in harmonic oscillator units,H = H/( Ω).
Within the MF approximation all particle correlations are neglected. Such a simplification allows for expressing the MB wavefunction of a binary system as a product state of the respective MF wavefunctions where x σ = x σ 1 , . . . , x σ Nσ denote the spatial σ = A, B species coordinates, N σ is the number of σ species atoms and φ σ (x σ i ; t) refers to the time-evolved wavefunction for the σ species within the MF approximation. Employing a variational principle, e.g. the Dirac-Frenkel one [55,56], for the ansatz of Eq. (2) we obtain the corresponding equations of motion in the form of the well-studied system of coupled Gross-Pitaevskii equations [1,2].
The binary BEC is a bipartite composite system residing in the Hilbert space H AB = H A ⊗H B , with H σ being the Hilbert space of the σ species. To incorporate correlations between the different (inter-) or the same (intra-) species, M distinct species functions for each species are introduced obeying M ≤ min(dim(H A ), dim(H B )). Then the MB wavefunction Ψ M B can be expressed according to the truncated Schmidt decomposition [57] of rank M The Schmidt weights λ k (t) in decreasing order are referred to as the natural species populations of the k-th species function Ψ σ k of the σ species. We remark that {Ψ σ k } forms an orthonormal N σbody wavefunction set in a subspace of H σ . To quantify the presence of interspecies correlations or entanglement we use the eigenvalues λ k of the species reduced density matrix ρ Nσ ( , · · · , x σ Nσ−1 ), and σ = σ . When only one (multiple) eigenvalue(s) of ρ Nσ is (are) macroscopic the system is referred to as non-entangled (species entangled or interspecies correlated). It is also evident from Eq. (3) that the system is entangled [58,59] when at least two distinct λ k (t) are finite, further implying that the MB state cannot be expressed as a direct product of two states stemming from H A and H B . In this manner, 1 − λ 1 (t) offers a measure for the degree of the system's entanglement. Moreover, a particular configuration of A species Ψ k ( x A ; t) is accompanied by a particular configuration of B species Ψ k ( x B ; t) and vice versa. Indeed, measuring one of the species states e.g. Ψ A k collapses the wavefunction of the other species to Ψ B k thus manifesting the bipartite entanglement [60,61]. Concluding, the above MB wavefunction ansatz Ψ M B constitutes an expansion in terms of different interspecies modes of entanglement, To include interparticle correlations we further expand each of the species functions Ψ σ k ( x σ ; t) using the permanents of m σ distinct time-dependent single particle functions (SPFs) namely ϕ 1 , . . . , ϕ m σ Here, c k,(n1,...,n m σ ) (t) are the time-dependent expansion coefficients of a particular permanent, P is the permutation operator exchanging the particle configuration within the SPFs, and n i (t) denotes the occupation number of the SPF ϕ i ( x; t). Following the Dirac Frenkel [55,56] variational principle for the generalized ansatz [see Eqs.
(3), (4)] yields the ML-MCTDHB equations of motion [43,44,62]. These consist of a set of M 2 ordinary (linear) differential equations of motion for the coefficients λ k (t), coupled to a set of where σ = σ ,x σ = (x σ 1 , x σ 2 , . . . , x σ Nσ−1 ), and ρ denotes the one-body density matrix of the i-th species function. Note here that the system is termed intraspecies correlated or fragmented if multiple eigenvalues of ρ (1),σ (x, x ) are macroscopically occupied, otherwise is said to be fully coherent or condensed.
The eigenfunctions of the one-body density matrix ρ (1),σ (x, x ) are the so-called natural orbitals φ σ i (x; t). Here we consider them to be normalized to their corresponding eigenvalues, n σ i (natural populations) It can be shown that when the corresponding natural populations obey n σ 1 (t) = N σ , n σ i =1 (t) = 0 and then the first natural orbital φ σ 1 (x σ ; t) reduces to the MF wavefunction φ σ (x σ ; t). Therefore, 1 − n σ 1 (t) serves as a measure of the degree of the σ species fragmentation [63,64].

III. INTERACTION QUENCH DYNAMICS
In the following the quench-induced phase separation dynamics of a binary repulsively interacting BEC is investigated both within and beyond the MF approximation. In particular, interspecies interaction quenches are performed from the miscible to the immiscible regime of interactions and vice versa. Recall [65] that species separation in the absence of a trap occurs for g 2 AB g AA g BB , while the two species overlap when the above inequality is not fulfilled [15]. It is relevant to note, however, that for sufficiently strong trapping-a scenario not considered here-, the above condition is suitably modified [66]. In that case, the g AB needed to induce immiscibility can become substantially larger, as it needs to overcome the restoring, and hence implicitly miscibility favoring, effect of the trap.
First we find the ground state of the system in both the MF and the MB case for fixed intra and interspecies interactions namely g AA = 1.004, g BB = 0.9544, and g AB = 0. To initialize the dynamics we then abruptly vary the interspecies coefficient within the interval g AB = [0, 2], in the dimensionless units adopted herein. Notice that e.g. g AB = 0 corresponds to two decoupled overlapping BECs formed around the center of the harmonic trap. With the above choice of parameters the critical point, i.e. the miscibility-immiscibility threshold, in the absence of the trap, is g AB ≈ 0.9789. The number of particles in each species is fixed to N A = N B = N/2 = 50, with N being the total number of particles of the system. Dynamical phase separation for smaller bosonic ensembles is addressed in Appendix B.

A. Quench dynamics to the immiscible regime
As a first step an interaction quench of an initially species uncorrelated (since g AB = 0) mixture towards the immiscible regime with g AB = 1.2 is performed, driving the system abruptly out-of-equilibrium and letting it dynamically evolve. As shown in Fig. 1, the initial ground state quickly becomes deformed and breaks into multiple filaments within the MF approach, depicted in Figs. 1 (a 1 ) and (a 2 ), as well as in the MB case shown in Figs. 1 (b 1 ) and (b 2 ). The dramatic phase separation observed between the two species, and depicted for t = 60 in the density profiles of Figs. 2 (a 1 ), (a 2 ), results in a different number of filaments formed, the latter being greater within the MF approximation. This suggests that the wavenumber associated with the emergence and growth of these filaments is larger in the MF regime. Notice that in both cases the filaments of the two species locate alternately while the total density does not change dramatically after the filament formation. Additionally here, the first species is found to be expelled further off of the trap center when compared to the second species since this configuration is energetically preferable by virtue of g AA > g BB . Besides the filamentation of its density, each species performs collective oscillations that result in an expansion and contraction of the bosonic cloud. Namely a breathing mode [31,67] possessing a frequency ω br = 2π/T ≈ 0.2 ≡ 2Ω. Finally we remark that for a stronger post quench repulsion, g AB , an increased number of filaments is observed and a more dramatic phase separation takes place, occurring much faster when compared to smaller g AB values.
In all cases, the dominant wavenumber associated with the above-observed unstable dynamics when entering the phase separated regime, is found to be higher in the MF approach when compared to the MB scenario. To quantify the distinct features of the manifestation of the phase separation dynamics within the two approaches we start by considering the stability properties of a homogeneous binary system of length L. Within the MF approximation the spectrum of quasi-particle excitations consists of two branches Ω ± , that in the case of equal masses between the bosons read [68] where n = N/2L denotes the linear atom density [69]. It turns out that if g 2 AB > g AA g BB , i.e. in the immiscible regime of interactions, Ω − becomes imaginary and gives rise to long wavelength modes that grow exponentially in time rendering the homogeneous binary system unstable [70]. For g 2 AB < g AA g BB both branches Ω 2 ± of Eq. (7) remain positive implying that the binary system is stable within this miscible regime. The two species remain then mutually overlapped and undergo a breathing dynamics. Turning to g 2 AB > g AA g BB , the most unstable k = k max modes, corresponding to max{Im(Ω − )}, are presented in Fig. 1 (c 1 ) for varying g AB . For the numerical identification of k max we calculate the spectrum,ρ (1) (k; ω), of the binary system in both the MF and the MB level. Among the modes that appear in this spectrum, we identify as the fastest growing one the mode that maximizes the growth rate ω = ω max . As is evident in Fig. 1 (c 1 ), our numerical findings are in very good agreement with the analytical predictions within the MF approximation (except for very small values of g AB ). Note here, that we have checked the validity of our calculations for different x -20 0 20 trapping frequencies within the local density approximation (see discussion below). However, the unstable modes identified within the MB approach involve considerably shorter k max values which result in longer spatial scales for the filament formation (and thus consist of fewer filaments formed). For example the wavelength obtained in the MF case depicted in Fig. 1 The observed difference of k max between the MF and MB evolution can be attributed to the participation of additional MB excitations which lie beyond the linear response theory as demonstrated, e.g., in [71] for single component setups. Additionally, having identified the wavenumber associated with the fastest growth, we can also infer the time at which the filament formation occurs. We have estimated this time, namely t = t F , by identifying the time at which the amplitude of this wavenumber, k = k max , starts to grow. The formation time, t F , is illustrated in Fig. 1 (c 2 ) for increasing g AB and is fitted by a biexponential function. It is evident that close to the miscibilityimmiscibility threshold (g AB ≈ 1) both approaches coincide, while deviations between the two become apparent as we increase the interspecies interactions. Note also that decreasing the trapping strength towards the homogeneous case alters the time scale at which the instability manifests itself, the more, the closest we are to the above threshold.
To quantify the degree of phase separation we evaluate the overlap integral [72,73] where, Λ(t) = 1 [Λ(t) = 0] denotes complete [zero] overlap of the two species upon abruptly driving the system out-of-equilibrium. As depicted in Fig. 1 (c 3 ) the transition to immiscibility is signalled at slightly earlier times in the MB approach with the overlap between the two species being of about 50% on average, while being almost 60% on average within the MF approximation.
Moreover, the abrupt quench protocol entails rapid oscillations in the MF case when compared to the smoother drop down towards immiscibility observed in the MB scenario. It is worth mentioning at this point, that the same overall phenomenology is observed even upon linearly quenching the system between the same initial and final g AB values (results not shown here for brevity). The key outcome in this case is that the filamentation process is signalled at times proportional to the ramping time used resulting to a larger Λ(t) when compared to the abrupt quench protocol.

B. Single-shot simulations
As a next step we elaborate on how the MB character of the dynamics can be inferred by performing in-situ single-shot absorption measurements [74]. Such measurements probe the spatial configuration of the atoms which is dictated by the MB probability distribution. An experimental image refers to a convolution of the spatial particle configuration with a point spread function. The latter describes the response of the imaging system to a point-like absorber (atom). Relying on the MB wavefunction being available within ML-MCTDHB we mimic the above-mentioned experimental procedure and simulate such single shot images for both species A [namely at each instant of the evolution (for more details see Appendix A) when we consecutively image first the A and then the B species. We remark that the employed point spread function (being related to the experimental resolution), consists of a Gaussian possessing a width w = 1 l ≈ 3.2. Figs. 2 (a 3 ), (a 4 ) illustrate the first and the second simulated in-situ single-shot images at t im = 60 for both species, namely A A (x; t im = 60), and A B (x |A A (x); t im = 60). It is evident that in both shots the two species exhibit a phase separated behavior resembling this way the overall tendency observed in the one-body density [see also Fig. 2 (a 2 )]. However, a direct observation of the one-body density in a single-shot   image is not possible due to the small particle number of the considered binary bosonic gas, N A = N B = 50, as well as the presence of multiple orbitals in the system. The MB state builds upon a superposition of multiple orbitals [see Eqs. (4)-(5)] and therefore imaging an atom alters the MB state of the remaining atoms and hence their one-body density. This is in direct contrast to a MF product state, composed from a single macroscopic orbital, where the imaging of an atom does not affect the distribution of the rest (see also the discussion below for the corresponding variance). Note also here that the above-mentioned single-shot images are reminiscent of the experimental images obtained in a two-dimensional (2D) geometry when examining the phase separation process [13]. To reproduce the one-body density of the system one needs to rely on an average of several singleshot images. Indeed, Fig. 2 (a 5 ) shows within the MB approach the obtained average,ρ (1),σ , over N shots = 1000 images for both species, namelyĀ A (x; respectively. As expected, a direct comparison of this averaging and the actual one-body density obtained within the MB approach [see Figs. 2 (a 2 ) and (a 5 )] reveals that they are almost identical. Finally, let us remark here that similar observations can be made when performing the single-shot procedure initially for the B and then for the A species.
Let us now investigate whether the presence of correlations can be deduced from the time evolution of the variance V(t) of a sample of single-shot measurements [75][76][77]. As before, we mainly focus on the scenario where the imaging is performed first on the A and then on the B species, but the same results can be obtained for the reverse consecutive imaging process. The variance of a concerning the A species reads In the same manner, one defines the variance of a set of single- Figs. 3 (a 1 ), (a 2 ) present both V A (t) and V B (t) with w = 1, and N shots = 1000 at the MF and the MB level respectively. As it can be seen, at the MF approximation V A M F (t) and V B M F (t) remain almost constant exhibiting small amplitude oscillations which essentially resemble the breathing motion that both species feature. However, when inter and intraspecies correlations are taken into account V A M B (t) and V B M B (t) show a completely different behavior. In particular, an increasing tendency is observed at the initial stages of the unstable dynamics, while after the filament formation (t F ≈ 27), V σ M B (t) undergoes large amplitude oscillations reflecting the global breathing of each bosonic cloud. More importantly, the aforementioned increasing tendency of the variance resembles the growth rate of the entanglement, [see 1−λ 1 (t) in Fig. 3 (b 1 )] and the corresponding discussion below]. The above resemblance can be explained as follows. In a perfect condensate, i.e. λ 1 (t) = 1 and n σ is almost constant during the dynamics as all the atoms in the corresponding single-shot measurement are picked from the same SPF ϕ σ (t) [see also Eq. (2)]. The only relevant information that is imprinted in V σ M F (t) concerns the global motion, here the breathing mode, of the entire cloud. It is also worth mentioning here that M F (t) during the MF evolution, testifying the absence of both inter and intraspecies correlations. The observed negligible differences between V A M F (t), and V B M F (t) [hardly visible in Fig. 3 (a 1 )] are caused by the slight deviations in the magnitude of the breathing motion that each species undergoes.
On the contrary, for a MB system where entanglement and fragmentation are present due to the inclusion of inter and intraspecies correlations, the corresponding MB state consists of an admixture of various mutually orthonormal species functions Ψ A k (t) and Ψ B k (t) respectively, k = 1, 2, ..., 15 [see Eq. (3)] each of them building upon different mutually orthonormal SPFs ϕ A i (t) and ϕ B i (t) respectively, i = 1, 2, 3 [see also Eq. (4)]. In this way, the corresponding single-shot variance is drastically altered from its MF counterpart as the atoms are picked from the above-mentioned superposition and thus their distribution in the cloud depends strongly on the position of the already imaged atoms [74,75,78], see also Appendix A. To fairly discern between the impact of the inter and intraspecies correlations on the variance we first inspect V σ (t) when neglecting the entanglement between the species [this approach will be referred in the following as species mean-field approximation (SMF)]. Namely we calculate V σ SM F (t) assuming that the N σ -body state of each species is described by only one species function (Ψ A k (t) = Ψ B k (t)=0 for k = 1) that builds upon distinct SPFs ϕ A i (t) and ϕ B i (t), i = 1, 2, 3. As shown in Fig. 3 (a 2 ) during the filamentation process V σ SM F (t) increases slightly and V A SM . This latter deviation is attributed to the different degree of fragmentation [1 − n σ 1 (t), see e.g. Fig. 3 (b 1 )] that each species possesses after the filamentation process t > 27. Having identified that the presence of fragmentation essentially causes a slight increase on the single-shot variance and more importantly gives rise to deviations between the V σ SM F (t)'s of the two species we can elaborate on the impact of the entanglement when also interspecies correlations are taken into account. In the MB case V σ M B (t) shows a remarkable increasing tendency during the filamentation process highlighting this way the presence of entanglement in the system. Indeed, the increase of entanglement [evident in 1 − λ σ 1 (t)] and consequently of the variance can be attributed to the build up of higherorder superpositions during the filamentation process. Since the absorption imaging destroys the entanglement between the species, we expect that the single-shot images heavily depend on the first few imaged atoms giving rise to pronounced V σ M B (t). We further remark that this increasing tendency of the variance becomes more pronounced (reduced) for larger (smaller) quench val-ues (results not included for brevity). Moreover, during the filamentation process . This latter deviation can be attributed to the different degree of fragmentation that builds up during evolution in each of the two species [compare 1 − n σ 1 (t) for t ≥ 40 illustrated in Fig.  3 (b 1 )]. We finally note that the above-described overall increasing behavior of V A M B (t) and V B M B (t) is robust also for smaller samplings of single-shot measurements, e.g. N shots = 100, or different widths, e.g. w = 0.5, (results not shown here for brevity).

C. Correlation dynamics
The degree of entanglement is encoded in the species functions of the binary system, i.e. Ψ σ k ( x σ ; t), with σ = A, B, being weighted by the λ k (t) coefficients. We remind the reader that if λ 1 (t) = 1 and λ i (t) = 0 (i = 2, ..., k) then the non-entangled limit is reached while if λ k (t) = 0 the more modes are occupied the more strongly entangled the binary system is [58]. In particular, by considering the evolution of the natural occupations λ k (t), depicted in Fig. 3 (b 2 ) it is observed that from the beginning of the quench induced dynamics the occupation of the initial single mode (non-entangled) wavefunction reduces rapidly and higher-lying modes become spontaneously populated. Notice that before the filament formation, e.g. at t ≈ 13, λ 1 ≈ 0.37 and λ 2 ≈ λ 3 ≈ 0.12, while after the breaking (t ≈ 27) the amplitude of the higher-lying modes drops below 0.1 and remains in this ballpark till the end of the propagation. The insets depict selected time instants during the phase separation process of the first three modes of entanglement: namely, just after the breaking [upper insets in Fig. 3 (b 2 )] and the consequent filamentation of the MB wavefunction, and for larger propagation times [lower insets in Fig. 3 (b 2 )] corresponding to Λ(t) ≈ 0.5 during evolution [see also Fig. 1 (c 3 )]. In all cases the leading order mode weighted by λ 1 , and the first two of the higher-lying modes that are predominantly occupied, weighted by λ 2 and λ 3 respectively, are shown for both the A and B species. As it is evident, the dominant mode clearly captures all the filaments formed for both species. The second mode for species A builds a hump at the location centered around the density dip of the first mode, while it also follows the outer filaments formed, and the corresponding third mode mostly supports the inner filaments. As far as the B species is concerned the above observed phenomenology is somewhat reversed. Notice that, the second mode mostly follows the outer filaments, and the third mode is found to be predominantly associated with the filaments developed closer to the trap center.
To further elaborate on the MB nature of the observed quench dynamics we next examine the population of the natural orbitals shown in Figs. 3 (b 3 ), (b 4 ). The occupations of the three natural orbitals used for each of the two species are significant from the early stages of the dynam- ics, with the two lower-lying orbitals being monotonically ordered, acquiring lower populations during evolution.
As already discussed in Sec. II the non-negligible population of both λ k and n σ k (k > 1) signifies the presence of inter-and intraspecies correlations respectively. To identify the degree of intraspecies correlations at the onebody level during the quench dynamics, we employ the normalized spatial first order correlation function [79,80] This quantity measures essentially the proximity of the MB state to a MF (product) state for a fixed set of coordinates x, x . ρ (1),σ (x, x ; t) is the one-body reduced density matrix of the σ species [see also Eq. ). An augmented character of |g (1),σ (x, x ; t)| for increasing distances (e.g. for fixed x ≈ 0, towards x ≈ 25 at t = 7) is also observed. For later evolution times, i.e. after the filamentation process, a significant build up of one-body correlations occurs for both species. Referring to |g (1),A (x, x ; t)|, see Figs. 4 (a 3 ), (a 4 ), we observe that each filament is perfectly coherent with itself (see the diagonal elements), while a small amount of correlations occurs between the inner filaments (|g (1),A (x ≈ 6, x ≈ −6; t = 33)| ≈ 0.9) or the outer ones (|g (1),A (x ≈ 14, x ≈ −14; t = 33)| ≈ 0.8). More importantly, strong correlations appear between neighbouring inner and outer filaments as well as among an inner (outer) filament and its long distance outer (inner) one (|g (1),A (x, x ; t)| ≈ 0.5) signalling their independent nature. Finally, significant losses of coherence are observed between the inner (outer) filaments and the central dip. Turning to |g (1),B (x, x ; t)|, see Figs. 4 (b 3 ), (b 4 ), it is evident that strong correlations appear among each outer and the central filaments (see e.g. x ≈ 10, x ≈ 0 at t = 33) as well as between the outer ones (x = −x ≈ 10 at t = 33). This latter behavior is manifested by the almost vanishing off-diagonal elements of |g (1),B (x, x ; t)| after the filamentation process, indicating a tendency of localization of each filament formed.
is the diagonal two-body reduced density matrix referring to the probability of measuring two particles located at positions x 1 , x 2 at time t.
is the bosonic field operator that creates (annihilates) a σ species boson at position x i . Regarding the same (different) species, i.e. σ = σ (σ = σ ), |g (2),σσ (x 1 , x 2 ; t)| accounts for the intraspecies (interspecies) two-body correlations and is also experimentally accessible via in-situ density density fluctuation measurements [81][82][83]. We remark here that a perfectly condensed MB state leads to |g (2),σσ (x 1 , x 2 ; t)| = 1 and it is termed fully second order coherent or uncorrelated. However, if |g (2),σσ (x 1 , x 2 ; t)| takes values smaller (larger) than unity the state is referred to as anti-correlated (correlated). Let us first comment on the intraspecies two-body correlated character of the dynamics. Focusing on |g (2),AA (x 1 , x 2 ; t)| we observe a consecutive formation of two-body correlations during the dynamics, see Figs. 4 (c 1 )-(c 4 ). Besides a bunching tendency (smaller for the inner filaments) of two bosons to lie within each filament (see the diagonal elements), a correlated behavior is observed among two parity symmetric outer ones (see e.g. x 1 = −x 2 = 14 at t = 33). In addition, an outer filament is anti-correlated both with an inner one (x 1 ≈ 14, x 2 ≈ 6 at t = 33) as well as with the central dip (x 1 ≈ 14, x 2 ≈ 0). Combining this latter behavior with the above suppression of |g (1),A (x, x ; t)| between the filaments, implies the formation of domain-wall-like structures between the area of central filaments and an outer one. Another interesting observation here is that the region between neighbouring inner and outer filaments (e.g. x 1 ≈ 16 at t = 33) is strongly correlated (anti-correlated) with its parity symmetric one. Similar observations can also be made for the |g (2),BB (x 1 , x 2 ; t)|, see Figs. 4 (d 1 )-(d 4 ). Evidently, it is preferable for two bosons to reside within each filament (see the diagonals) or one in each of the outer filaments (e.g. x 1 = −x 2 ≈ 10, t = 33). The central filament is anti-correlated with the outers throughout the dynamics and since |g (1),BB (x, x ; t)| → 0 in the same region, the formation of a domain-wall-like structure between a central and an outer filament can be inferred.
Here, an outer A species filament (x 1 ≈ 14 at t = 33) is anti-correlated (correlated) with the corresponding B species outer located at x 2 ≈ 10 (central at x 2 = 0). However, an inner A species filament (x 1 ≈ 5 at t = 33) is correlated (anti-correlated) with the respective B species outer (central) one. Moreover, we find that the central dip of the A species exhibits a correlated (anti-correlated) behavior with the outer (central) B species filaments. Summarizing the outcome of |g (2),AB (x 1 , x 2 ; t)| is twofold. The fact |g (2),AB (x 1 , x 2 ; t)| = 1 indicates the entangled character of the MB binary system. Additionally, the presence of anti-correlations between the inner and outer filaments of A and B species respectively (or vice versa) supports the phase separation process being imprinted as domain-walls at the two-body level.

D. Reverse quench dynamics
Up to now we explored cases which involve transitions from the miscible to the immiscible phase, by initializing the dynamics from the species uncorrelated (g AB = 0) case and abruptly switching on the interspecies repulsion. Our aim here, is to consider the reverse process, namely initialize the system from a species correlated ground state with g AB = 1.4, i.e. deep in the immiscible regime of interactions, and suddenly reduce g AB . A characteristic example of an immiscible to the immiscible transition with post-quench value g AB = 1.0 is realized in Figs. 5 (a 1 )-(a 4 ). Notice that the phase separated species remain as such at all times with species A forming two humps symmetrically placed around the center of the trap. Closer inspection of the central almost zero density region, suggests that two hardly visible density dips are spontaneously formed in the regions indicated by dashed rectangles in Figs. 5 (a 1 ) and (a 3 ) for the MF and the MB case respectively. These density dips inter-  . We also note that the second and the third natural orbitals of species B are multiplied by a factor of 8 to provide better visibility. (b) Temporal evolution of the overlap integral calculated in both approaches and for both transitions depicted in Fig. 5 (see legend). Other parameters used are the same as in Fig. 1. act with the density peaks created in this species right at their phase boundary, and via this interaction multiple interference fringes can be seen around the center of the trap in both approaches. It is these events which are more pronounced in the MF than in the MB approach, that result in the differences measured in the overlap between the two species. In particular as shown in Fig. 6 (b), Λ M F (t) ≈ 0.35 on average, while Λ M B (t) 0.05 during evolution, which is significantly smaller. The location of these dips is also the location of a "giant" density hump formed in species B. It is also worth mentioning at this point that the evolved phase separated state formed here, consists the 1D analogue of the so-called "ball" and "shell" state that forms in higher-dimensional binary BECs [12].
However a far more rich dynamical behavior of the binary system is observed when the two immiscible species are abruptly quenched towards the miscible regime, with the post-quench value g AB = 0.5. Such a situation is illustrated in Figs. 5 approach. The quench dynamics leads to the formation of multiple DAD solitary waves [48,49] both in the MF and in the MB approach. In the former case, the DAD structures are directly discernible and can be seen to interact and perform oscillations, splitting and recombining within the parabolic trap, in a way remi-niscent of the one-component dark solitons in the experiments of [86,87]. To verify the nature of these structures we further depict as an inset in Fig. 5 (b 1 ) the spatiotemporal evolution of the phase, where the phase jumps corresponding to the location of each dark soliton shown in the density can be easily seen. In contrast to that, in the MB scenario the dynamical evolution of these DAD structures is less transparent, since the system in this case is strongly correlated and the background at which the solitons are formed is highly excited. Recall that dark-bright states are prone to decay in the presence of quantum fluctuations [50] into faster (travelling towards the periphery of the cloud) and slower (remaining closer to the trap center) solitary waves. A similar dynamical phenomenology is also observed here for the abovementioned DAD states. Indeed, at the early stages of the dynamics several decay events occur. Two case examples of such a decay are marked with circles in Figs  reach at different times the periphery of the cloud and thus multiple collision events occur at different times during evolution. A case example of such a collision is indicated with arrows in Figs. 5 (b 3 ), (b 4 ).
To expose the multi-orbital nature of the above dynamics, both the one-body density as well as the different orbital contributions are depicted in Figs. 6 (a 1 )-(a 6 ) at initial (t = 15), intermediate (t = 27) and larger evolution times (t = 40). Notice that at initial times the two species are still phase separated, while the first orbital predominantly describes the MB dynamics of the system. Here, we can easily measure the number of DAD solitary waves that are initially formed, illustrated with two-directional arrows in Figs. 6 (a 1 ) and (a 4 ), by observing that each density dip created in species A, Fig. 6 (a 1 ), is filled by a density hump (on top of the BEC background) developed in species B, Fig. 6 (a 4 ), and vice versa. Furthermore, it is found that consecutive orbitals within the same species also follow the abovedescribed phenomenology with a clearly visible domainwall [4,45] formed between the second and the third orbital of species B [see arrows in green in Figs. 6 (a 4 )-(a 6 )]. For intermediate times the merging of the most inner solitary states discussed above is indicated with circles in Figs. 6 (a 2 ), (a 5 ). Notice the pronounced density hump that occurs in species B around the center of the trap, being supported by all three orbitals developed in this species. Additionally, also the faster DAD solitary waves are monitored in this time slice, where again it is observed that these states are supported by all orbitals used in each of the two species being marked with dashed rectangles. However, at larger propagation times and since we "kicked" the system towards miscibility, multiple interference events more pronounced in species B, result to a dephasing of these matter wave patterns and most of these states are lost as can be seen in Figs. 6 (a 3 ), (a 6 ), rendering the two species mostly overlapped. Notice the increasing tendency towards miscibility with the overlap integral [see again here Fig. 6 (b) for g AB = 0.5] reaching its maximum value, Λ M B (t ≥ 60) ≈ 0.95, at large propagation times, when compared to the MF approximation. In the latter case, Λ M F (t) ≈ 0.85 is reached from the early stages of the dynamics remaining on average almost the same as time progresses.
To conclude our investigation, let us also briefly com-ment on the manifestation of the MB correlated character of the quench-induced dynamics with the aid of in-situ single-shot measurements. Figs. 7 (a), (b) present the first and the second simulated in-situ single-shot images at t im = 15 for both species, with the DAD structures being clearly imprinted in both shots. Notice that the two species are almost completely overlapped resembling the overall tendency observed in the averaged, over N shots = 1000, one-body density illustrated in Fig. 7 (c). By inspecting the corresponding variances [see also Eqs. (9) and (10)] during the evolution shown in Fig. 7 (d), we observe that within the MF We should bear in mind that the initial pre-quenched state is both strongly fragmented and entangled on the MB level. Therefore, in this strongly correlated scenario both fragmentation as well as entanglement are greatly manifested in the evolution of the variance of a set of single-shot measurements.

IV. CONCLUSIONS
In the present work we explored the quench-induced phase separation dynamics of an inhomogeneous repulsively interacting binary BEC both within and beyond the MF approximation including multiple orbitals. To achieve such a miscible to immiscible transition (positive quench case) the intraspecies interactions are held fixed and the system is abruptly driven out-of-equilibrium by switching on the interspecies repulsion. Quench dynamics leads to the filamentation of the density of each of the two species and also in both approaches (MF and MB) while the filaments formed perform collective oscillations of the breathing-type. The wavenumbers associated with the observed growth are identified to be shorter in the MB case for all g AB values that we have checked, whilst our numerical findings at the MF level are in very good agreement with the analytical predictions available in this limit, as regards the instability growth rate. It is found that increasing the interspecies repulsion, not only accelerates the filamentation process but also increases the number of filaments formed in both approaches, occurring faster on the MB level. Additionally, stronger interspecies repulsion leads to almost complete phase separation being more pronounced in the MB scenario. We further note, that upon fixing the interspecies repulsion while decreasing significantly the system size (few boson case) phase separation is absent in the MB case while still present at the MF limit.
Detailed correlation analysis at the one-and the twobody level bear the signature of the phase separation process as the miscibility-immiscibility threshold is crossed. On the one-body level significant losses of coherence are observed, verifying the fragmented nature of the system, between filaments residing around the center of the trap with the longer distant ones lying at the periphery of the bosonic cloud. At the two-body level domain-wall-like structures are revealed, since the inner filaments in both species are found to be anti-correlated with their respective outer ones. These domain-walls support the fact that for smaller interspecies interactions, but well inside the immiscible regime, we never observe perfect de-mixing of the two species. Furthermore, and even more importantly, the presence of both entanglement and fragmentation are related to the variance of single-shot images, that are utilized for the first time in the current effort for binary systems, offering a direct way for the experimental realization of the observed dynamics. In particular, it is found that the growth rate of the variance resembles the growth rate of the entanglement. The fragmentation of the binary system is captured by the deviations in the variance measured in the course of the dynamics with respect to each of the two species.
Interestingly enough, when considering the reverse (negative) quench scenario, namely quenching from the immiscible towards the miscible regime multiple darkantidark solitary waves are spontaneously generated in both approaches and they are found to decay in the MB case [50]. The evolution of the variance of single-shot measurements reveals enhanced entanglement, since the system in this case is strongly correlated on the MB level. Finally, for transitions inside the immiscible regime we retrieve the 1D analogue of the so-called "ball" and "shell" structure that appears in higher-dimensional binary BECs [12,88].
There are multiple directions that are of interest for future work along the lines of the current effort. A systematic study of the dynamical phase separation process following a time-dependent protocol (e.g. a linear quench) presents one of the major computational chal-lenges for further study. In particular, in such a scenario one can explore the domain formation crossing the critical point with different velocities and thus testing the Kibble-Zurek mechanism [69] in the presence of quantum fluctuations. However, to examine the latter, a major challenge that it is imperative to overcome is that of considering low atom numbers, in order to explore the associated thermodynamic limit, avoiding the potential influence of finite size effects. Another straight forward direction is to consider the corresponding already experimentally realized [13] 2D setting, and examine how the MF properties are altered in the presence of quantum fluctuations. Also of great interest would be to consider the quench dynamics of spinor BECs, for which phase separation processes are of ongoing interest at the MF limit [89] and also investigate the relevant MB aspects. As in the single component case, the single-shot simulation procedure relies on a sampling of the MB probability distribution [74,75,78]. The latter is available within the ML-MCTDHB framework. However, in a twospecies BEC and when inter and intraspecies correlations are taken into account, the entire single-shot procedure is significantly altered when compared to the single component case. Here, the role of entanglement between the species manifested by the Schmidt decomposition [see Eq.
(3)] and in particular the Schmidt coefficients λ k 's play a crucial role concerning the image ordering.
For instance, to image first the A and then the B species we consecutively annihilate all the N A particles. Focusing first on a certain imaging time instant, t im , a random position is drawn according to the constraint ρ (1) N A (x 1 ) > l 1 where l 1 refers to a random number within the interval [0, max{ρ (1) N A (x; t im )}]. Then we project the (N A + N B )-body wavefunction to the (N A − 1 + N B )body one, by employing the operator 1 N (Ψ A (x 1 ) ⊗Î B ), whereΨ A (x 1 ) denotes the bosonic field operator that annihilates an A species boson at position x 1 and N is the normalization constant. The latter process directly affects the λ k 's (entanglement weights) and thus despite the fact that the B species has not been imaged yet, both ρ (1) This can be easily understood by employing once more the Schmidt decomposition. Indeed after this first measurement the MB wavefunction reads where the normalization factor andλ i, are the Schmidt coefficients that refer to the (N A − 1 + N B )-body wavefunction. The above-mentioned procedure is repeated for N A − 1 steps and the resulting distribution of positions (x 1 , x 2 ,...,x N A −1 ) is convoluted with a point spread function leading to a single-shot for the A species. Herex refers to the spatial coordinates within the image and w is the width of the point spread function. It is worth mentioning also at this point that before annihilating the last of the N A particles, the MB wavefunction has the form where |Φ A i,1 (t im ) denotes a single particle wavefunction characterizing the A species. Then, it can be easily shown that annihilating the last A species particle the MB wavefunction reads where x|Φ A j,1 is the single particle orbital of the j-th mode. After this last step the entanglement between the species has been destroyed and the wavefunction of the B species |Ψ N B M B (t im ) corresponds to the second term of the cross product on the right hand side of Eq. (A3). In this way, it becomes evident that |Ψ N B M B (t im ) obtained after the annihilation of all N A atoms is a non-entangled N Bparticle MB wavefunction and its corresponding singleshot procedure is the same as in the single species case [74]. The latter is well-established (for details see [74,75]) and therefore it is only briefly outlined below. Referring to t = t im we first calculate ρ N B (x; t im )]. Next, one particle located at a position x 1 is annihilated and ρ (1) . Following this procedure for N B − 1 steps we obtain the distribution of positions (x 1 , x 2 ,...,x N B −1 ) which is then convolved with a point spread function resulting in a single-shot A B (x |A A (x)).
We remark here that the same overall procedure can be followed in order first to image the B and then the A species. Such an imaging process results in the corresponding single-shots A B (x) and A A (x |A B (x)). ditionally, the interparticle repulsion between the species leads to breathing-type oscillations of the particle densities.
As the number of particles is increased, namely for N = 40, the one-body density evolution of the A species shown in Figs. 8 (e), (g) for the MF and the MB scenario respectively also differ. In particular, while in both approaches four filaments are formed, they are found to be significantly broader in the MB case. This broadening together with the breathing that the cloud undergoes, leads to an attraction, collision, and repulsion of the inner filaments in a periodic manner, being more pronounced in the MB case when compared to the single merging, and repulsion observed at around t ≈ 70 in the MF approach of Fig. 8 (e). Moreover, the disparity between the two approaches becomes rather transparent when further inspecting the spatio-temporal evolution of the density of species B illustrated in Figs. 8 (f ), (h) for the MF and the MB case respectively. Interestingly here, in the MB scenario only two filaments are formed located alternately in regions that correspond to density dips of species A, restoring the phase separation process absent for smaller particle numbers. However, the central filament created in the MF approach [see for comparison Fig. 8 (f )] is clearly absent in the MB case, resulting in this way in a larger overlap between the two gases at the MB level.

Appendix C: Remarks on Convergence
Let us first briefly comment on the main features of our computational methodology, ML-MCTDHB, and then showcase the convergence of our results. ML-MCTDHB [43,44] constitutes a flexible variational method for solving the time-dependent MB Schrödinger equation of bosonic mixtures. It relies on expanding the total MB wavefunction with respect to a time-dependent and variationally optimized basis, which enables us to capture the important correlation effects using a computationally feasible basis size. Finally, its multi-layer ansatz for the total wavefunction allows us to account for intraand interspecies correlations when simulating the dynamics of bipartite systems. For our simulations, we use a primitive basis consisting of a sine discrete variable representation containing 800 grid points. To perform the simulations into a finite spatial region, we impose hard-wall boundary conditions at the positions x = ±50. Note that the Thomas-Fermi radius of each bosonic cloud is of the order of 20 and we never observe appreciable densities beyond x = ±30. Therefore the location of the imposed boundary conditions is inconsequential for our simulations. The truncation of the total system's Hilbert space, namely the order of the considered approximation, is indicated by the used numerical configuration space C = (M ; m A ; m B ). Here, M = M A = M B refers to the number of species functions and m A , m B denote the amount of SPFs for each of the species. In the limit M = m A = m B = 1 the ML-MCTDHB ex-pansion reduces to the MF ansatz. Finally, in order to guarantee the accurate performance of the numerical integration for the ML-MCTDHB equations of motion the following overlap criteria | Ψ|Ψ − 1| < 10 −10 and | ϕ i |ϕ j − δ ij | < 10 −10 have been imposed for the total wavefunction and the SPFs respectively.
Next, we demonstrate the order of convergence of our results and thus the level of our MB truncation scheme. To show that our MB results (more specifically the quantities and observables considered here) are numerically converged, we inspect for the σ species the overlap ∆ σ CC (t) = 1 − δ σ CC (t) between the one-body densities ρ evolution times. In the same manner, convergence occurs for a varying number of SPFs in both species. For instance, ∆ A CC (t) [∆ B CC (t)] between C = (20; 3; 3) and C = (20; 4; 4) shows a maximum deviation of the order of 7% for large evolution times. Similar observations can be deduced also for the case of even smaller particle numbers, and the reverse quench scenario (not included here for brevity reasons). To summarize, according to the above systematic investigations, the considered orbital configurations provide adequate approximations for the description of the non-equilibrium correlated dynamics.