Dark-Bright Soliton Dynamics Beyond the Mean-Field Approximation

The dynamics of dark-bright solitons beyond the mean-field approximation is investigated. We first examine the case of a single dark-bright soliton and its oscillations within a parabolic trap. Subsequently, we move to the setting of collisions, comparing the mean-field approximation to that involving multiple orbitals in both the dark and the bright component. Fragmentation is present and significantly affects the dynamics, especially in the case of slower solitons and in that of lower atom numbers. It is shown that the presence of fragmentation allows for bipartite entanglement between the distinguishable species. Most importantly the interplay between fragmentation and entanglement leads to the splitting of each of the parent mean-field dark-bright solitons, placed off-center within the parabolic trap, into a fast and a slow daughter solitary wave. The latter process is in direct contrast to the predictions of the mean-field approximation. A variety of excitations including dark-bright solitons in multiple (concurrently populated) orbitals is observed. Dark-antidark states and domain-wall-bright soliton complexes can also be observed to arise spontaneously in the beyond mean-field dynamics.


I. INTRODUCTION
Over the last two decades, atomic Bose-Einstein condensates (BECs) have offered an ideal testbed for the exploration of coherent structures relevant to nonlinear wave systems [1,2]. Most of these entities have not only been theoretically predicted but also in many case examples experimentally verified [3][4][5]. Some of the most notable ones are the bright [6][7][8], dark [9] and gap [10] matter-wave solitons. Higher dimensional analogues of these structures have also been examined, such as vortices [11,12], solitonic vortices and vortex rings [13].
Although the main emphasis of study has been the case of single component systems, in recent years, the study of multi-component variants has also gained a significant traction [2]. There, new fundamental waveforms may arise such as the symbiotic dark-bright (DB) soliton in self-defocusing media (where a bright wave would not exist, yet it survives being confined by the dark component) [14]. Interestingly, the experimental study of such states was initiated considerably earlier in nonlinear optics, per the observation of DB solitary wave structures and their molecules [15,16]. More recently, atomic BECs enabled a wide variety of relevant studies initially motivated by the theoretical study of [17]. In particular, the experimental realization of DBs [18] was followed by a series of experiments exploring the dynamics and properties of these states including in-trap oscillations of the solitons, their spontaneous generation (e.g. via counterflow experiments) and their interactions both with other DBs and with external potential barriers [19][20][21][22][23][24].
A natural question that has emerged in connection with solitary waves concerns the fate of these types of excitations in the presence of quantum fluctuations [25][26][27][28]. In atomic BECs, this has been considered for bright solitons [29][30][31][32][33][34][35][36][37][38], and dark solitons [39][40][41][42][43], in trapped settings, as well as in the presence of optical lattices [44,45]. A number of these works, including [46] among others, has also explored the role of quantum fluctuations and higher orbitals in the presence of solitonic interactions. Many of the resulting findings suggest that slower solitons and smaller atom numbers result in significant deviations from the mean-field, Gross-Pitaevskii (GP) limit.
However, to the best of our knowledge, such studies have not been performed in a systematic fashion in multi-component settings and for associated solitary wave structures. In that light, herein we explore the case of DB solitons and their dynamics, as well as collisions both at and beyond the mean-field limit. To incorporate the quantum fluctuations stemming from the correlations in the DB soliton dynamics, we employ the Multi-Layer Multi-Configuration Time-Dependent Hartree Method for bosons (ML-MCTDHB) [47,48] designed for simulating the quantum dynamics of bosonic mixtures. We consider both the oscillation of a single DB solitary wave in a trap, as well as the interaction of two symmetric solitary waves inside a parabolic trap. We compare and contrast the findings of the mean-field case (where a single orbital is effectively used in the dark-and bright-components) with cases where multiple orbitals are used. In all cases it is found that the initial mean-field DB solitons split into daughter DB solitary waves, in contrast to the well-known mean-field predictions. The robustness of the presented results is ensured by exploring settings that involve a higher number of orbitals thereby supporting the validity of our approximation and of the observed beyond mean-field excitations (see also Appendix B). Within the employed multi-orbital approximation bipartite entangle-ment (see [49,50] and references therein) between the distinguishable species resulting from the spontaneous fragmentation of the DBs is generally present. More importantly it is the interplay between fragmentation and the resulting entanglement which gives rise to the observed dynamical structures.
The dynamics of a single DB soliton being initialized off-center in the parabolic trap, is first analyzed. It is shown that from the early stages of the dynamics the initial DB solitary wave becomes both fragmented and entangled, and consequently splits into two DB solitonic fragments signalling its decay. This is in direct contrast with the pure mean-field description where the DB soliton remains intact for large propagation times. After the splitting process novel structures of multi-orbital nature emerge, such as dark-antidark states [51,52], with the latter corresponding to density humps on top of the BEC background, and domain-walls (DW) [1,2]. The corresponding decay time increases for larger initial velocities thus approaching the mean-field limit. Furthermore the lifetime of the DB entity decreases for larger particle number imbalance (i.e. varying the number of atoms in the dark soliton component upon keeping fixed the number of atoms in the respective bright counterpart). As a next step the collisional dynamics between two DBs is investigated. We find that in this case too, both fragmentation and entanglement are evident from the beginning of the dynamics. Perhaps more intriguingly, a variety of unprecedented (and some unique to the multi-component setting) states are identified, even if in a transient form during our dynamical evolutions. These include among others DB solitons in higher orbitals alone, DB solitons involving both the lower and higher orbitals, DWs, as well as dark-antidark type structures between the orbitals of the same (dark) component. It is important to note that in the cases reported significant fractions of populations arise in higher orbitals even among the 500 or 1000 atoms present in the dark component.
Our presentation proceeds as follows. Section II provides an overview of the setup and background, in which we discuss the preparation of the DB solitons and related states. Next, we present and discuss our numerical results for the single DB soliton dynamics (section III) and the DB soliton collisions (section IV). Finally, in section V, we provide our conclusions. Appendix A focuses on the initial state preparation, while Appendix B is dedicated to a discussion of our computational approach and the convergence of our simulations.

II. SETUP AND BACKGROUND KNOWLEDGE
Dark-bright solitons are non-linear excitations observed in one-dimensional (1D) population-imbalanced binary BECs. Such excitations are characterized by a density depletion for the majority component (that we will hereafter dub species D) of the BEC, accompanied by a density accumulation for the minority component (hereafter referred to as species B). In the mean-field approximation the prototypical 1D model where such states can be found to arise [21], is a system of coupled GP equations, i.e. a vector variant of the nonlinear Schrödinger equation [53][54][55] with cubic nonlinearity. Each DB soliton is characterized by its position x 0 j , its inverse width d j , and by the so-called soliton's phase angle a j . The latter is associated with the soliton's velocity u j /c = sin a j , with c = gn/m, being the speed of sound. Here, g denotes the interparticle interaction, n refers to the local particle density, and m is the particle mass.
In the following we assume N S DB solitons (with N S = 1 and N S = 2 corresponding to a single and two DB solitons respectively) embedded in a background density,φ 0 (x). The wavefunction of the majority (dark) component [56] reads The corresponding wavefunction for the minority (bright) component and upon considering in-phase bright counterparts reads Note that, µ is the chemical potential of the bright component, being fixed so as to correspond to a total number of atoms for the bright counterpart N B = 5, that we also keep fixed throughout this work. Furthermore, η j refers to the amplitude of the bright soliton. In the absence of a confining potential and for N S = 1, this amplitude is related to the amplitude of the dark soliton, √ µ cos a j , where µ is the background chemical potential, and the inverse width d j by The above expressions namely Eqs. (1)-(2), represent an approximate initial profile of the mean-field (i.e. single orbital) setting. The approximate nature of the profile stems from the effective multiplication with the equilibrium backgroundφ 0 (x) (at least in the case whereφ 0 (x) is not a constant), which in the Thomas-Fermi limit reads:φ 0 (x) = 1/ √ g µ − V (x). Here, V (x) refers to the external trapping potential. We note also here, that in the case of collisions to be presented below, the two DB solitons (N S = 2) should be far enough apart to be individual entities, i.e., their width (proportional to 1/d j ) should be much smaller than their relative separation |x i − x j |.
In addition to the above approximation, the realm of the mean-field itself leads to the following two assumptions (irrespectively of configuration): (a) the two species of the binary BEC are uncorrelated and (b) the constituting particles of each component are uncorrelated. Therefore, the total many-body wavefunction within the mean-field approximation is expressed in terms of the mean-field wavefunctions where N D (N B ) refers to the number of atoms for the D (B) species, while ; t) denotes the time-evolved wavefunction for the species D (B) respectively within the mean-field approximation. The equations of motion for the mean-field ansatz [see Eq. (4)] yield the well-studied GP equation.
The binary BEC consisting of species D, B (with Hilbert spaces H D , H B respectively) is a bipartite composite system with the corresponding Hilbert space H DB = H D ⊗ H B . To examine the system of N S solitons beyond the mean-field approximation, we incorporate the inter-and intra-species correlations by introducing M distinct species functions for each component of the binary BEC. Then the many-body wavefunction Ψ M B can be expressed according to the Schmidt decomposition where the coefficient λ k (t) is referred to as the natural occupations of the species function k. We remark that due to M < min(dim(H D ), dim(H B )) the above expansion (see Eq. (5)) corresponds to a truncated Schmidt decomposition of rank M [49]. It is also important to mention that a state of the bipartite system (see Eq. (5)) cannot be expressed as a direct product of two states from the two subsystem Hilbert spaces H D , H B if at least two coefficients λ k (t) are nonzero. In the latter case the system is referred to as entangled [57,58]. Note that a particular particle configuration of species D (represented by Ψ k ( x D ; t)) is accompanied by a particular particle configuration of species B (denoted by Ψ( x B ; t)) and vice-versa [59]. A corresponding measurement of one of the species states e.g. Ψ D k collapses the wavefunction of the other species to Ψ B k thus manifesting the bipartite entanglement [60,61]. In this way, the above ansatz constitutes an expansion for the many-body wavefunction Ψ M B in terms of inter-species modes of entanglement. In the following we shall refer to λ k (t)Ψ D k ( x D ; t)Ψ B k ( x B ; t) as the k-th mode of entanglement. We remark that in the case of λ 1 (t 0 ) = λ 2 (t 0 ) = 1/2 the many-body state of the system Ψ M B x D , x B ; t 0 at time t 0 is referred to as a Schrödinger cat state. The particle correlations are included by constructing each of the species functions Ψ σ k ( x σ ; t) using the permanents of m σ distinct time-de-pendent single particle functions (SPFs, ϕ 1 , . . . , ϕ m σ ) where P is the permutation operator exchanging the particle configuration within the SPFs, c k,(n1,...,n m σ ) (t) denote the time-dependent expansion coefficients of a particular permanent. N σ refers to the particle number within species σ and n i (t) is the occupation number of the SPF ϕ i ( x; t). Therefore the N σ -body Hilbert space of species σ is approximated by H σ ∼ S H ⊗Nσ ϕ σ (t) , where S is the symmetrization operator and H ϕ σ (t) the single particle Hilbert space spanned by the SPFs. Following a time-dependent variational principle, e.g. the Dirac Frenkel [62,63] or the McLachlan variational principle [64] for the generalized ansatz [see Eqs. ; t) and m D + m B nonlinear integrodifferential equations for the SPFs ϕ k ( x; t). To the best of our knowledge analytical solutions of the many-body ansatz [Eqs. (5), (6)] that contain DB solitons are not known while systematic numerical studies in this direction are still lacking. Here, we utilize the many-body approach that ML-MCTDHB provides and embed at t = 0 the mean-field wavefunction [see Eq. (4)] within the many-body ansatz [see Eqs. (5), (6)]. To achieve the latter, we consider λ 1 (0) = 1, λ i =1 (0) = 0 for the natural occupations of the species functions, c 1,n1=N D (0) = 1, (1), (2)] for the expansion coefficients of species D (B) and ϕ D 1 (x; 0) = 1 for the SPFs of the majority (minority) component. The mean-field ground-state is used as the background density,φ 0 (x). Summarizing, we initialize the many-body quantum dynamics problem employing the mean-field initial state, aiming to examine how the single-orbital population will spontaneously give rise to higher orbital dynamics.
For binary mixtures the one-body density can be expanded in different modes stemming from the species functions expansion [see Eq. (5)]. For instance, the one-body density of the majority species D reads where ρ (1),D i (x, x ; t) refers to the one-body density matrix of the i-th species function andx D = The one-body density for the minority species ρ (1),B (x, x ; t) can be defined in a similar manner. The natural orbitals, φ σ i (x; t), are defined as the eigenfunctions of the one-body density matrix ρ (1),σ (x, x ). For our purposes we consider them to be normalized to their corresponding eigenvalues, n σ i (natural occupations) as In this way, we can ensure that when our many-body ) the corresponding natural occupations obey n σ 1 (t) = N σ , n σ i =1 (t) = 0. In the latter case the first natural orbital φ σ 1 (x σ ; t) reduces to the mean-field wavefunction φ σ (x σ ; t). To examine the beyond-mean-field dynamics of a DB soliton in a setting relevant to most of the recent experiments, we consider a binary bosonic gas trapped in a 1D harmonic oscillator potential. The many-body Hamiltonian consisting of N D , N B bosons with mass m D , m B for the species D, B respectively, reads Within the ultracold s-wave scattering limit [53], we model both the inter-and intra-species with contact interaction by a delta function with respect to the relative coordinate of two bosons and the strength being denoted by g DB , g DD , g BB respectively. We also assume that the bosons of both species possess the same mass, i.e. m D =m B =m and experience the same external potential, i.e. ω D =ω B =Ω, and effective 1D coupling strengths, i.e. g DD =g BB =g DB =g. We remark that the effective 1D coupling strength [66] becomes denotes the Riemann zeta function at x = 1/2. Here, the transversal length scale is a ⊥ = /mω ⊥ , with ω ⊥ the frequency of the transversal confinement, while a σσ s denotes the free space s-wave scattering length within or between the two species. Recall at this point, that the miscibility/immiscibility condition in the absence of a trap reads a 12 s ≤ a 11 s a 22 s [67]. The latter refers to the absence or presence of phase separation between the species, and our choice of equal coupling strengths corresponds to the miscibility/immiscibility threshold in the above expression. As a first approximation and motivated by the very nature of the DB state (which is somewhat intrinsically phase separated as one component operates as a well for the atoms of the other) and by the insensitivity of the DB features near this threshold in the mean-field limit [68], we operate at the Manakov limit of equal values of all a σσ s 's. The interaction strength can be tuned via a σσ s with the aid of Feshbach resonances [69,70]. Finally, note that a scaling transformation has been performed in Eq. (8) setting the length scale equal to the longitudinal characteristic oscillator length a = /mω , the energy scale to ω and the scaled interaction strength is g = g 1D m/ 3 ω . The assumption of equal interactions, as also discussed in [11,14], is tantamount to considering, e.g., the case of hyperfine states of 87 Rb; in this case also the masses are equal. The coefficients of interactions are not exactly equal, but are very proximal to that limit [19,21] (see also the more recent work of Ref. [71]). Finally, for reasons of computational convenience, we shall set = m = g = 1, and therefore all quantities below are given in dimensionless units.
To initialize the beyond mean-field dynamics we first trace the mean-field ground stateφ 0 (x), using Newton's method. The latter is applied to the following steadystate problem: On top of this relaxed mean-field background for fixed number of atoms N D , we then embed the N S solitons of Eqs. (1)-(2) at t = 0 with fixed amplitude as well as number of atoms N B for the bright component. Then the inverse widths d i are obtained using Eq. (3). We also note that the corresponding free parameters are the velocity u i and the position x 0 j of the DB solitons. We remark that by following the above-mentioned procedure we minimize the sound wave emission during the dynamics. For more details on the selection of the soliton and background density parameters we refer the reader to Appendix A.
Given our interest in the qualitative characteristics of the emerging structures, we study their persistence with increasing number of allowed modes of entanglement [specified by the number of species functions M in the Eq. (5)] for fixed number of single particle configurations [specified by the number of SPFs for species D, m D and Note that the 1-(1, 1) approximation corresponds to the mean field case [see also Eqs. (4), (5), (6)]. For further discussion on convergence issues we refer the reader to Appendix B.

III. SINGLE DB SOLITON DYNAMICS
First we explore the oscillation of a DB solitary wave inside the parabolic trap. We consider a binary bosonic gas consisting of N D = 300 atoms and N B = 5 atoms confined in a 1D harmonic trap with frequency Ω = 0.1. The dynamics is initialized with a single DB soliton with velocity u 1 /c = 0.5 and initial position x 0 = −2.5. The chemical potential of the background density, the amplitude and the inverse width of the DB soliton are chosen as µ = 6.42, η = 1.88 and d = 1.42, respectively (for details see Appendix A).
In the mean-field case the numerically obtained (quarter-)period of oscillation for the single DB soliton is measured as T osc /4 ≈ 32.5 (see Fig. 1 (a), (b) and the corresponding insets). This result is in good agreement with the analytical predictions of Ref. [17] (see also Ref. [21]) which for the parameters used herein reads T osc /4 ≈ 30. In contrast to that, in the correlated 15-(2,4) case, shown in panels (c) and (d) of Fig. 1, the single DB soliton decays (t ≈ 5) into a faster, and a slower DB soliton (alias soliton fragments). By inspecting the de- To gain more insight into the decay dynamics of a single DB soliton we next study the population of the natural orbitals. We remind the reader that a state with n i (t) = 1 is referred to as fully condensed, while for n i (t) = 1 fragmentation phenomena arise [72,73]. The occupations of the two natural orbitals used for species D   Fig. 3 (e)]. After this initial relaxation stage, the splitting of the dark soliton leads to the emergence of two solitons manifested by density dips appearing in the first orbital [see also solid ciel line at t 2 = 7 and around the center of the trap in Fig. 3 (f )]. Interestingly enough, within the cloud and in particular around x = −20 for evolution times up to t = 3, regions with additional density depletion for the first orbital and a corresponding density accumulation in the second orbital can also be observed, indicated by dashed blue circles in Fig. 3(a) and (c), respectively. These are reminiscent of the recently experimentally observed dark-antidark states [51].  [51] appear, e.g., in the same panel and around x = −20, or at later times and e.g., around x = 10 in Fig. 3 (f ). In both of the above cases the location of the formation of these matter wave entities is denoted by dashed black circles. Furthermore, in this latter panel, Fig. 3  around x = 3 a multi-orbital DB entity is formed within the region indicated by black arrows. The location of the formation of this structure corresponds also to the location where the two orbitals of species B both being bright solitons are aligned. It is also important to comment here that the two bright soliton fragments of the first orbital (denoted by solid blue lines) are in-phase while the respective bright solitons formed in the second orbital of the minority species are out-of-phase (see nodal structures denoted by solid gray lines in the same panel). Intriguingly, at later times and as far as the slow solitons are concerned, the initial dark-antidark structure evolves into a DW between the first and the second orbital of the majority species D [see around the center of the trap for t 3 = 11 in Fig. 3 (g) also indicated by a black arrow]. The formation of this DW corresponds to the location where species B localizes both of its orbitals at this time. Hence, an unusual domain-wall-bright (DWB) soliton pattern appears to arise. Notice that at the same time, the fast DB structure is also visible around x ≈ 9 (indicated by black arrow). For larger propagation times the aforementioned DWB structure remains robust [see top and bottom panels of Fig. 3 (h)], while the multiorbital fast DB structure is highly pronounced as indicated by a black arrow in the same figure.

(f ), and in particular
Let us now elaborate in more detail on the decay mechanism. As stated above, see Figs. 1 (e)-(h), the splitting dynamics, i.e. the decay of the initial mean-field soliton into two distinct solitary waves, is more transparent on the level of the species functions. It is already known [39][40][41][42][43] that a single mean-field dark soliton, when exposed to quantum fluctuations, decays due to the filling of the initially density depleted region by localized states. In our case, similar dynamics can be observed in the density evolution of the first and second natural orbital of the dark component, see Fig. 3 (e). Here, the filling process progresses very rapidly as Fig. 2 (b) suggests, in particular see the curvature of n 1 (t), n 2 (t) at the initial time instants. Next, we discuss how the above-mentioned established dynamics triggers the splitting mechanism in our multispecies case by inspecting the species function dynamics being directly related to the natural orbitals, see also Eq. (6). At times before the manifestation of the splitting the first two species functions Ψ D i (t) acquire similar occupations. However due to their mutual orthogonality they possess different weights [i.e. c 1,(n1,N −n1) (t) = c 2,(n1,N −n1) (t) in Eq. (6)] on the two occupied natural orbitals. The latter is manifested in the upper panel of Fig. 3 (i) as a shift in the position of the corresponding density minima of ρ (t). Thus, the decay mechanism is also manifested in the bright component. Summarizing, the splitting mechanism is a manifestation of the dynamical build up of the first two species functions due to the interspecies coupling and it is triggered by the filling of the first orbital by the second one in the dark component.
To infer about the degree of first order coherence during the DB soliton dynamics, we employ the normalized spatial first order correlation function [74] Here, ρ (1),σ (x; t) refers to the one-body density matrix of σ species at position x. We also note that |g (1),σ (x, x ; t)| 2 is bounded taking values within the range [0, 1]. Indeed, a spatial region with |g (1),σ (x, x ; t)| 2 = 0 is referred to as perfectly incoherent, while if |g (1),σ (x, x ; t)| 2 = 1 is said to be fully coherent. Especially, in our case we are interested for the appearance of Mott-like correlations. By this we mean the localization of the one-body correlations within a certain spatial region R (i.e. |g (1),σ (x, x ; t)| 2 ≈ 1 x, x ∈ R) and perfect incoherence between different spatial regions R, R (i.e. |g (1) (1),D (x, x ; t)| 2 for different time instants during the dynamics, namely before and after the splitting of the initial mean-field DB soliton. As it can be seen, there is a gradual formation of Mott-like correlations, in particular after the decay, being separated by the location where the DWB develops. Similar type of correlations is also observed for the bright component, see the insets in Figs. 3 (m)-(p). Indeed, a localization of the |g (1),B (x, x ; t)| 2 is manifested within the spatial region in which the density of the DWB is dominant.
In order to examine the validity of a form of the particle picture (which is known to hold in the mean-field limit) in terms of the above-mentioned solitary fragments of the initially imprinted DB soliton, we employ the position expectation value of the first (corresponding to the fast DB structure) and the second (attributed to the slow DB structure respectively) species function that reads Here, t) refers to the particle number and ρ (1),B i (x; t) denotes the one-body density of the i-th B species function. We shall show that the aforementioned decay process results into two DB solitons possessing the same mass, while the momentum of the initial DB "particle" is transferred to the resulting DB fragments so as for the total momentum to be conserved. To this end, we also calculate the average mean position i.e. x B aver = 1 2 ( x B 1 + x B 2 ). Fig. 4(a) illustrates the position expectation value of the first or the second species function and their corresponding fits, as well as the average position obtained from these two. The linear fits to x B 1 , x B 2 and x B aver are shown by dashed lines. For t < t 0 (no splitting of the DB soliton) the linear extrapolation of x B 1 , x B 2 is also shown, where the time instant t 0 is identified by the condition aver . The fitting is chosen to reflect the time for which different entities exist. The ob- aver − x B,f it aver is larger. The latter can be attributed to the coupling of the first two species functions with the remaining of the species functions. Among the two, the somewhat more significant contribution to suggesting that the second species functions possesses a stronger coupling with the remaining species functions.
Following the above-described procedure for the range of initial velocities u 1 /c ∈ {0, 0.8} we show in Fig. 4 we also present the imprinted velocities after the initialization of the dynamics, denoted by v init [75] (i.e. the slope of x B,f it aver ) and the corresponding average velocity v aver = v f ast +v slow 2 . It is observed that for u 1 /c < 0.6 the particle picture makes sense, as v aver = v f ast +v slow 2 = v init , indicating the conservation of momentum. In this case, each component of the DBs splits into two fragments possessing the same mass [see also the relevant occupations in Fig. 2 (a)]. However, for u 1 /c > 0.6 larger deviations from the particle picture are evident. This can be attributed to the fact that one of the emergent fragments within the DB entity is "heavier" from the other one due to a rapid mass redistribution caused by the large velocities.
Next, let us examine how the decay time t D observed during the dynamics depends on both the initial velocity of the DB structure and the particle number of the dark component (N D ). To quantitatively consider the decay time we impose the criterion x B 1 − x B 2 > 2.5w, where w = 1/d refers to the width of the DBs [see also the inset of Fig. 4 (c)]. As it can be seen, in Fig. 4(c), the decay time shows a gradual growth for increasing initial velocity, indicating that in the limit of large velocities we tend to the mean-field picture. On the other hand, Fig. 4(d) presents the corresponding decay time t D for different particle number imbalances N D /N B (here de- Let us also remark here that our results of Figs. 4 (b), (c) in conjuction with the decay mechanism presented before [see also Figs. 3 (i)-(l)] suggest that the resulting solitonic structures do not undergo further splittings at least within the presented evolution time scales. On the one hand, the slow velocity fragment is not a meanfield DB and therefore the observed splitting process of a purely mean-field DB is not expected to occur for such a structure. The slow velocity fragment evolves from darkantidark and bright states to a DWB. The latter possesses no mean-field analogue for a binary system and its stability is presently not known; this constitutes a particularly interesting topic for future study. However, it is expected to be robust as it is a correlated structure produced by quantum fluctuations, and therefore it is not likely that the same mechanism can lead to its decay. In that light, we interpret our results as an indication of stability of this beyond mean-field structure. On the other hand, the fast soliton fragment is more proximal to a mean-field solitary wave and as such it is expected to be more prone to decay. Fig. 4 (b) suggests a velocity u 1 /c > 0.7 which in combination to Fig. 4 (c) shows a decay time t D > 8, being much larger than the decay time of the parent mean-field soliton (t D ≈ 5). There- fore even if a splitting process is permitted it can not be observed within the presented evolution times. Note here that the fast solitary fragment is fully formed only at times t > 11, see Figs. 3 (g), (h). However, we have to be cautious since this fragment is also not of a pure mean-field nature, and therefore its stability is highly non-trivial.

IV. TWO DB SOLITON DYNAMICS
We now turn to the examination of two-soliton dynamics involving collisions. In particular, both fast collisions where the solitons rapidly go through each other (dominated by kinetic energy), as well as slow collisions where the inter-particle interaction (and the potential energy landscape) dictate the dynamics are investigated. Fur-thermore, in all cases to be presented below, the two DB solitons have common inverse widths, and are initialized at x 0 1 = −x 0 2 = −2.5 travelling towards one another with velocity, u 1 /c = −u 2 /c.

A. Fast Soliton Collisions
The results for the case of the fast solitons are summarized in Figs. 5-7. It can be seen that the mean-field picture differs substantially from the beyond-mean-field one. The reason for this is that while the mean-field scenario displays simply the collision of the two DBs [see . It becomes evident that the collision occurring at earlier times entails a rapid population redistribution stage. One can clearly infer that the different modes of entanglement acquire comparable occupations which subsequently become roughly constant throughout the propagation. Notice also the negligible population of the species functions λ 8 (t) and higher, shown in the inset of Fig. 6(a). To further demonstrate the deviations from the mean-field approximation, the evolution of the natural populations is employed [see Fig. 6(b) and (c)]. In particular a rapid depletion is observed up to the collision event for both the dark and the bright soliton components. The fragmentation process is indicated by the decrease of population in the first orbital and the simultaneous emergence of population in the second natural orbital (and in all four orbitals) for the majority species D (minority species B). Closely following the space-time evolution of each of the natural orbitals of the majority species D [see Fig. 7 (a) and (c)] the existence of a density depletion in the first orbital, with an associated density accumulation in the second orbital is observed. Similarly the first two orbitals associated with the minority species B feature bright solitons in both locations, although arguably the first one is (in comparison with the second) slightly more populated at the outer bright, while the second one is definitely more populated (in comparison with the first) at the inner one [see Fig. 7 (b) and (d)]. Moreover, since the picture in the context of the orbitals is somewhat less transparent, in panels (e)-(h) of Fig. 7, profile snapshots of the density evolution of each of the aforementioned natural orbitals is depicted. At initial times [see t 1 = 1 in Fig. 7 (e)], and before the collision takes place, the second orbital of the majority species develops a bright structure [see dashed orange line in panel (e)] within the region at the center [see solid ciel line in panel (e)]. The latter, seemingly, develops into two anti-dark solitons [see t 2 = 2 in Fig. 7 (f ) top] with the respective bright counterpart [see t 2 = 2 in Fig. 7 (f ) bottom] featuring a particle concentration aligned with the merger of the two initial bright solitons, at the collision point. Just after the collision multiple DWB structures are realized for the inner DB solitary wave [see orange and gray dashed and solid lines respectively t 3 = 5 in Fig. 7 (g)]. These persist for larger propagation times [see inner solitonic structure denoted by black arrows at t 4 = 14 in Fig. 7 (h)]. Furthermore, the previously suppressed contribution of the third natural orbital of the minority species B becomes dominant at this larger propagation times, featuring the outer DB soliton pattern [see dashed red lines in the bottom panel of Fig. 7 (h) also indicated by black arrows]. This latter observation leads to further support of the multi-orbital nature of the solitonic structures. We also remark that the same qualitative results remain valid for larger number of atoms N D (results not shown here for brevity).

B. Slow Soliton Collisions
We now turn to an example of a slow collision, summarized in Figs. 8-10, following a similar format as before. Here, at first glance and for times up to t ≈ 3 the dynamics of the one-body densities appears to be similar between the mean-field (shown in panels (a) and (b) of Fig. 8) and the beyond-mean-field scenario (shown in panels (c) and (d) of Fig. 8). However, at later times and as is clearly shown in panels (e)-(j) of Fig. 8, the meanfield approximation is far from an adequate description of the dynamics. At the mean-field level, the in-phase dynamics of the bright components leads to an effective net repulsion [16,21] between them that occurs at large propagation times (t 20) as indicated by the insets in panels (a) and (b) of Fig. 8. This repulsion along with the restoring force of the trap is essentially responsible for the observed phenomenology of the dynamics. On the other hand (and for t > 3), within the beyond mean-field approximation, repulsion is evident at the one-body density level depicted in panels (c) and (d). In particular, it can be clearly discerned that the first mode (k = 1 in Eq. (5)) features a DB solitonic structure which repels and moves outward [see panels (e) and (f )], while the second mode [see panels (g) and (h)] acquires a DB with demonstrable attraction. This inner DB solitonic structure after a first collision at about t ≈ 3.5 seems to settle at a nearly constant distance from each other in an effective "bound state". It is relevant to mention that in the single-component mean-field realm, such a bound state is only relevant in the case of out-of-phase bright structures in the minority component [16,21]. To further illustrate this effective bound state, but also to reveal its multi-orbital nature, in panels (i) and (j) of Fig. 8 a summation over the higher species function contributions, i.e. ing an occupation higher than 0.1%) during the evolution, see the corresponding coefficients λ 1 (t), and λ 2 (t). This is in contrast to the fast collision scenario presented above, where three modes were significantly occupied. The higher-lying species functions possess negligible occupations [see the inset in Fig. 9(a)]. Furthermore, the fragmentation in this velocity regime is less pronounced for the majority species D but of about 50% for the minority species B, shown in Fig. 9 (b), (c) respectively. The respective evolution of the two natural orbitals for both components is shown in Fig. 10 (a)-(d). Notice, that once more the picture in terms of orbitals is less transparent. However, the first orbital can be associated with the outer DB solitary wave. This observation is clearly evident in the bright counterpart [see Fig. 10 (b)] and less pronounced in the corresponding dark component [see Fig. 10 (a)]. The excitation of the second orbital predominantly leads to the formation of the inner soliton. In the latter, traces of attractivity can be identified over time both in the gradually approaching dark component, as well as in the periodically merging (and separating again) bright component. To further elaborate on the resulting dynamics, different instants during propagation are shown in Fig. 10 (e)-(h). It is observed that initially the second natural orbital of the majority species D features a density hump located at the center of the trap. The respective four orbitals of the minority species B are all aligned with the two dark solitons of the first natural orbital of the majority species D [see t 1 = 2 top and bottom panels of Fig. 10 (e) and the relevant inset]. Notice also that the second orbital of the minority species develops two bright solitons. In particular, since the first orbital consists of in-phase bright counterparts the second orbital (due to orthogonality) generates out-of-phase brights. It is this latter effect that naturally leads to the formation of the effective "bound state" discussed above [see also Fig. 8(i)-(j)]. We note here, that an analogous result was also reported in Ref. [33]. The respective third natural orbital develops at the same time density dips. Thus a multi-orbital DB soliton emerges. At the collision point, shown in Fig. 10 (f ), the outer DB soliton is clearly visible (indicated by black arrows), and is accompanied by the inner DB structure formed by the first and second natural orbitals of the majority species D and predominantly by the second orbital of the respective bright soliton component. After the collision event [see t 3 = 4.5 Fig. 10 (g)], the outer DBs can be seen to move outwards. However, a more complicated picture is drawn by the higher natural orbitals. Namely, a dark soliton formed in the second orbital of species D creates an effective potential that traps atoms of the third orbital of species B around the center of the trap, as shown in the inset of Fig. 10 (g). The inner DB structure, indicated by black arrows in the same panel, is aligned with the two dark solitons formed in the third orbital of species B, verifying its multi-orbital nature. Finally, at later times, depicted in Fig. 10 (h), both the outer and inner DBs are moving towards the boundaries (all four are indicated by black arrows, although the dark solitons of the first orbital in the majority component are less clearly discernible). Furthermore, a DW structure developed by the two natural orbitals of the majority species traps the bright solitons formed in the fourth orbital of species B around the center of the trap. Thus, even in this case a DWB structure arises.

V. CONCLUSIONS AND FUTURE CHALLENGES
In the present work, we examined the dynamics of one, as well as of two DB solitons comparing the single-orbital mean-field approximation with the beyond-mean-field approximation featuring multiple inter-species modes of entanglement (and multiple orbitals). A predominant conclusion is that both the oscillation and the interaction between solitons appear to lead to their fragmentation and the population of additional modes (and orbitals). In the typical case scenario, the fragmentation leads to the decay of the initial DB solitons and the simultaneous emergence of an outer DB pertaining to the dominant mode, and an inner one associated chiefly with the first subdominant mode. In the process, we could also identify more complex patterns including, e.g., DWs even within the orbitals of the same (majority D) species, as well as patterns reminiscent of dark-antidark structures.
Furthermore, for the single DB soliton dynamics, it has been shown that the considerations involving momentum redistribution at the particle level between the modes of entanglement capture quite accurately the beyond meanfield decay. Additionally the mean-field DBs were found to be more robust against decaying upon increasing the initial velocity. It was also found that the lifetime decreases with the particle number imbalance between the dark and the bright soliton components.
In the case of collisions, as an effect of the interaction between the constituent DBs, a more involved excitation dynamics when compared to the single DB case is observed. In particular it is found that both fast and slow soliton collisions appear rather disparate from their mean-field analogues as far as the one-body density is concerned. In the former case, this is due to the inner DBs featuring a multi-orbital state more pronounced at the one-body density level. Yet, arguably, slow collisions are fairly dramatic too in repartitioning the density to inner, realizing nearly a bound state, and outer solitary waves. Particularly this bound state is created due to the spontaneous emergence of out-of-phase bright solitons building up in the next-to-dominant order species function for the bright component.
Let us also comment on the corresponding experimental realization of our setup. DB soliton structures have already been realized in [18] by utilizing the hyperfine structure of 87 Rb. Moreover, ensembles consisting of a small number of atoms (∼ 500) have been realized in [76]. Focusing on 87 Rb atoms our dimensionless parameters can be expressed in dimensional form by assuming a transversal trapping frequency ω ⊥ = 2π × 200 Hz. Then, all time scales should be rescaled by 4.09 s and all length scales by 54.7 µm. This yields an axial trapping frequency Ω ≈ 2πHz giving an aspect ratio = Ω ω ⊥ = 5 × 10 −3 . These parameters lie within the range of applicability of the 1D GP theory according to the criterion [ 1, where α ⊥ , α z denote the oscillator length in the transversal and the axial direction respectively, while α is the 3D s-wave scattering length. As a conclusion due to the time scales and Fig. 4 showing the suppression of the decay time with increasing number of atoms, the aforementioned effects might be observed via time of flight imaging in contemporary experiments that report a BEC lifetime of the order of tens of seconds [18].
These findings suggest a number of interesting possibilities for future studies. On the one hand, even at the level of a single species, identifying states involving multiple orbitals such as some of the case examples presented here (e.g., the DWs) would be of particular interest. On the other hand, extending considerations to other contexts, involving, e.g., a higher number of species, or a higherdimensional setting would hold considerable promise. In particular, in the former spinor setting, there have been argued to exist dark-dark-bright and dark-bright-bright solitary waves [77], in which it would be useful to explore not only the spin-independent collisional features (as is the case here), but also the spin-dependent ones. On the other hand, in the higher dimensional setting, the role of dark solitons is played by vortices [78,79]. Such "vortex-bright" solitons are on the one hand rather robust at the mean-field level [78], while also being topologically protected against splitting (of the type observed here) [79]. Hence, it would be interesting and relevant to compare/contrast the dynamical scenarios in the latter case in connection with the former one, and also in light of recent investigations of vorticity features in repulsive condensates beyond the mean-field limit [80,81]. Some of these directions are currently in progress and will be reported in future studies.

Appendix A: Initial state preparation
To initialize the DB soliton dynamics we embed the mean-field solution to the ML-MCTDHB ansatz [see Eqs. (5), (6)]. The number of particles for each species N D , N B , the initial positions x 0 j , the velocities u j and the relative amplitudes ηj η1 for the DB solitons are kept fixed to obtain the corresponding mean-field state. We initialize the algorithm assuming ansatz values for the background chemical potential µ (0) and inverse widths d (0) j of the solitons. The structure of the algorithm is as follows: (a) we obtain the mean-field solution for the background density of the GP equation for µ (0) and µ (0) + δµ using the Newton-Krylov algorithm (b) we calculate N D (µ) = dx|φ D µ (x)| 2 , approximate dN D dµ and update the chemical potential according to the Newton-Raphson method (c) we iterate (a)-(b) until the number of particles converges to N D , thus obtaining a new value for the chemical potential µ (1) (d) we optimize the amplitude of the first soliton η 1 utilizing the Newton-Raphson algorithm (the remaining amplitudes η j =1 and the inverse widths d j are fixed due to the specification of the relative amplitudes is reached. The solutions obtained by the above process are properly normalized and embedded as the first SPF for each species D and B. The remaining SPFs are initialized in a randomly-generated parity-even (for the case of collisions) state and are orthonormalized according to the Gram-Schmidt algorithm. The first species function for each species is initialized in the state where all the particles reside in the corresponding first single particle function. The remaining species functions allow for excitations in the randomly generated orbitals. Finally, the total wavefunction is initialized to the state where only the first species functions for each species is occupied, i.e. λ 1 = 1, and λ i =1 = 0.

Appendix B: Remarks on Convergence
In the present Appendix we briefly discuss the main features of our computational method (ML-MCTDHB) and then demonstrate the convergence of our results.
Within ML-MCTDHB the total many-body wave function is expanded with respect to a time-dependent variationally optimized many-body basis. The latter allows us to obtain converged results with a reduced number of basis states compared to expansions relying on a timeindependent basis. Also, the symmetry of the bosonic species is explicitly employed. Finally, the multi-layer ansatz for the total wave function is based on a coarsegraining cascade, where strongly correlated degrees of freedom are grouped together and treated as subsystems, which mutually couple to each other. The latter enables us to adapt the ansatz to intra-and inter-species correlations and makes ML-MCTDHB a highly flexible tool for simulating the dynamics of e.g. bipartite systems.
The primitive underlying basis for both the Newton-Krylov algorithm and the ML-MCTDHB simulations corresponds to a sine discrete variable representation (DVR) with 1200 grid points. The number of species functions is always kept equal to M D = M B = M . To obtain the corresponding order of approximation we use a specific number of species functions M and single particle functions m D = m B (given within the main text). To track the numerical error and guarantee the accurate performance of the numerical integration for the ML-MCTDHB equations of motion we impose the following overlap criteria | Ψ|Ψ − 1| < 10 −10 and | ϕ i |ϕ j − δ ij | < 10 −10 for the total wavefunction and the SPFs respectively. For the same reasons, also, for the two DB soliton dynamics we ensured that | Ψ M B |P|Ψ M B − 1| < 10 −6 , whereP denotes the many-body parity operator, i.e.
Next, let us elaborate in more detail on the conver-gence of our simulations. To demonstrate the level of our many-body wavefunction truncation scheme we first show the behaviour of the center-of-mass variance calculated both analytically (see below) and numerically. Then, we present the relative error at the core of the dark component (i.e. minimizing the effect of the background) between the 10-(2,4) and 10-(3,4) approximations i.e. by adding one more SPF in the dark component. We remark that the harmonic oscillator potential allows for the separation of the center of mass (CM), R = 1 N i x i , from the relative coordinates r i = x i+1 − x i . This separation enables to reduce the N -body (N = N D + N B ) interacting problem to an interacting N − 1-body problem in the relative coordinates, and a non-interacting one for the CM coordinate. Then the many-body Hamiltonian readsĤ where is the single particle harmonic oscillator Hamiltonian of the center of mass, and H ri is the interaction Hamiltonian in the relative coordinates. However, our calculations within ML-MCTDHB have been performed in the lab frame and as a consequence do not utilize the aforementioned separation of variables. Note that the ML-MCTDHB ansatz (as well as the corresponding mean-field ansatz) does not trivially respect the separation between the CM and relative frame. However, as we shall show below our results can capture the decoupling of the CM motion for the entire N D + N B bosonic cloud. To judge about relative deviations of the ML-MCTDHB propagation with the full Schrödinger equation (and consequently about convergence) we compare the ML-MCTDHB obtained evolution of the center of mass coordinate to the analytical one. Since we are interested only in the one and two-body densities we are able to calculate the first and second statistical moments of the center of mass. The second moment of the CM position (position variance) reads where x i denotes the mean position of the i-th species and x i x j refers to the two-body spatial moment of the position of particles of species i and j respectively [82,83]. By using the Ehrenfest theorem on the center of mass Hamiltonian H R we obtain the exact evolution of the CM position variance R, P denote the spatial coordinate and the momentum operators acting on the CM degree of freedom. This latter expression offers the opportunity to directly measure the deviation between the correlated approach, the meanfield ansatz, and the analytical result [84]. To expose this deviation we numerically calculate the position variance (using different number of orbitals or species functions) and compare with Eq. (B3). This comparison is shown, for the case of a collision between two DBs, in Fig. 11 (a), (b) for fast (u/c = 0.5) and slow (u/c = 0.01) velocities respectively. As it can be seen all the correlated obtained results follow the behaviour of the analytical result and therefore can be considered trustworthy. The observed deviation at long propagation times is of the order of 0.1% lying within the numerical error of the lattice discretization 0.6%. Note that increasing the number of SPFs the correlated approach (see e.g. the 15-(2,4) approximation) tends to a better agreement with the analytical result. As expected, the mean-field result (when compared to the correlated approximation) possesses a larger deviation with the analytics. Also the case of slow velocities shows a comparatively larger deviation (than the case of fast velocities) as it possesses a larger amount of depletion. Next, we illustrate the convergence of our numerical results for another two-body observable, namely the position variance for the CM of the bright soliton Note here that due to the coupling of the CM of the B component with the entirety of the D component no analytic form can be extracted for σ 2 X B (t). Fig. 11 (c) presents σ 2 X B (t) for an increasing number of species functions and/or number of SPFs of both the dark and bright component. As it can be seen the CM for the bright soliton expands during the evolution both within the mean-field (see also the inset) and within the correlated approximation. However, in the correlated case the expansion possesses a larger amplitude when compared to the mean-field approximation. On the convergence side, as shown for increasing number of entanglement modes M the corresponding σ 2 X B (t) graphs are almost indistinguishable (see brown dashed and light blue lines respectively). It is important to note here that the species mean-field case 1-(2,4) differs significantly from the case 15-(2,4). The latter signals the necessity for the inclusion of entanglement modes M in order to fairly capture the DB soliton dynamics (see also below). At the same time, the near-coincidence of the different graphs with the number of modes used is a strong indication of the convergence and the trustworthiness of our results.
Note that within the current calculations we do not ensure the convergence of the ML-MCTDHB simulations at the level of the (N D + N B )-body wavefunction, as such investigations are computationally prohibitive for the large particle numbers considered here. However, we showcase the robustness of the emerging structures in the beyond mean-field dynamics for lower order properties e.g. the one-body density matrix and related one-body quantities. To demonstrate the order of approximation of our results we present the evolution of the relative error integrated within the "core" of the dark soliton. In particular, we calculate the relative one-body density difference of the dark component , where the domain R is such that x ∈ [−10, 10] and refers to the "core" of the dark component; L = 20 is the length of the integrated domain R. More precisely, here, we compare the approximations 10-(2,4) and 10-(3,4), i.e. we include one more SPF in the dark components, thus enlarging the available Hilbert space (H DB ) for the simulation by taking into account 462370 coefficients (in contrast to the 10800 used in the 10 − (2, 4) case). Fig. 11 (c) shows ∆(t) for the oscillation of a DB solitary wave and the cases of both slow and fast collisions between two DBs respectively. As can be seen, for the oscillation of a DB solitary wave ∆(t) before the decay (i.e. t < t D ) lies below 2%, while for t > t D it increases in a linear manner up to 6%. On the other hand, for the case of a slow collision between two DBs ∆(t) increases during the evolution reaching a maximum value of the order of 8.5% during the decay process, while after the decay ∆(t) decreases to 5%. For the fast soliton collision the respective error reaches a maximum value of about 5%. Finally, in Figs. 11 (d), (e), (f ) we present the one-body density evolution of the dark component for the case of a fast oscillation of a DB solitary wave for different approximations, namely 1-(3,4), 10-(3,4) and 10-(2,4) respectively. It is observed that in the species mean-field case, i.e. 1-(3,4) approximation, the two DB structures after the decay produce a jet-like structure and all DB entities are lost. The latter indicates the inescapable necessity of taking more than one species functions for the described excitation dynamics, due to the emergence of the inter-species modes of entanglement. The convergence of our results for a given number of SPFs (m D , m B ) and increasing number of M is observed, verifying that the nature of excitations does not change for increasing number of available inter-species modes of entanglement M . On the other hand, by inspecting the dynamics within the approximations 10-(2,4) and 10-(3,4) we observe a similar qualitative behaviour (confirming also the above described relative error [see Fig. 11 (c)]) before and after the decay process. Similar investigations have also been performed with respect to the bright component SPFs (m B ) not included here for brevity reasons.
Finally, we remark that most of the presented results have been performed within the 15-(2,4) approximation (i.e. using 12780 coefficients), which according to the aforementioned criteria is an adequate approximation for the description of the induced dynamics at the one-body density level. Summarizing, we have observed the same emerging structures during the dynamics for all of the presented orbital/species functions combinations except for the species mean field case (namely the 1-(3,4) approximation).