Correlated Quantum Dynamics of a Single Atom Collisionally Coupled to an Ultracold Finite Bosonic Ensemble

We explore the correlated quantum dynamics of a single atom, regarded as an open system, with a spatio-temporally localized coupling to a finite bosonic environment. The single atom, initially prepared in a coherent state of low energy, oscillates in a one-dimensional harmonic trap and thereby periodically penetrates an interacting ensemble of $N_A$ bosons, held in a displaced trap. We show that the inter-species energy transfer accelerates with increasing $N_A$ and becomes less complete at the same time. System-environment correlations prove to be significant except for times when the excess energy distribution among the subsystems is highly imbalanced. These correlations result in incoherent energy transfer processes, which accelerate the early energy donation of the single atom and stochastically favour certain energy transfer channels depending on the instantaneous direction of transfer. Concerning the subsystem states, the energy transfer is mediated by non-coherent states of the single atom and manifests itself in singlet and doublet excitations in the finite bosonic environment. These comprehensive insights into the non-equilibrium quantum dynamics of an open system are gained by ab-initio simulations of the total system with the recently developed Multi-Layer Multi-Configuration Time-Dependent Hartree Method for Bosons.


Introduction
Many physically relevant quantum systems are in fact open. Intriguing effects in e.g. condensed matter physics [1], quantum optics [2], molecules or light harvesting complexes [3,4] are intimately related to environmentally induced dissipation and decoherence and thus require a careful treatment beyond the unitary dynamics of the time-dependent Schrödinger equations. As one can often neither experimentally control nor theoretically describe all the environmental degrees of freedom, a variety of theoretical methods has been developed for an effective description of the reduced dynamics of the open quantum system of interest over the last decades [5].
Due to their high degree of controllability and, in particular, isolatedness [6], ensembles of ultracold atoms or ions serve as ideal systems in order to systematically study the dynamics of open quantum systems allowing for various perspectives on this subject [7]. Open quantum system dynamics has been implemented by digital quantum simulators [8,9] or by partitioning an ultracold atomic ensemble into carefully coupled, distinguishable subsystems [10]. As a matter of fact, many impurity problems can also be viewed as open quantum system settings: enormous experimental progress such as [11][12][13][14] allows to prepare a few impurities or even a single one in an ensemble of atoms in order to study thermalization and atom loss mechanisms [15], transport and polaron physics [16][17][18] or the damping of the breathing mode [19,20]. Such implementations of open quantum systems offer the unique flexibility to tune both the character and strength of the system-environment coupling [19] and the environmental properties [18]. Moreover, dissipation and environment engineering can also be employed for controlling many-body dynamics [21], state preparation [22] as well as quantum computation [23,24].
In this work, we theoretically study the open quantum system dynamics of a single atom with a weak spatio-temporally localized coupling to a finite bosonic environment. By considering a binary mixture of neutral atoms interacting via contact interaction both within and between the species, the local character of the coupling is realized. In order to localize the coupling also in time, species-selective one-dimensional trapping potentials as well as a particular initial condition are considered: the single atom is initially displaced from the centre of its harmonic trap, i.e. resides in a coherent state, while the bosonic ensemble is prepared in the ground state of a harmonic trapping potential being shifted from the trap of the single atom. Thereby, only during a certain phase of the single atom oscillation period both subsystems couple. In a different context, this kind of coupling has already been realized effectively experimentally [25]. The basic ingredient of such a coupling, a single binary collision, can result in entanglement between the collision partners [26][27][28][29], which has been shown to depend on the details such as the scattering phase shift, the mass ratio as well as the relative momentum of the atoms [30]. As a consequence, significant correlations between the single atom and the finite bosonic environment can occur after many collisions despite of possibly weak interactions (cf. also [27,31,32] in this context), undermining a mean-field approximation for the two species on a longer time scale. We note that the scenario under consideration might seem to be reminiscent of the quantum Newton's cradle [33,34]. Yet instead of investigating the (absence of) thermalization in a closed system of indistinguishable constituents we are concerned with unravelling the interplay of energy transfer between distinguishable subsystems and the emergence of correlations when systematically increasing the number of environmental degrees of freedom N A .
Rather than simulating the reduced dynamics of the single atom only, we employ the recently developed ab initio Multi-Layer Multi-Configuration Time-Dependent Hartree Method for Bosons (ML-MCTDHB) [35,36] for obtaining the numerically exact non-equilibrium quantum dynamics of the whole system for various numbers of environmental degrees of freedom. Such a closed system perspective on an open quantum system problem gives the unique opportunity to investigate not only the dynamics of the open system but also its impact on the environment and, moreover, to systematically uncover correlations between the two subsystems (cf. also e.g. [37,38]). Full numerical simulations are supplemented with analytical and perturbative considerations.
Having introduced our setup and discussed the possibility of controlling the energy transfer channels by tuning the separation of the species-selective traps in sect. 2, the energy transfer between the single atom and the finite bosonic environment is treated in sect. 3. Here, we show that the energy transfer between the subsystems is accelerated with increasing N A due to a level splitting of the involved excited many-body states. At the same time, the energy transfer becomes less complete so that after an initial slip the energy distribution among the subsystems features fluctuations about quite a balanced one for larger N A . In sect. 4, we investigate how the subsystem states are affected by the system-environment coupling: By inspecting the Husimi phase space distribution of the single atom, we conclude that the energy transfer is mediated by non-coherent states manifesting themselves in a drastic deformation of the phase space distribution during relatively short time slots. Concerning the environment, depletion oscillations are observed, whose amplitude decreases with increasing N A when fixing the initial displacement of the single atom. In sect. 5, we unravel the oscillatory emergence and decay of inter-species correlations showing that whenever the excess energy distribution among the subsystems is highly imbalanced correlations have to be strongly suppressed. The maximal attained correlations turn out to be independent of the considered number of environmental degrees of freedom N A . By analysing the subsystem energy distribution among the so-called natural orbitals, we show how inter-species correlations result in incoherent energy transfer processes, which accelerate the early energy donation of the single atom. Most importantly, we ultimately uncover the interplay between subsystem excitations and correlations by means of a Fock space excitation analysis in sect. 6, characterizing the environmental excitations and unravelling correlations between subsystem excitations with a tailored correlation measure. Thereby, we can show that inter-species correlations (dis-)favour certain energy transfer channels depending on the instantaneous direction of transfer in general. The conclusions and outlook in sect. 7 contain also a proposal for an experimental realization of this open quantum system.

Setup and initial conditions
The subject of this work is a bipartite system confined to a single spatial dimension and consisting of two atomic species: a single atom is collisionally coupled to a finite bosonic environment (cf. fig. 1 (a)). We model the intra-and inter-species shortrange interactions with a contact potential as it is usually done in the ultracold s-wave scattering limit [39]. In addition, we assume two external harmonic potentials for a species-selective trapping of the atoms. In the following, we adopt a shorthand notation labelling the bosonic environment with A and the single atom with B. Throughout this work, all quantities are given in terms of natural harmonic oscillator units of B, which base on the mass of the single atom m B and the frequency of its trap ω B . Lengths and energies are thus given in terms of /(m B ω B ) and ω B , respectively, times are measured in units of ω −1 B . The Hamiltonian of the composite system takes the form with the Hamiltonian of the environmentĤ A , the open systemĤ B and the intersubsystem couplingĤ AB : The dimensionless representation of the Hamiltonian involves the mass ratio β = m A /m B and the frequency ratio α = ω A /ω B . R > 0 denotes the separation between the traps. We refer to the intra-and inter-species interaction strengths as g A and g AB , respectively. We emphasize that the variety of parameters specifying the Hamiltonian (1) accounts for the versatility of the considered system, which to explore in its complexity goes far beyond this work. Rather than this, we fix g A , g AB to experimentally feasible values by considering two 87 Rb hyperfine states of comparable intra-and interspecies scattering lengths, implying β = 1. Assuming ω B = 10 Hz and a transverse to longitudinal trapping frequency ratio of about 68 for both species, we are left with g A = g AB = 0.08 for the one-dimensional coupling constants after dimensional reduction [40]. We initialize B in a coherent state |z : where |u σ n denotes the n th harmonic oscillator (HO) eigenfunction of the species 2 > 0 such that |z equals theĤ B ground state displaced by a distance d to the right from the B species trap centre. The environment, A, in turn is prepared in the ground state ofĤ A . In order to realize the desired spatio-temporally localized coupling, whose impact on the inter-species energy transfer and emergence of correlations shall be investigated, we require R and d to be such that there is (i) no inter-species overlap at t = 0 and (ii) finite overlap after half of a B atom oscillation period, i.e. at t ≈ π. For more details on the initial state preparation and subsequent numerically exact propagation with our ab initio ML-MCTDHB method, see Appendix A and Appendix B.
In order to characterize the spatio-temporally localized coupling, we expressĤ AB in terms of the HO eigenbasis {|u σ n }, which serves as a natural choice for the low-energy and weak-coupling regime we are interested in, withâ ( †) i denoting the bosonic annihilation (creation) operator corresponding to |u A i and the interaction-matrix elements given by which obviously depends on α and R. Eventually, we are interested in rather small displacements d for which the matrix element v 1001 is of major importance for the dynamics. For e.g. α = 1, we find that v 1001 = (1 − R 2 ) exp(−R 2 /2)/(2 √ 2π). Hence, the scattering channel "|u A 0 |u B 1 → |u A 1 |u B 0 " is completely suppressed at R = 1. As it turns out, the choice R = 1.2 is well suitable for our purposes since the most important energetically allowed scattering channels shall not be suppressed.
Although we have fixed some of the parameters, the system still features a high sensitivity to α, d, N A and thus controllability. For N A = 1, analytical expressions for the energy spectrum and the wave functions can be worked out and trap-induced shape resonances between molecular and trap states are detectable [41]. Since this two-body problem has already been treated in detail, we shall rather focus on the impact of the environment size, considering N A = 2, ..., 10.

Energy transfer
In order to understand which kind of energy transfer processes between the species are feasible with the spatio-temporally localized coupling and how efficient these are, we examine the energies of the subsystems, identified with Ĥ A t and Ĥ B t , separately. Since the interaction energy Ĥ AB t cannot be attributed to the energy content of a single species, the aforementioned identification is problematic, in general. However, the spatio-temporally localized coupling always allows to find times during a B atom oscillation at which the inter-species interaction is essentially negligible so that Ĥ A t and Ĥ B t measure how the energy is distributed among the open system and its environment. Because of the weak and spatio-temporally localized inter-species coupling, the short-time dynamics taking place on the time-scale of the free oscillation period of the B atom, i.e. T = 2π, is clearly separated from the long-time dynamics of the energy transfer between the two species. When quantifying times, we will use the term "B oscillations" always with reference to free harmonic oscillations of the B atom. Due to this time-scale separation, we present the expectation value of any physical observableÔ as locally time-averaged over one free B oscillation period, Thereby, we average out fluctuations on the short time-scale of a B oscillation, giving a clearer view on the long-time dynamics. The fast short-time dynamics is captured by the variance of the local time-average, In some cases, it turns out to be insightful to encode this information in the figures by accompanying the lines corresponding toŌ(t) with tubes indicating the standard deviation. We prepare the system such that the two species have no spatial overlap at t = 0. Therefore, the initial energy of the single atom reads By adjusting the displacement d, we can thus vary the amount of deposited excess energy ε ≡ d 2 /2, which is available for energy transfer processes between the subsystems. These energy transfer processes are captured in terms of the intra-species excess energies ε σ t = Ĥ σ t −E σ gs with E σ gs denoting the ground state energy ofĤ σ , σ = A, B. We employ the normalized excess energy of B, ∆ B t = ε B t /ε, in order to measure how balanced the deposited energy is distributed among the subsystems: At instants for which Ĥ AB t is negligible, one obviously finds ε A t /ε = 1−∆ B t . Thus ∆ B t ≈ 1(0) implies then that almost all excess energy is stored in the B (A) species, corresponding to a maximally imbalanced excess energy distribution, while ∆ B t ≈ 0.5 refers to a balanced distribution. Due to the initial imbalance ∆ B 0 = 1, the B atom will first donate energy to its environment implying an overall decrease of ∆ B t until ∆ B t reaches a minimum. We call the energy donation of B the more efficient the closer to zero this minimum is. Afterwards, the B atom will generically accept energy from its environment, resulting in an overall increase of ∆ B t until it reaches a maximum. We call such an energy transfer cycle the more complete the closer to unity this maximum is.
Whether B can induce excitations in the A species is predominantly determined by the ratio of ε and the excitation gap ofĤ A , which equals α when neglecting the intra-species interaction. If /α 1, the B atom cannot excite A bosons and thus the environment constitutes a weak additional potential for B without experiencing any back-action. Non-trivial open-quantum system dynamics can only be expected for ε α. In the following, we concentrate on the threshold case ε ≈ α so that not too In order to investigate the efficiency of the B energy donation in dependence of α, d, we inspect the maximal fraction of the total excess energy transferred to A, i.e. max t (1 − ∆ B t ), for N A = 4. Fig. 1 (b) reflects that the two species can most efficiently exchange energy when the traps are in resonance, i.e. α = 1. In this case, approximately 68% of ε is transferred to A independently of d. For α = 11/12, 13/12, the energy transfer is much less efficient thus indicating a sharp efficiency resonance around α = 1. For N A = 2, 3, even more excess energy is donated (plots not shown). On these grounds, we fix α = 1. If not stated otherwise, d = 1.5 is assumed realizing the threshold case ε ≈ α, so that the only free parameter is now N A . Fig. 2 shows the normalized excess energy∆ B (t) for N A = 2, 4, 7, 10. As a starting point, we discuss the case N A = 2. The initial excess energy distribution is maximally imbalanced, i.e. ∆ B t=0 = 1. As B donates energy to its environment,∆ B (t) monotonously decreases until reaching a global minimum at t ≈ 289, implying a directed energy transfer from B to A. The energy transfer decelerates at t ≈ 145 and reaccelerates at t ≈ 190. When∆ B (t) attains its minimum at t ≈ 289, the reverse energy transfer process is initiated and B accepts the energy that was donated to A in the first place. The phase of decelerated energy transfer can also be observed between t ≈ 395 and t ≈ 440. At t ≈ 584,∆ B (t) attains a local maximum and the excess energy is almost completely restored in B (H B (t = 584) 0.97 Ĥ B (t = 0)).
Now we turn to the impact of N A on∆ B (t). The oscillatory evolution of∆ B (t) is similar for N A = 2 and N A = 4, 7, 10. For example, we identify the minimum of∆ B (t) for N A = 2 at t ≈ 289 with the minimum of∆ B (t) for N A = 4 at t ≈ 195. Analogously, the maximum of∆ B (t) for N A = 2 at t ≈ 584 indicating that B accepted much of the excess energy ε from the environment corresponds to the maximum of∆ B (t) for N A = 4 at t ≈ 345. Already after a few B oscillations,∆ B (t) is considerably smaller for higher N A than for N A = 2. Thus, all results shown in fig. 2 indicate that the overall energy transfer process accelerates with increasing N A . Classically, this can be understood in terms of an increasing number of collision partners, making excitations in A more likely to occur within a B oscillation. This manifests itself in the effective mean-field coupling strength g AB N A . For a quantum-mechanical explanation, we firstly extract the energy transfer time-scale T cycle by applying compressed sensing [42,43] to the Ĥ B t data, which gives an accurate, sparse frequency spectrum despite of the short signal length. In fig.  2 (b), we depict T cycle = 1/ω 2 with ω 2 denoting the angular frequency of the dominant peak aside from the DC peak. We observe excellent agreement of T cycle with 2π/∆E 1 , where ∆E 1 denotes the energy gap between the first and second stationary excited state ofĤ obtained with ML-MCTDHB. In order to obtain a pictorial understanding of this level splitting ∆E 1 increase with N A , we have repeated our analysis for the case g A = 0, allowing for the analytical expression within first-order degenerate perturbation theory w.r.t.Ĥ AB . This perturbative result is in qualitative agreement with the numerics in fig. 2 (c). Denoting |n 0 , n 1 , ... A HO as a bosonic number state with n i bosons of type A occupying |u A i , it is immediately clear that the inter-species interaction raises the unperturbed state |N A , 0 A HO ⊗ |u B 1 more in energy than the unperturbed state Moreover, the results shown in fig. 2 (a) for N A > 2 feature less pronounced energy minima and maxima. In our terminology, the energy donation becomes less efficient and the whole transfer cycle less complete as N A is increased. For N A ≥ 7, this culminates in the following behaviour: After the initial phase of energy donation from B to A lasting a few tens of B oscillations,∆ B (t) fluctuates around the value for a balanced distribution, i.e. ∼ 0.5...0.6.
We note that the above mentioned phases of decelerated energy transfer are not observed for N A > 4. Instead, additional structure emerges in the form of not only deceleration but a temporary phase in which B accepts energy. This culminates in the appearance of an additional minimum and maximum, e.g. for N A = 7 at t ≈ 70 and t ≈ 100, respectively.

State analysis of the subsystems
In this chapter, we aim at a pictorial understanding of the subsystem dynamics. For this purpose, we study the impact of the spatio-temporarily localized coupling on the short time-evolution of the B atom density and the reduced one-body density of the A species (sect. 4.1). Afterwards, we employ the Husimi phase-space representation of the state of B for investigating to which extent its initial coherence is affected by the environment (sect. 4.2). Likewise, we shall also characterize the environment whilst inspecting its depletion from a perfectly condensed state (sect. 4.3). Further insights into the environmental dynamics are given in sect. 6.2 focusing on intra-species excitations.

One-body densities
Any predictions about measurements upon the single B atom alone can be inferred from its reduced density operatorρ B t , which is thus regarded as the state of the B atom: Here, Tr A denotes a partial trace over all A bosons. Analogously, we defineρ A t as a partial trace over the B atom and all A bosons but one.
Due to inter-species interactions, the single atom leaves a trace in the spatial density profile of A. Fig. 3 shows the reduced densities, ρ σ

Coherence analysis
However, the spatial density yields no information about how coherent the stateρ B t of the single atom remains. Instead, the Husimi distribution serves as a natural quantity to unravel deviations from a coherent state characterized by an isotropic Gaussian distribution. Due to its positive-semi-definiteness, Q B t (z, z * ) allows for an interpretation as the probability density to find the system in the coherent state |z . The real and imaginary parts of z ≡ re iϕ can be identified with position and momentum in a phase space representation, i.e. (z) =x/ √ 2, (z) =p/ √ 2. Since B performs harmonic oscillations in its trap, which entail rotations of Q B t in phase space, the subsequent phase space analysis is performed in the co-rotating frame In fig. 4 (a),∆ B (t) is shown as well as snapshots ofQ B t (z, z * ) at characteristic points in time for N A = 2 and d = 1.5, 2.5. The initial Husimi distributions are Gaussians centred aroundx = d as desired.
As a matter of fact,Q B t (z, z * ) reflects the dissipative dynamics of the B atom. The distancer t of the mean valuez t ofQ B t (z, z * ) from the origin decreases (increases) with decreasing (increasing) energy of the B atom (cf. fig. 4 (b)).
On the other hand, the shape ofQ B t (z, z * ) provides us with information about the quantum state of B. In all our investigations for low excess energies,Q B t (z, z * ) resembles a Gaussian, undergoing a breathing into the ϕ-and r-directions with a relatively constant mean valuez on the time-scale of a B oscillation, during most of the phases of the dynamics. Drastic shape changes are only observed during short phases lasting a few B oscillations and are accompanied with drifts of the phaseφ t ofz t . During these phases,Q B t (z, z * ) features a less symmetric shape, e.g. of a squeezed state as for d = 1.5 at t = 76.0. As can be inferred from fig. 4 (a), this short phase coincides with a rapid energy flux from B to A. This suggests that the directed inter-species energy transfer is mediated through non-coherent B states.
In our frame of reference, the time-dependence ofφ t is due to the collisional coupling to the environment. From the Husimi distribution and thez t trajectory in fig. 4 (b), we infer thatz accumulates a collisional phase shift. At t = 282.9, i.e. when∆ B (t) is minimal, a phase shift ofφ t ≈ π is observed for d = 1.5. Remarkably, at t = 590.2, the phase shiftφ t ≈ 2π indicates that the initial state is almost fully recovered.
In the case of a displacement d = 2.5, the excess energy ε is almost three times larger than for d = 1.5. First of all, we note that the extrema of∆ B (t) are less pronounced than for d = 1.5 such that the recovery of the initial energy is very incomplete. Ultimately, one finding is similar for the higher displacement: most of the time,Q B t (z, z * ) resembles a Gaussian. Again, only within short phases the shape changes drastically.Q B t (z, z * ) then differs significantly from a Gaussian, showing a delocalization of a strongly squeezed state as observed for d = 2.5 at t = 389.0. We remark that for displacements d ≥ 2.5, our results are of lower accuracy and therefore of qualitative character.
In order to quantify the coherence of the B quantum state, we employ an operator norm to measure the distance ofρ B t from its closest coherent state: We note that C t ∈ [0, 1) and C t = 0 iff B is in a coherent state. Focussing firstly on N A = 2 and d = 1.5, the implications drawn from the evolution ofQ B t (z, z * ) persist as depicted in fig. 4: At instants of extremely imbalanced energy distribution among A and B, B is very close to a coherent state as indicated byC(t) ≈ 0 at t ≈ 289 and t ≈ 584. For d = 1.5 and N A = 2, we can relate the phases around local minima ofC(t), e.g. at t ≈ 145, to the phases of decelerated energy transfer. For larger d and N A = 2, the initial coherence is less restored at instants of extremal excess energy imbalance and in between the state of B deviates more strongly from a coherent one. The recovery of the coherence also becomes less complete and the maximum ofC(t) increases as the size of the environment is increased while d is kept fixed (cf. N A = 4, 10). Moreover, as we will see in sect. 5.1, the coherence measure C t strongly resembles the time evolution of the inter-species correlations. Therefore, we conclude that the temporal deviations from a coherent state are caused byρ B t becoming mixed.

State characterization of the bosonic environment
In this section, we investigate how the initially condensed state of the finite bosonic environment evolves structurally, thereby finding oscillations of the A species depletion, whose amplitude is suppressed when increasing N A while keeping ε fixed. We stress that the term "condensed" is not used in the quantum-statistical sense but shall refer to situations when all N A bosons approximately occupy the same single-particle state. For this analysis, we employ the concept of natural orbitals (NOs) and natural populations (NPs) [44], which are defined as the eigenvectors and eigenvalues of the reduced density operator of a certain subsystem in general. The spectral decomposition of the reduced one-body density operators for the species, σ = A, B, in particular reads: where m σ denotes the number of considered time-dependent single-particle basis functions in the ML-MCTDHB method (cf. Appendix A). Due to our normalization of reduced density operators, the λ σ i (t) ∈ [0, 1] add up to unity. In the following, we label the NPs in a decreasing sequence λ σ i (t) ≥ λ σ i+1 (t) if not stated otherwise. For characterizing the state of the bosonic ensemble, we use the NP distribution corresponding to the density operator of a single A bosonρ A t : If λ A 1 (t) ≈ 1, the bosonic ensemble is called condensed [45], whereas slight deviations from this case indicate quantum depletion. Since it is conceptually very difficult to relate the NP distribution λ A i (t) to intra-species (and also inter-species) correlations, we may only regard the λ A i (t) as a measure for how mixed the state of a single A atom is.
For d = 1.5, only two NOs are actually contributing toρ A t -the other NPs are smaller than 8.7 · 10 −3 . Therefore, we depict only the first two NPs for N A = 2, 4, 7, 10 in fig. 6. Due to the weak intra-species interaction strength, the initial depletion of the bosonic ensembles is negligibly small, 1 − λ A 1 (0) < 10 −3 . One clearly observes that the A species becomes dynamically depleted and afterwards approximately condenses again in an oscillatory manner. The instants of minimal depletion coincide with the instants of maximal excess energy imbalance between the subsystems for N A = 2. Therefore, the two bosons behave collectively in the sense that both approximately occupy the same single-particle state not only at phases when the A species is effectively in the ground state ofĤ A but also when 1 −∆ B (t) becomes maximal. For larger N A , the depletion minima turn out to be less strictly synchronized with the∆ B (t) extrema. We will encounter the same finding concerning the minima of inter-species correlations in sect. 5.1 and discuss the details there.
Strikingly, the maximal depletion is significantly decreasing with increasing N A . In Appendix C, we show that when neglecting the intra-species interaction the higher order NPs are bounded by λ A i (t) d 2 /(2N A ), i ≥ 2, for sufficiently large N A , which is a consequence of the gapped excitation spectrum, the extensitivity of the energy of the A species and the fixed excess energy d 2 /2. Although this line of argument neglects intra-species interactions, which will become important for larger N A ‡, the numerically obtained depletions satisfy the above bound well.

Correlation analysis of the subsystems
As we consider a bipartite splitting of our total system, a Schmidt decomposition shows that the eigenvalue distributions ofρ B t and the reduced density operator for the A species, coincide. Here, we study both the emergence of inter-species correlations on a short time-scale and their long-time dynamics in sect. 5.1. In particular, we show analytically that these correlations essentially vanish whenever the excess energy ε σ t of one species σ = A, B is much lower than the excitation gap ofĤ σ . In sect. 5.2, we then unravel the energy transfer between the species by inspecting the energy stored in the respective NOs |Φ A i (t) , |ϕ B i (t) of the subsystems. Thereby, we identify incoherent energy transfer processes, which are related to inter-species correlations and shown to accelerate the early energy donation of the B atom to the A species.

Inter-species correlations
The NPs λ B i (t) are directly connected to inter-subsystem correlations: Since the B species consists of a single atom being distinguishable from the A atoms, and since the bipartite system always stays in a pure state, the initial pure state of the B atom characterized by λ B 1 (0) = 1 can only become mixed if inter-species correlations are present. Deviations from λ B 1 (t) = 1 indicate entanglement between the single atom and the species of A bosons. In order to quantify these inter-subsystem correlations, we employ the von Neumann entanglement entropy: ‡ The situation is rather involved for bosons in a one-dimensional harmonic trap, as the ratio of g A and the (local) linear density, which depends on g A and N A , determines the effective interaction strength. which vanishes iff inter-species correlations are absent, i.e. |Ψ t = |Φ A 1 (t) |ϕ B 1 (t) . For all considered N A , the von Neumann entropy of both the ground state ofĤ and our initial condition |Ψ 0 with d = 1.5 does not exceed 3.5 · 10 −3 , which is negligible compared to the dynamically attained values depicted in fig. 7. Thus, we may conclude that (i) the spatio-temporarily localized coupling requires a certain amount of excess energy for the subsystems to entangle and (ii) inter-species correlations are dynamically established via interferences involving excited states.
In fig. 7 (a), we present the short-time dynamics of the S vN (t) for N A = 2. The entanglement between the single B atom and the two A bosons is built up in a stepwise manner with each collision. Since the local maxima of S vN (t) are delayed w.r.t. the corresponding maxima of Ĥ AB t , we conclude that a finite interaction time is required in order to enhance correlations. After a few B oscillations, the effective interaction time is increased since (i) dissipation moves the left turning point of the B atom from x = −d < −R towards the centre of the A species at x = −R and (ii) the A / B species density oscillations become synchronized (cf. fig. 3). As a consequence, the overall slope of S vN (t) increases during the first few B oscillations. We note that such a step-wise emergence of correlations has been reported also in [27] for two atoms colliding in a single harmonic trap.
Concerning the long-time dynamics, we first focus on the N A = 2 case in fig. 7  (b). The correlation increase continues until a maximum at t ≈ 78, which is followed by a minimum around t ≈ 144 of modest depth, a further maximum at t ≈ 210 and a deep minimum at t ≈ 288. This alternating sequence of maxima and minima is repeated thereafter. In passing, we note that to some extent similar entropy oscillations with a long period comprising many collisions have been reported for two atoms in a single harmonic trap interacting with repulsive and attractive contact interaction in [27] and [31], respectively. As indicated by the vertical lines, the deep minima ofS vN (t) at t ≈ 288 and t ≈ 580 coincide very well with instants of both a maximally imbalanced excess energy distribution among the species (cf. fig. 2) and minimal A species depletion (cf. fig. 6).
This intimate relationship between maximal excess energy imbalance and absence of correlations can be understood analytically by expressing Ĥ B t as a function(al) of the NPs and NOs, and identifying Σ B with the energy content of the i-th NO. In Appendix D, we show that ε B t being much smaller than the excitation gap At those instants, the presence of an excitation gap prevents correlations between the subsystems due to a lack of orthonormal states The same line of arguments analogously holds also for instants when ε A t is much smaller than the excitation gap ofĤ A . In between instants of maximal excess energy imbalance, more states are energetically accessible such that inter-species correlations can be established. Increasing N A while keeping d = 1.5 just above the excitation threshold alters theS vN (t) dynamics in the following way: (i) Inter-species correlations initially emerge faster because binary collisions become more likely within a B oscillation (cf. N A = 2, 4, 10 curves for t < 35).
(ii) The depth of shallower minima § decreases (cf. the N A = 4 curve at e.g. t ≈ 83) such that neighbouring pairs of local maxima (cf. the N A = 4 curve at e.g. t ≈ 60 and t ≈ 107) degenerate to a single maximum for N A > 4 (all entropy minima for e.g. N A = 10 correspond to the deep minima observed for smaller N A ).
(iii) We find the coincidence of the deep entropy minima with extrema in the excess energy imbalance∆ B (t) to be less established. For e.g. N A = 4, the first deepS vN (t) minimum at t ≈ 169 is attained before the first∆ B (t) minimum at t ≈ 195. Since the overall energy donation of the B atom is less efficient for N A = 4 (compared to N A = 2), S vN (t) is not forced to attain a deep minimum at the first minimum of∆ B (t) in the sense of the above analytical line of argument. For N A = 10,∆ B (t) essentially fluctuates around 0.6 for t > 60 such that the excitation gaps do not impose any restrictions on S vN (t) and now∆ B (t) minima coincide withS vN (t) maxima (for 60 < t < 300).
(iv) The minimal value ofS vN (t) for t > 0 increases, which goes hand in hand with the energy donation of the B atom becoming less efficient and the energy transfer cycle becoming less complete.
(v) The maximal values ofS vN (t) are all comparable for N A ≤ 10. We note that the persistence of inter-species correlations when increasing N A does not contradict the § For N A = 1, these minima even attain a depth similar to the deep minima (plot not shown).
We note that the deepS vN (t) minima remain synchronized with the A species depletion minima for all considered N A (for N A = 10 and t > 400 the depletion minima, however, are quite washed out).
fact that the A species becomes simultaneously more condensed (cf. sect. 4.3), which can be easily seen by a minimal example provided in Appendix E.

Energy distribution among natural orbitals and incoherent transfer processes
In this section, we unravel the impact of the previously identified inter-species correlations on the energy transfer and identify incoherent transfer processes. For this purpose, we will evaluate the contribution of Σ B i (t) to Ĥ B t and complement this energy decomposition by analysing how the energy of the A species is distributed among the A species NOs |Φ A i (t) : Here, we have introduced the energy content of the i-th A species NO: . We remark that the presence of intra-species interaction prevents us from expressing Ĥ A t as a function(al) of the NPs and NOs ofρ A t in a similar fashion as in (15). The highly challenging task of diagonalizing the reduced density operator η A t of the whole A species for obtaining the N A -body states |Φ A i (t) , however, can be faced efficiently with the ML-MCTDHB method due to its beneficial representation of the many-body wave function (cf. Appendix A). We note that the NO energy contents Σ A/B i (t) constitute system-immanent, basis-independent quantities as they result from the diagonalization of reduced density operators.
For the considered displacement d = 1.5 and all considered N A , we find that only two NOs contribute toρ B t andη A t (λ B i 6.4·10 −3 for i > 2) so that we only analyse their energy contents in fig. 8. Larger displacements generically result in more populated NOs and thus a much more involved dynamics, whose analysis goes beyond the scope of this work. In fig. 8, we bring the evolution of the second NPλ B 2 (t) face to face with the energy contents of the dominant and the second dominant NOsΣ does not contain information about how much energy of the σ species is actually stored in the corresponding NO, we moreover quantify the contribution of the second NO to the subsystem energy in terms of the ratio λ B 2 Σ σ 2 / Ĥ σ (t). We clearly observe that the energy content of the dominant NO,Σ A/B 1 (t), qualitatively resembles the subsystem energyH A/B (t) evolution (cf. fig. 2). However, at phases when inter-species correlations are present, i.e. whenλ B 2 (t) 3 · 10 −2 , a significant fraction of the subsystem energy is stored in the respective second dominant NO. This fraction can become as large as 23%. Since essentially all subsystem energy is stored in the respective dominant NO whenever∆ B (t) is extremal (cf. at t ≈ 289 and t ≈ 584) and since the evolution of λ B 2 Σ A 2 / Ĥ A (t) appears to be synchronized with the dynamics of λ B 2 Σ B 2 / Ĥ B (t), we conclude that the energy transfer is mediated via incoherent processes in the following sense: Rather than keeping all its instantaneous energy in the dominant NO, the B atom shuffles energy from the dominant to the second NO and back while donating (accepting) energy to (from) the A species. The A species, in turn, accepts (donates) energy from (to) the B atom while shuffling energy from the dominant to the second NO and back. As predicted by the analytical line of argument in Appendix D, one can furthermore witness how the intra-species excitation gap of the σ species makes the second NO to energetically separate drastically from the dominant one when ε σ t becomes too small (cf.Σ B 2 (t) at t ≈ 289 andΣ A 2 (t) at t ≈ 584). Finally, we observe thatΣ (t) holds most of the time. Consequently, if the B atom has been detected in the NO of lower (higher) energy, the A species is found almost certainly in the NO of lower (higher) energy and vice versa, according to the above Schmidt decomposition.
As discussed at the end of sect. 5.1 in detail, the inter-species correlation dynamics becomes less strictly synchronized with the∆ B (t) evolution when increasing N A . Nevertheless, the above overall picture of the NO energy content evolution remains valid for all N A ≤ 10. Concerning N A = 2, we summarize peculiarities, which are diminished or absent for larger environment sizes: Firstly, the decelerated energy transfer around t ≈ 160 coincides with a phase of significant contribution of the second NO toH B (t) so that the incoherent processes can be made partially responsible for the delay. Secondly, pairs of subsequent in between consecutive deep minima (e.g. at t ≈ 91 and t ≈ 214) are related to the observed pairs of localS vN (t) maxima and, thus, merge to a single maximum for larger N A . Finally, we analyse how the overall energy transfer is influenced by these incoherent energy transfer processes. The ML-MCTDHB method allows to manually switch off inter-species correlations in order to clarify their role for the dynamics (cf. Appendix A). By comparison, we have found for all considered N A that the presence of interspecies correlations accelerates the energy donation of the B atom at first, i.e. for t 50, which is plausible since more energy transfer channels are open. Whether the overall energy donation is accelerated, however, depends on N A : For N A = 2, we observe an acceleration by a factor of two, while for N A = 7 no overall acceleration is found (plots not shown).

Excitations in Fock space and their correlations
Operating in the weak interaction regime, we naturally define intra-system excitations as occupations of excited single-particle eigenstates corresponding to the respective HO Hamiltonian. The joint probability for finding (n 0 , n 1 , ...) A bosons in the respective HO eigenstates and the B atom in its j-th HO eigenstate is given by: In this chapter, we firstly show that the actual long-time excitation dynamics takes place only in certain active subspaces being mutually decoupled (sect. 6.1). Secondly, the A species excitations are shown to be governed by singlet and delayed doublet excitations (sect. 6.2). The findings of these two sections are explained by means of a tailored time-dependent perturbation theory. Finally, correlations between the intra-species excitations of A and B are shown to (dis-)favour certain energy transport channels depending, in general, on the direction of energy transfer (sect. 6.3).

Decoupled active subspaces
For a fixed total single-particle energy E SP (k) = k + (N A + 1)/2, k ∈ N 0 , we have evaluated the probability P SP (k; t) to find the bipartite system in the subspace spanned by configurations |n 0 , n 1 , ... A HO ⊗ |u B j with total single-particle energy E SP (k). So P SP (k; t) is obtained by summing over all P t (n 0 , n 1 , ...; j) with j + r≥1 r n r = k. In fig. 10 (a), we depict the probabilitiesP t (n 0 , n 1 , ...; j) for the configurations of total single-particle energy E SP (1), |N A , 0, ..  happens solely within the various subspaces of fixed total single-particle energy E SP (k), k ≥ 1 + , while inter-subspace transitions are suppressed.

Intra-species excitations
The joint probability (17) describing the distribution of excitations in the total system defines the following two marginal distributions for the A species and the single B atom, p A t (n 0 , n 1 , ..) = j P t (n 0 , n 1 , ...; j) = A HO n 0 , n 1 , ...|η A t |n 0 , n 1 , ... A HO , (18) In order to classify the collisionally induced excitations in the environment, we moreover inspect the probabilities for having no (p A;0 t ), a singlet (p A;1 t ) and a doublet excitation where n 0 = (N A , 0, 0, ...) and n i 0 ( n ij 0 ) refers to occupation number vectors with N A − 1 (N A − 2) atoms in the harmonic oscillator ground state and one (two) atom(s) in the i-th (i-th and j-th) excited state(s). In fig. 9, we compare these classes of A species excitations for N A = 2, 4, 10: One can clearly see that the excitations attain a maximum (minimum) whenever∆ B (t) (cf. fig. 2) becomes minimal (maximal). Moreover, singlet excitations dominate over doublets (and also not shown higher order excitations), which occur only after a certain delay with respect to the appearance of singlets. As expected from energetical considerations, the singlet excitations involve dominantly the |N A − 1, 1, 0, ... A HO configuration and second dominantly |N A − 1, 0, 1, 0... A HO , which can occur as likely as the dominant doublet contribution |N A − 2, 2, 0, ... A HO . In sect. 3, the increasing number of collision partners has been argued to be responsible for the acceleration of the energy transfer with increasing N A . Fig. 9 clearly shows this enhancement of excitation probability with increasing N A during the first few B oscillations (t 80). The larger N A , the faster does the probability for singlet excitations increase. For N A = 10, this effect is less visible since the intraspecies interaction results in a finite initial singlet probability. Concerning the long-time dynamics, we observe that the maximal singlet / doublet excitation probability reduces for larger environment sizes, which goes hand in hand with the energy donation to A becoming less efficient.
For obtaining analytical insights into the delayed doublet emergence and the decoupled active subspaces, we have applied a stroboscopic version of time-dependent perturbation theory to the simplest case of g A = 0. The total wave function |Ψ r after r oscillations is propagated for one free oscillation period T = 2π with linear perturbation theory with respect toĤ I AB (t) (with I denoting the interaction picture): where the time-averaged coupling Hamiltonian reads The energy conservation enforced by the Kronecker δ i+j,q+p , being reminiscent of Fermi's golden rule, is an immediate consequence of the equidistant single-particle spectra and α = 1. In order to avoid norm loss, we raise the above map to a unitary one, |Ψ I r+1 ≈ exp(−i 2πH I AB )|Ψ I r , which is equivalent to solving a temporally coarse-grained Schrödinger equation with the effective HamiltonianH I AB . We have observed good agreement between full numerical simulations (with g A = 0) and the perturbation theory (plots not shown) concerning the observables∆ B (t) and p A;i (t) during the first 18 (12) collisions for N A = 2 (N A = 7). Yet the entanglement entropy turns out to be overestimated a bit earlier by the perturbative approach. Comparing our full numerical data for g A = 0 and g A = g AB , we find qualitative agreement during the first r = 10, ..., 20 collisions such that the following insights from the perturbation theory should also apply in the presence of intra-species interaction. Firstly, one can easily show thatH I AB commutes with the projector onto the subspace of configurations with fixed E SP (k). Due to the detailed energy conservation in binary collisions, the temporally coarse-grained dynamics decouples thus subspaces of different total single-particle energy. Secondly,H I AB involves only one-body processes within the A species. Therefore, doublets can only emerge after the second collision with -as second order processes w.r.t. toH I AB -smaller amplitude than singlets, which explains the suppression and delay of the doublet creation. It goes without saying that the perturbative approach breaks down after fewer collisions for larger N A since (i) the intra-species interaction becomes more relevant and (ii) the enhanced inter-species coupling g AB N A leads to less separated time scales making a temporal coarse-graining questionable.

Correlations between intra-species excitations
Ultimately, we aim at unravelling correlations between excitations of the B atom and the A species. For this purpose, we introduce the following correlation measure: which is reminiscent of the diagonal of higher order coherence functions introduced by Glauber [46]. While g t (n 0 , n 1 , ...; j) = 1 means that the intra-species excitations take place statistically independently, a value larger (smaller) than unity implies that the corresponding joint excitation event (n 0 , n 1 , ...; j) is measured with an enhanced (decreased) probability, i.e. (anti-)bunches. In fig. 10, we present the evolution of both the joint probabilityP t (n 0 , n 1 , ...; j) and the Fock space correlation measurē g t (n 0 , n 1 , ...; j) for N A = 4 and selected configurations, which have a significant contribution to the total wave function. We emphasize that we do not intend to present a comprehensive overview over all energy transfer channels but rather concentrate on some important ones.
As expected, the considered excitations turn out to be uncorrelated wheneverS vN (t) attains a deep minimum, which is partly related to extremal values of∆ B (t) as discussed in sect. 5 (2) subspace. We conclude that, depending on the direction of the energy transfer, inter-species correlations stochastically favour or disfavour certain subsystem configurations and thereby energy transfer channels.

Conclusions and outlook
In this work, we have analysed the correlated energy transfer between a single atom and a finite, interacting bosonic environment mediated by a spatio-temporally localized coupling, focusing on an initial excess energy of the single atom just above the environment excitation threshold. By numerically exact simulations of the whole system, we have shown that the energy transfer happens the faster the larger the number of environmental degrees of freedom N A is, which is explained classically by an increased number of collision partners and quantum-mechanically by a with N A increasing level splitting of the degenerate two first excited states of the uncoupled system. At the same time, the energy transfer between the subsystems becomes less complete when increasing N A so that the excess energy distribution oscillates around quite a balanced one for larger N A , which might be regarded as a precursor of temporal relaxation induced by dephasing and has important consequences for inter-species correlations: While correlations are analytically shown to be suppressed whenever the excess energy distribution is maximally imbalanced, more balanced distributions allow for stronger correlations, which are dynamically established indeed. By inspecting the energy distribution among the natural orbitals of the B atom as well as the whole A species, incoherent transfer processes being intimately connected to inter-species correlations are unravelled. While donating (accepting) energy, each subsystem simultaneously shuffles a significant fraction of its energy into the second dominant natural orbital and back. These incoherent processes accelerate the early energy donation of the B atom and might be experimentally detectable by correlating the outcomes of single shot energy measurements upon both subsystems. Moreover, by analysing correlations between subsystem excitations, we show that inter-species correlation (dis-)favour certain energy transfer channels depending on the instantaneous direction of transfer.
The correlated energy transfer dynamics manifests itself in the subsystem dynamics as follows: The B atoms experiences a temporal loss of coherence in terms of short phases during which the Husimi function undergoes significant shape changes. These deviations from a coherent state are more pronounced for larger N A . The bosonic environment features depletion oscillations with a maximal depletion decreasing faster with increasing N A than its upper, energetically allowed bound ∝ 1/N A , which is derived for negligible intra-species interaction. In the considered low-excess energy regime, the excitations of the environment induced by the coupling to the B atom mainly consist of singlet and delayed doublet excitations, whose delay is explained by a temporally coarse-grained Schrödinger equation.
As an outlook, our work on systems with spatio-temporally localized couplings allows various interesting extensions: (i) Having understood the low-excess energy regime, a natural step consists in a careful increase of the excess energy in order to systematically open more and more energy transfer channels. (ii) It is important to address the questions whether the excess energy distribution among the subsystems temporally relaxes due to dephasing and whether inter-species correlations are persistent when significantly increasing N A to values O(100). However, such investigations will require to either separate the traps further or switch off / rescale the intra-species interaction strength in order to keep the initial inter-species overlap vanishing. (iii) The demonstrated tunability of the inter-species coupling matrix elements (5) in terms of the trap separation R should be exploited to control the correlated energy transfer.
Finally, we briefly comment on an experimental realization of this open quantum system problem. Having loaded an ensemble of |F = 1, m F = −1 polarized 87 Rb atoms in an optical dipole trap with a deep transverse optical lattice, the site-selective spin flip technique based on an additional longitudinal pinning lattice [13] or the doping technique in [15] allows for creating a single |F = 1, m F = 1 impurity. By adiabatically ramping up a linear longitudinal magnetic field gradient, the species are spatially displaced in opposite directions. Then the m F = 1 to m F = 0 transition can be selectively addressed by a RF field using the quadratic Zeeman effect such that the B atom (m F = 0) is initialized in a displaced coherent state, realizing d = R and α = 1. ∆ B t can be inferred from the B oscillation turning points via in situ density measurements [11,12] or tunnelling measurements [47].

Achnowledgements
The authors would like to thank Oriol Vendrell, Hans-Dieter Meyer, Walter Strunz, Lushuai Cao and Antonio Negretti for inspiring discussions, Christof Weitenberg for valuable support concerning the experimental realizability of this system as well as Maxim Efremov for his expertise in phase-space analysis. S.K. acknowledges a scholarship of the Studienstiftung des deutschen Volkes and P.S. acknowledges funding by the Deutsche Forschungsgemeinschaft in the framework of the SFB 925 "Light induced dynamics and control of correlated quantum systems".
One key feature of the ML-MCTDHB ansatz lies in its flexibility (iii): It covers both the mean-field limit M = m σ = 1 and the numerical exact limit. The latter is achieved by considering as many single-particle particle basis functions as one employs time-independent basis states, e.g. grid points, in order to represent them and choosing M = (N σ − m σ − 1)!/[N σ !(m σ − 1)!]. Moreover, in the presence of quite some intra-species but weak inter-species correlations, as it is often the case for systemenvironment problems, one might get converged results with M being much smaller than (N σ − m σ − 1)!/[N σ !(m σ − 1)!]. In our particular situation with N B = 1, one has to choose m B = M . So we have to ensure convergence with respect to numbers of grid points n and time-dependent basis functions, i.e. m A and M . By setting the number of time-dependent optimized species layer basis functions M = 1 while still bringing the simulation to convergence with respect to the number m A of time-dependent optimized A boson single-particle basis functions, the ML-MCTDHB equations of motion give a variationally optimized solution for the many-body problem under the constraint of vanishing inter-species correlations, which can illuminate the role of inter-species correlations as discussed in sect. 5.2.

Appendix B. Initial state preparation and convergence
In order to prepare our initial state, we first shift the harmonic trap for the single atom to the right by replacing (x B ) 2 /2 by (x B − d) 2 /2 in H B and propagate an initial guess for (A.2) with the ML-MCTDHB equations of motion in imaginary time. Thereby, we relax the binary mixture to its many-body ground state. Due to the almost vanishingly small spatial overlap of the two species, the result of this procedure is a many-body wave function factorizing to good approximation into the ground state of N A interacting bosons in their harmonic trap and the single atom being in the coherent state |z = d/ √ 2 . Afterwards, we propagate this state with the unshifted Hamiltonian (1).
For representing the time-dependent single-particle basis functions, we employ a harmonic oscillator discrete variable representation [52,53] with n = 401 grid points. Concerning the size of the time-dependent basis, we ensured convergence by checking that incrementing m A or M does not affect physical observables such as energy and the natural populations on the species and particle layers. For our purposes, m A = 3 and M = 4 turned out to lead to convergence for N A ≤ 10 and d = 1.5. Higher displacements d ≥ 2.5 lead to results which are still variationally optimal but of lower accuracy indicated e.g. by a small but non-vanishing least dominant NP λ A/B i .

Appendix C. Bound on the A species depletion
Neglecting the intra-species interaction, we show here that the depletion of the A species vanishes as 1/N A in the large N A limit as a consequence of having a fixed excess energy and a gapped excitation spectrumĤ A . Setting g A = 0 and assuming the initial inter-species interaction energy to be negligible Ĥ AB 0 Ĥ 0 (which is the case to a very good approximation in fact), energy conservation implies: given the excess energy ε = d 2 /2. SinceĤ A contains only the single-particle Hamiltoniansĥ A i = [(p A i ) 2 + (x A i + R) 2 ]/2 by assumption, we may express Ĥ A (t) as a function(al) of theρ A t NPs / NOs sorted in increasing sequence with respect to their instantaneous excess energyε A i (t) = φ A i (t)|ĥ A 1 |φ A i (t) − E A 0 being indicated by a tilde and obtain:ε From now on, we omit the tilde and the time-dependencies in the notation. The excess energy ε A 1 of the energetically lowest lying NO is thus bounded from above by ε/N A as a consequence of Ĥ A being extensive. Assuming now N A ε/(E A 1 − E A 0 ), similar arguments concerning the orthonormality of the NOs as explicated in appendix Appendix D apply, showing ε A i E A 1 − E A 0 = 1 ε/N A for i > 1. Finally, we may employ (C.2) in order to prove: Thus, the A species becomes condensed as N A → ∞.
Appendix D. Absence of correlations for strongly imbalanced excess energy distribution between subsystems In this appendix, we sketch the proof for the fact that the excess energy ε B t of the B atom being much smaller than the excitation gap (E B 1 − E B 0 ) implies the absence of significant inter-species correlations: Reordering the NOs in an increasing sequence with respect to their energy content Σ B i (t), one obtains the representation: From now on, we omit the tilde denoting this particular ordering as well as all timedependencies in the notation. With this expression, the excess energy ε B of B can be decomposed into the excess energies of the NOs obeying the ordering ε B i+1 ≥ ε B i : which leads to ε B 1 ≤ ε B . Furthermore, we will show ε B i > ε B 1 for all i ≥ 2 below implying an upper bound for the i-th NP: Expanding |ϕ B i in terms of the HO eigenstates |u B 0 , |u B 1 and a normalized vector |v i ⊥ orthogonal to the former ones, according to the Ritz variational principle. This identity can be employed to see that

Appendix E. Persistence of inter-species correlations under N A increase
Neglecting complications due to the intra-species interaction, we have shown in Appendix C that the A species becomes condensed as N A → ∞. This might seem to contradict the apparently persistent inter-species correlations discussed in sect. 5.1 since λ A 1 = 1 would necessarily imply λ B 1 = 1 given that the total system is in a pure state. In order to disprove this seeming contradiction, we provide a minimal example illustrating that inter-species correlations can survive the N A → ∞ limit: Let us assume that the state of the A species is given by: where the corresponding basis vectors are given by the following bosonic number states with |u A 0 and |u A 1 as the underlying single-particle states: |N A , 0 A

HO
(1, 0) T and |N A − 1, 1 A HO (0, 1) T . Equation (E.1) defines a density operator for all a ∈ R, b ∈ C such that the Bloch vector norm obeys |n| 2 ≡ a 2 + |b| 2 ≤ 1 and is consistent with the typical energies of the A species considered in this work. The corresponding NPs are given by λ B 1/2 = (1 ± |n|)/2 and, thus, any inter-species entanglement entropy S vN ∈ [0, ln 2] can be realized by appropriately choosing |n|. For the above stateη A , the depletion of the A species, vanishes in the large N A limit independently of the inter-species correlations.