Non-Equilibrium Quantum Dynamics of Ultra-Cold Atomic Mixtures

# Non-Equilibrium Quantum Dynamics of Ultra-Cold Atomic Mixtures: the Multi-Layer Multi-Configuration Time-Dependent Hartree Method for Bosons

## Abstract

We develop and apply the multi-layer multi-configuration time-dependent Hartree method for bosons, which represents an ab initio method for investigating the non-equilibrium quantum dynamics of multi-species bosonic systems. Its multi-layer feature allows for tailoring the wave function ansatz in order to describe intra- and inter-species correlations accurately and efficiently. To demonstrate the beneficial scaling and the efficiency of the method, we explore the correlated tunneling dynamics of two species with repulsive intra- and inter-species interactions, to which a third species with vanishing intra-species interaction is weakly coupled. The population imbalances of the first two species can feature a temporal equilibration and their time-evolution significantly depends on the coupling to the third species. Bosons of the first and of the second species exhibit a bunching tendency, whose strength can be influenced by their coupling to the third species.

###### pacs:
03.75.Kk, 05.30.Jp, 03.65.-w, 31.15.-p

## 1 Introduction

Due to the high degree of controllability and isolatedness, trapped ultra-cold atoms serve as an ideal system for observing many-body quantum phenomena [1] and can even be employed to simulate quantum systems of quite a broad physical context [2]. In particular, there is a growing interest in the meanwhile accessible regime where a mean-field description [3, 4] as given by the Gross-Pitaevskii equation fails. Such states can be realized in e.g. optical lattices [5]. Feshbach [6] or confinement induced resonances [7, 8, 9] can be employed to tune the inter-atomic interaction strength. In particular, quasi-one-dimensional trapping geometries can enhance correlation effects in the strong interaction regime, leading to fascinating novel phases [10, 8, 11, 12] and quantum phase transitions [13]. As the mean-field theory becomes exact for weak interactions and large particle numbers [14, 3, 4], beyond mean-field physics can also be expected for small ensembles and finite interaction strengths. The latter regime is experimentally explored in e.g. arrays of decoupled one-dimensional tubes typically containing two to 60 atoms [13]. Therefore, the transition from the few- to many-body behaviour is, in particular for the strongly correlated quantum dynamics, a subject of immediate interest.

Moreover it is meanwhile also experimentally routinely achievable to trap and manipulate different components or species1, which allows for studying distinguishable subsystems with indistinguishable constituents. Such mixtures can be realized e.g. by preparing alkali atoms in different hyperfine states [15] or by trapping different elements [16]. Due to the interplay between the intra- and inter-species interaction strengths, these systems show a number of intriguing features such as phase separation [17] including symbiotic excitations like interlacing vortex lattices with mutually filled cores [18] and dark-bright solitons [19], spin-charge separation [20], various tunneling effects [21, 22, 23, 24, 25, 26, 27, 28], collective excitations [29, 30] and counterflow and paired superfluidity [31, 32, 33]. The purpose of the present work is to develop a broadly applicable and efficient ab initio method for the quantum dynamics of such mixtures in order to explore the fundamental dynamical processes in trapped ultra-cold multi-species setups and to study the few- to many-body transition.

Simulating the quantum dynamics of an interacting many-body system, however, is a tough task in general due to the exponential2 scaling of the state space with the number of particles. Besides e.g. the time-dependent density matrix renormalization group approach [34, 35, 36, 37], a promising concept to soften this scaling is based on a many-body wave function expansion with respect to a time-dependent, with the system comoving basis. This idea has been incorporated in the multi-configuration time-dependent Hartree method (MCTDH) [38, 39]. Being based on time-dependent Hartree products as the many-body basis, MCTDH is designed for distinguishable particles, but has also been applied to bosonic few-body systems (e.g. [40, 28]). Later the MCTDH theory has been generalized and extended in several ways: There is the multi-layer MCTDH (ML-MCTDH) method [41, 42, 43], which takes correlations between various subsystems into account and is thus particularly suitable for system-bath problems with distinguishable degrees of freedom (e.g. [44]). Taking the fermionic or bosonic particle exchange (anti-) symmetry in the time-dependent many-body basis into account, MCTDH has been specialized to treat larger fermionic (MCTDHF) [45] or bosonic systems (MCTDHB) [46, 47]. Furthermore, a direct extension of MCTDHB and MCTDHF to treat bose-bose, bose-fermi and fermi-fermi mixtures has been developed [48], including the possibility of particle-conversions [49]. An alternative approach to systems of indistinguishable particles is the so-called ML-MCTDH method in second-quantization representation, which employs the factorization of the many-body Hilbert space into a direct product of Fock spaces [50].

In this work, we derive and apply a novel ab initio approach to the non-equilibrium dynamics of ultra-cold correlated bosonic mixtures, which takes all correlations of the many-body system into account. We call this method the multi-layer multi-configuration time-dependent Hartree method for bosons (ML-MCTDHB). The multi-layer structure of our many-body wave function ansatz allows us to adapt our many-body basis to system specific inter- and intra-species correlations, which leads to a beneficial scaling. Moreover, the bosonic exchange symmetry is directly employed for an efficient treatment of the indistinguishable bosonic subsystems. We apply ML-MCTDHB to simulate the correlated tunneling dynamics of a mixture of three bosonic species in a double well trap. It is shown that the dynamics of the population imbalances of the species significantly differ for ultra-weak and vanishing inter-species interaction strengths. In particular, bosons of different kind show a bunching tendency and the inter-species interaction strengths allow for tuning these correlations up to a certain degree.

This paper is organized as follows: In section 2, the derivation and properties of the ML-MCTDHB method for bosonic mixtures are presented. ML-MCTDHB is then applied to simulate the complex tunneling behaviour of a mixture of three bosonic species in section 3. Finally, we summarize our results and embed the presented ML-MCTDHB theory for mixtures into a more general framework in section 4.

## 2 The ML-MCTDHB method

Let us consider an ensemble of bosonic species. In the ultra-cold regime, the interaction between neutral atoms can be modelled by a contact interaction [51, 3, 4]. For simplicity, we restrict ourselves to one-dimensional settings, which can be prepared by energetically freezing out the transversal degrees of freedom [3]. The Hamiltonian of such a mixture with bosons of species reads:

 ^H=S∑σ=1(^Hσ+^Vσ)+∑1≤σ<σ′≤S^Wσσ′. (1)

Here, denotes the one-body Hamiltonian of the species containing a in general species dependent trapping potential :

 ^Hσ=Nσ∑i=1(^pσi)22mσ+Uσ(^xσi)), (2)

and , refer to the intra-species interaction of species and to the inter-species interaction between and bosons, respectively:

 ^Vσ=gσ∑1≤i

Please note that the intra- and inter-species interaction strengths , have to be properly renormalized with respect to their 3d values as a consequence of dimensional reduction [7]. We remark that the Hamiltonian may be explicitly time-dependent for studying driven systems.

### 2.1 Wavefunction ansatz

The ML-MCTDHB method is an ab initio approach to the time-dependent Schrödinger equation for systems like (1). To reduce the number of basis states necessary for a fair representation of the total wave function , we employ a time-dependent, with the system comoving basis and restrict ourselves to the following class of ansatzes: For each species , we take time-dependent orthonormal species states (), i.e. states of all the bosons of species , into account. Due to the distinguishability of bosons of different species, the total wave function is expanded in terms of Hartree products of these many-body states:

 |Ψ(t)⟩=M1∑i1=1...MS∑iS=1Ai1,...,iS(t)|ψ(1)i1(t)⟩...|ψ(S)iS(t)⟩. (5)

Each species state refers to a system of indistinguishable bosons and should therefore be expanded in terms of bosonic number states :

 |ψ(σ)i(t)⟩=∑→n|NσCσi;→n(t)|→n⟩σt, (6)

where we allow each boson to occupy time-dependent single particle functions (SPFs) , indicated by the time-dependence of bosonic number states . The integer vector contains the occupation number of the -th SPFs such that all ’s sum up to , indicated by the symbol “” in the summation.

Summarizing, our wave function ansatz consists of three layers: The expansion coefficients form the top layer. Then we have the ’s on the species layer which allow the species states to move with the system and, finally on the particle layer, the SPFs allow for rotations of the single particle basis. It is crucial to notice that, in contrast to the standard method for solving the time-dependent Schrödinger equation by propagating expansion coefficients while keeping the basis time-independent, ML-MCTDHB is based on an expansion with respect to a comoving basis with a two-fold time-dependence in terms of the species states and the SPFs . This two-fold time-dependence allows for significantly reducing the number of basis states leading to a very efficient algorithm. Please also note that our ML-MCTDHB approach to mixtures conceptually differs from ML-MCTDH in second-quantization representation [50] by the facts that we only employ two layers, one for the whole species and one for the single bosons, but allow for a time-dependent single particle basis.

Having the number of grid points for representing the SPFs fixed, the numbers of species states and SPFs serve as numerical control parameters: Taking to be equal to the number of grid points and equal to the number of number state configurations, i.e. , the ansatz (5,6) proves to be numerically exact. Opposite to this full CI limit, the choice leads to the mean-field or Gross-Pitaevskii approximation [3, 4]. In between these two limiting case, any choice with equal to or smaller than the number of grid points and is possible, which allows us to adapt our ansatz to system specific intra- and inter-species correlations. If, for instance, the inter-species interactions are relatively weak compared to the intra-species interactions, a “species mean-field” ansatz with but might be sufficient.

### 2.2 Equations of motion

Our final task is to find appropriate equations of motion for the ansatz constituents , and , whose time dependence we will omit in the notation from now on. In order to find the variationally optimal wave function within our class of ansatzes for given , , we can employ the McLachlan variational principle, which enforces the minimization of the error of our equations of motion with respect to the exact Schrödinger equation [52]. In practice, however, it is easier to work with the Dirac-Frenkel variational principle with being a variation within our ansatz class () [53, 54], which turns out to be equivalent to McLachlan’s variational principle on our manifold of wave function ansatzes[55].

The variation of the top layer coefficients gives us the usual linear equation of motion known from matrix mechanics:

 i∂tAi1,...,iS=M1∑j1=1...MS∑jS=1⟨ψ(1)i1...ψ(S)iS|^H|ψ(1)j1...ψ(S)jS⟩Aj1,...,jS, (7)

where the Hamiltonian matrix with respect to Hartree products of species states becomes time-dependent due to the coupling to the coefficients and to the SPFs. Its explicit form is given in [56].

Varying the species state expansion coefficients , we obtain the following equations of motion on the species layer:

 i∂tCσi;→n = σ⟨→n|(\mathds1−^P1;σ)∑→m|Nσ(mσ∑j,k=1[hσ]jk^a†σj^aσk|→m⟩σCσi;→m+ (8) + 12mσ∑j,k,q,p=1[vσ]jkqp^a†σj^a†σk^aσq^aσp|→m⟩σCσi;→m+ + ∑σ′≠σMσ∑s,t=1Mσ′∑u,v=1mσ∑j,k=1[η−11,σ]is[η2,σσ′]sutv[wσσ′]jkuv^a†σj^aσk|→m⟩σCσt;→m),

where denotes the bosonic annihilation (creation) operator corresponding to the SPF , obeying the canonical commutation relations and . and represent the matrix elements of the one-body Hamiltonian and the intra-species interaction potential with respect to the SPFs, respectively. The inter-species interaction leads to the mean-field matrix coupling both SPFs and species states. The reduced density matrix of the species and the reduced density matrix of the subsystem constituted by the species and () enter (8) as and , respectively (cf. (15), (16)). The orthonormality of the species states is ensured by the projector . Formulas for the above ingredients are given in A and an efficient scheme for applying the annihilation and creation operators to the number states can be found in [56] (see also [57] in this context).

Finally, the variation of the SPFs leads to the following non-linear integro-differential equations:

 i∂t|ϕ(σ)i⟩=(\mathds1−^P2;σ) ( ^hσ|ϕ(σ)i⟩+ (9) + mσ∑j,k,q,p=1[ρ−11,σ]ij[ρ2,σσ]jkqp[^vσ]kq|ϕ(σ)p⟩+ + ∑σ′≠σmσ∑j,q=1mσ′∑k,p=1[ρ−11,σ]ij[ρ2,σσ′]jkqp[^wσσ′]kp|ϕ(σ)q⟩).

Here, denotes the reduced density matrix of a boson and , () refer to the reduced two-body density matrix of two bosons, a and a boson, respectively (cf. (19)-(21)). corresponds to the one-body Hamiltonian and the intra- and inter-species interactions enter these equations of motion in the form of the mean-field operator matrices and , respectively. All these ingredients are explicated in B. The projector again ensures the orthonormality of the SPFs. So we have arrived at a set of highly coupled evolution equations (7, 8, 9), whose general properties we analyse in the following section.

### 2.3 Properties of the ML-MCTDHB theory

Derived from the Dirac-Frenkel variational principle, the ML-MCTDHB evolution equations preserve both norm and energy [39]. Moreover, one can show that for a Hamiltonian with a (single particle) symmetry, ML-MCTDHB respects both the symmetry of the SPFs and the symmetry of the many-body state, given that initially the SPFs and the many-body state have a well-defined symmetry [56].

In the full CI limit, i.e. equal to the number of grid points and , the projectors in (8,9) turn into unit operators such that both the species states and the SPFs become time-independent. In this numerically exact limit, ML-MCTDHB becomes equivalent to the standard method of solving the time-dependent Schrödinger equation by propagating only the -coefficients. The full CI limit, however, is numerically only manageable for extremely small particle numbers, whereas ML-MCTDHB being based on a smaller but with the system optimally comoving basis can treat much larger ensembles. In the opposite mean-field limit , the time-dependence of the - and the -coefficients is given by trivial phase factors. With all the various reduced density matrices being equal to the c-number one, the equations (9) just differ from the coupled Gross-Pitaevskii equations of the mean-field theory for mixtures [3, 58] by a physically irrelevant phase factor as a consequence of the projector .

A converged ML-MCTDHB calculation takes all correlations into account. These correlations can be studied by means of reduced density matrices of various subsystems, which the ML-MCTDHB method provides for free. Single-particle coherence as well as correlations between two bosons of the same or of different species can be unravelled with the help of , and . The entropy of a species as well as correlations between two species can be deduced from and , for example. Moreover, an analysis of the natural populations and orbitals of various subsystems both serves as an internal convergence check (see below)[39] and can give physical insights [59].

In the case of just one species and in the full CI limit on the species level , the ML-MCTDHB theory becomes equivalent to MCTDHB [46, 47] and its generalization to mixtures [48], respectively. If, however, less species states are sufficient for a converged simulation, ML-MCTDHB proves to have a better scaling. With being the number of grid points, one has to pay:

 S∏σ=1Mσ+S∑σ=1(Mσ(Nσ+mσ−1mσ−1)+mσn) (10)

complex coefficients for storing a ML-MCTDHB wave function, which should be compared with the costs for a corresponding MCTDHB expansion:

 S∏σ=1(Nσ+mσ−1mσ−1)+S∑σ=1mσn. (11)

For a detailed scaling comparison of the MCTDH type methods, we refer to [56].

## 3 Application to correlated tunneling dynamics

Let us now explore the tunneling dynamics of three bosonic species, refered to as the A, B and C species in the following, in a double well trap. This setup both unravels interesting correlation effects and illustrates the beneficial scaling of ML-MCTDHB by introducing the extra species layer.

In the following, we assume that the three species are realized as different hyperfine states of an alkali element resulting in equal masses for all the bosons. Furthermore, each species shall consist of bosons and shall experience the very same trapping potential made of a harmonic trap superimposed with a Gaussian at the trap centre, i.e. in harmonic oscillator units . We choose and for the height and width of the barrier, respectively, which leads to three bands below the barrier, each consisting of two single particle eigenstates. The lowest band is separated by an energy difference of from the first excited band and its level spacing amounts to leading to a tunneling period of for non-interacting particles. For the contact interaction strengths, we take , and . Furthermore, the C bosons are assumed to have no intra-species interaction, i.e. , but an attractive, vanishing or repulsive coupling to the bosons of species A and B: . Anticipating the results, we will show that this very weak interaction of strength has a significant impact on the correlation between the A and B bosons.

As the particle numbers are the same for all species and because of the not too different interaction strengths, we provide for each species the same number of species states, , and SPFs3, . For preparing the initial state of the mixture, we block the right well by means of a high step function potential. All bosons are then put into the ground state of the resulting single particle Hamiltonian and, afterwards, we let the interacting many-body system relax to its ground state by propagating the ML-MCTDHB equations of motion in imaginary time. Ramping down the step function potential instantaneously, the resulting many-body state is finally propagated in real time in the original double well trap. Afterwards, we infer the probability of a particle to be in some well and the probability of finding two particles of the same or different species in the same well from the corresponding reduced one-body and two-body density matrices.

Here we would like to point out that we do not aim at an exhaustive study of this setup. Rather than showing a systematic parameter scan, we would like to present one striking example of multi-species non-equilibrium dynamics hardly being accessible in this precision by other methods but ML-MCTDHB, thereby illustrating the beneficial scaling and efficiency of the method. As we shall see, this setup shows very interesting correlation effects.

### 3.1 Short time tunneling dynamics

Let us firstly focus on the tunneling dynamics for an attractive coupling of the C bosons to the bosons of the A and B species, i.e. , up to time . From figure 1, we see that the A, B and C bosons exhibit Rabi tunneling with respect to the tunneling period on this time interval. The amplitude of the probability oscillations, however, decreases in the course of time for the A and the B bosons. This decrease can be interpreted as a temporal equilibration of the occupation probability of the left well as one can infer from the inset of figure 1 showing a somewhat lower accuracy long time propagation (see below). We also clearly see that the decrease of the probability amplitude is a genuine many-body property, not present in the mean-field description via coupled Gross-Pitaevskii equations. In contrast to this, the tunneling amplitude of the C bosons is not damped and its dynamics in the many-body calculation coincides with the mean-field description. This is a consequence of the vanishing intra-species and the very weak inter-species interaction strength.

A further phenomenon not captured in the mean-field picture is unravelled in figure 2: The probabilities for finding two bosons of the same species in the same well oscillate between 0.5 and 1.0 with the frequency in the mean-field calculation. In the many-body calculation, however, the probability for finding two A or B bosons in the same well features damped oscillations leading to a saturation of 0.73, which indicates a bunching tendency, while the probability of finding two C bosons stays oscillating between 0.5 and 1.0.

For discussing the convergence of the simulation, the ML-MCTDHB calculation for is compared with the results for in figures 1 and 2. The single particle probabilities show an excellent agreement. Only for the joint probability of finding two particles of the same or different (not shown) kind in the same well, there are marginal deviations. Hence, we can definitely regard the simulation as being converged. This judgement is also supported by the time evolution of the natural populations: From figure 3, we infer that most of the time only two natural orbitals significantly contribute to the reduced density matrix of the whole species A and, hence, to the total wave function. For times larger than , a third species state gains a weight more than 1%. Thus, much less species states than , i.e. the full CI limit on the species layer, are enough for a fair representation of the total wave function. Figure 3 shows that the initially fully condensed state of the A bosons evolves into a two-fold fragmented state. Increasing the number of particle SPFs from to just leads to a reshuffling of the third-highest natural population without affecting the results. The natural populations corresponding to a B boson show a similar behaviour due to the similar intra-species interaction strengths (not shown). In contrast to this, the C bosons stay in a condensed state and become depleted only by in the long time propagation up to (not shown). Please note that the extra species layer is crucial for this convergence check: Our simulation lasted roughly a week4, while a corresponding MCTDHB calculation would require to propagate 146 times more coefficients.

### 3.2 Long time propagation and build up of correlations

Now let us explore the tunneling on longer time scales with a somewhat lower accuracy calculation choosing and . A comparison with a , simulation shows only very small, quantitative deviations in the observables under consideration (plots not shown). In the inset of figure 1, the time-evolution of the probability to find an A boson in the left well is shown comparing four different situations, namely with , or and . Although all the inter-species interaction strengths are much smaller than and , their concrete values have a strong influence on the tunneling dynamics: For no inter-species interactions, there happens to be a partial revival of the tunneling oscillation after a temporal equilibration. In the case of , only for a vanishing or attractive coupling of the C species to the other species one can observe such a temporal equilibration state with subsequent partial tunneling revival of, however, smaller amplitude in comparison to the former case. A repulsive coupling between the C and the other species does not result in a complete temporal equilibration to a probability of but rather leads to a reduction of the amplitude of the oscillations around the equilibration value to . While the B bosons show a dynamics similar to the A bosons, the C bosons tunnel almost unaffected by the inter-species interactions and exhibit mean-field tunneling oscillations (plots not shown).

In order to measure the correlations between different species, we compare the conditional probability of finding a boson in e.g. the left well given that a boson has already been found there with the marginal probability of finding a boson in the left well: Let denote the probability for finding a and a boson in the left well and let () be the probability for finding a () boson in the left well. Then the above mentioned correlation measure reads and we similarly define for the right well. Please note that is a straightforward extension of the diagonal elements of the coherence / correlation measure [60] to spatially discrete systems with distinguishable components. The dynamics of the centre of mass positions of the and the species has an impact on and , of course. In order to diminish this impact, we finally construct our correlation measure for finding a and a boson in the same well as . If the and the bosons tunnel independently, will be unity. A value of greater (smaller) than one indicates an overall bunching (anti-bunching) tendency.

In figure 4, we find that a bunching tendency between an A and a B boson clearly builds up with a maximal correlation measure of up to above unity. For , this bunching tendency turns out to be most intense for the repulsive coupling of the C bosons to the other species and becomes least intense for an attractive coupling. In the absence of inter-species interactions with the C species, i.e. if is the only non-vanishing inter-species interaction strength, the correlation measure lies mostly inbetween these two curves. For a large fraction of the propagation time, the coupling of the C bosons to the other two species can thus control the inter-species correlations between species A and B up to a certain degree. Due to the fact that the C bosons approximately perform Rabi-oscillations with respect to the occupation probability of the left well, one might come to the conclusion that the C bosons provide a time-dependent potential for the other two species hardly experiencing a backaction on the considered time-scale. That this descriptive picture can only be approximately valid up to a certain time, can be inferred from the second largest natural population of , which monotonously increases up to () and (), respectively.

## 4 Conclusion and outlook

We have presented a novel ab initio method for simulating the non-equilibrium dynamics of mixtures of ultra-cold bosons. In particular, ML-MCTDHB is suitable for dealing with explicitly time-dependent systems, which will be explored in future works. Being based on an expansion in terms of permanents and on a multi-layer ansatz, our ML-MCTDHB method takes optimally and efficiently the bosonic exchange symmetry within each species into account and allows for adapting the ansatz to system specific intra- and inter-species correlations. Hereby, the numbers of provided single particle functions and species states serve as control parameters for ensuring convergence. For any choice of these numbers of basis functions, ML-MCTDHB rotates the species states and single particle functions such that one obtains a variationally optimal representation of the many-body wave function at any instant in time. This allows to achieve convergence with a much smaller basis than methods being based on a time-independent basis. Moreover, if the inter-species interactions are not too strong, i.e. do not require to consider as many species states as there are number state configurations for a given number of single particle functions, ML-MCTDHB proves to have a much better scaling than the best state-of-the-art method MCTDHB [46, 47, 48]. In the case of only a single species state and one single particle function, ML-MCTDHB reduces to coupled Gross-Pitaevskii equations.

Employing ML-MCTDHB for a tunneling scenario of three species, we have entered a parameter regime which is hardly accessible by other methods in such a controlled precision. Our simulations show that the imbalances of the populations can feature a temporal equilibration with subsequent revival of the population oscillations, where the duration of and the fluctuations around the equilibration state as well as the degree of completeness of the revival crucially depend on the inter-species interaction strengths. In our setup, we have furthermore found two-body bunching correlations between the first two species. The strength of this correlation can be tuned by a weak attractive or repulsive coupling of the third species - with no intra-species interaction - to the first two species without significantly altering the tunneling dynamics of that third species.

In this paper, ML-MCTDHB has been formulated for systems confined by quasi-one-dimensional traps and interacting via contact interactions. A direct generalization to arbitrary dimensions and interaction potentials is possible, of course. Moreover, it is also feasible to generalize ML-MCTDHB further by applying the multi-layering concept on the level of the single particle functions, which allows for optimally describing bosons in quasi-one- or -two-dimensional traps embedded in three-dimensional space with or without internal degrees of freedom [56]. Incorporating internal degrees of freedom on the level of SPFs then allows for taking particle converting interactions into account.

The authors would like to thank Hans-Dieter Meyer and Jan Stockhofe for fruitful discussions on MCTDH methods and symmetry conservation. Particularly, the authors would like to thank Jan Stockhofe for the DVR implementation of the ML-MCTDHB code. S.K. gratefully acknowledges financial support by the Studienstiftung des deutschen Volkes. L.C. and P.S. gratefully acknowledge funding by the Deutsche Forschungsgemeinschaft in the framework of the SFB 925 “Light induced dynamics and control of correlated quantum systems”.

## Appendix A Ingredients for the evolution equations of the species states

The matrix elements entering (8) read:

 [hσ]jk = ⟨ϕ(σ)j|[^p2σ2mσ+Uσ(^xσ)]|ϕ(σ)k⟩, (12) [vσ]jkqp = gσ⟨ϕ(σ)jϕ(σ)k|δ(^xσ1−^xσ2)|ϕ(σ)qϕ(σ)p⟩, (13) [wσσ′]jkuv = gσσ′∑→l|Nσ′−1mσ′∑q,p=1⟨ϕ(σ)jϕ(σ′)q|δ(^xσ1−^xσ′2)|ϕ(σ)kϕ(σ′)p⟩× (14) ×Q→l(q,p)(Cσ′u;→l+^q)∗Cσ′v;→l+^p,

where . “” refers to summation over all occupation numbers summing up to , and represents an occupation number vector with vanishing entries except for the -component being set to one. The reduced density matrix corresponding to the species can be calculated as:

 [η1,σ]is=∑Jσ(AJσi)∗AJσs, (15)

where the summation runs over all indices except for the index, which is fixed to be in the multi-index . For inversion, has to be regularized [39]. In analogy, the reduced density matrix of the subsystem constituted by the and species is given as:

 [η2,σσ′]sutv=∑Jσσ′(AJσσ′su)∗AJσσ′tv, (16)

where . Here, the summation runs over all indices except for the and index, which are fixed to be and in .

## Appendix B Ingredients for the evolution equations of the SPFs

In the particle layer, the mean-field operator matrices for the intra- and inter-species interaction are given as:

 [^vσ]kp = gσ∫\textupdx(ϕ(σ)k(x))∗ϕ(σ)p(x)δ(x−^xσ), (17) [^wσσ′]kp = gσσ′∫\textupdx(ϕ(σ′)k(x))∗ϕ(σ′)p(x)δ(x−^xσ). (18)

The one-body reduced density matrix of a boson, which also has to be regularized [39], can be calculated as:

 [ρ1,σ]ij=1NσMσ∑u,v=1[η1,σ]uv∑→l|Nσ−1Q→l(i,j)(Cσu;→l+^i)∗Cσv;→l+^j. (19)

For the reduced density matrices of two bosons and of a and a boson (), one has the following expressions:

 [ρ2,σσ]jkqp=1NσMσ∑u,v=1[η1,σ]uv∑→l|Nσ−2P→l(j,k)P→l(q,p)(Cσu;→l+^j+^k)∗Cσv;→l+^q+^p, (20) [ρ2,σσ′]jkqp=1NσMσ∑s,t=1Mσ′∑u,v=1[η2,σσ′]sutv∑→l|Nσ−1∑→m|Nσ′−1Q→l(j,q)Q→m(k,p)× (21) ×(Cσs;→l+^j)∗Cσt;→l+^q(Cσ′u;→m+^k)∗Cσ′v;→m+^p,

with and denoting the Kronecker delta function. We remark that the ML-MCTDHB code employs a different strategy than MCTDHB [61] for evaluating the various density matrices and the action of annihilation and creation operators on number states (cf. (8)) and refer to [56] for the details.

## References

### Footnotes

1. In the following, we will use the term species irrespectively of whether it refers to different elements, isotopes or internal states of an isotope.
2. For distinguishable particles. Indistinguishable particles lead to a binomial scaling.
3. The SPFs are represented by means of a harmonic discrete variable representation (DVR) [39].
4. For gridpoints on an Intel® Xeon® CPU E5530 with 2.40GHz.

### References

1. Bloch I, Dalibard J, and Zwerger W. Rev. Mod. Phys., 80:885, 2008.
2. Bloch I, Dalibard J, and Nascimbéne S. Nat. Phys., 8:267, 2012.
3. Pethick CJ and Smith H. Bose-Einstein Condensates in Dilute Gases. Cambridge University Press, 2nd edition, 2008.
4. Stringari S and Pitaevskii LP. Bose-Einstein Condensation. Oxford University Press, 2003.
5. Greiner M, Mandel O, Esslinger T, Hänsch TW, and Bloch I. Nature, 415:39, 2002.
6. Chin C, Grimm R, Julienne P, and Tiesinga E. Rev. Mod. Phys., 82:1225, 2010.
7. Olshanii M. Phys. Rev. Lett., 81:938, 1998.
8. Kinoshita T, Wenger T, and Weiss DS. Science, 305:1125, 2004.
9. Moritz H, Stöferle T, Günter K, Köhl M, and Esslinger T. Phys. Rev. Lett., 94:210401, 2005.
10. Girardeau M. J. Math. Phys., 1:516, 1960.
11. Paredes B, Widera A, Murg V, Mandel O, Fölling S, Cirac I, Shlyapnikov GV, Hänsch TW, and Bloch I. Nature, 429:277, 2004.
12. Haller E, Gustavsson M, Mark MJ, Danzl JG, Hart R, Pupillo G, and Nägerl H-C. Science, 325:1224, 2009.
13. Haller E, Hart R, Mark MJ, Danzl JG, Reichsöllner L, Gustavsson M, Dalmonte M, Pupillo G, and Nägerl H-C. Nature, 466:597, 2010.
14. Lieb EH and Seiringer R. Phys. Rev. Lett., 88:170409, 2002.
15. Myatt CJ, Burt EA, Ghrist RW, Cornell EA, and Wieman CE. Phys. Rev. Lett., 78:586, 1997.
16. Modugno G, Ferrari G, Roati G, Brecha RJ, Simoni A, and Inguscio M. Science, 294:1320, 2001.
17. Hall DS, Matthews MR, Ensher JR, Wieman CE, and Cornell EA. Phys. Rev. Lett., 81:1539, 1998.
18. Schweikhard V, Coddington I, Engels P, Tung S, and Cornell EA. Phys. Rev. Lett., 93:210403, 2004.
19. Becker C, Stellmer S, Soltan-Panahi P, Dörscher S, Baumert M, Richter E-M, Kronjäger J, Bongs K, and Sengstock K. Nat. Phys., 4:496, 2008.
20. Kleine A, Kollath C, McCulloch IP, Giamarchi T, and Schollwöck U. Phys. Rev. A, 77:013607, 2008.
21. Sun B and Pindzola MS. Phys. Rev. A, 80:033616, 2009.
22. Juliá-Díaz B, Guilleumas M, Lewenstein M, Polls A, and Sanpera A. Phys. Rev. A, 80:023616, 2009.
23. Satija II, Balakrishnan R, Naudus P, Heward J, Edwards M, and Clark CW. Phys. Rev. A, 79:033616, 2009.
24. Naddeo A and Citro R. J. Phys. B: At. Mol. Phys., 43:135302, 2010.
25. Pflanzer AC, Zöllner S, and Schmelcher P. J. Phys. B: At. Mol. Opt. Phys., 42:231002, 2009.
26. Pflanzer AC, Zöllner S, and Schmelcher P. Phys. Rev. A, 81:023612, 2010.
27. Chatterjee B, Brouzos I, Cao L, and Schmelcher P. Phys. Rev. A, 85:013611, 2012.
28. Cao L, Brouzos I, Chatterjee B, and Schmelcher P. New Journal of Physics, 14:093011, 2012.
29. Maddaloni P, Modugno M, Fort C, Minardi F, and Inguscio M. Phys. Rev. Lett., 85:2413, 2000.
30. Modugno M, Dalfovo F, Fort C, Maddaloni P, and Minardi F. Phys. Rev. A, 62:063607, 2000.
31. Kuklov AB and Svistunov BV. Phys. Rev. Lett., 90:100401, 2003.
32. Hu A, Mathey L, Danshita I, Tiesinga E, Williams CJ, and Clark CW. Phys. Rev. A, 80:023619, 2009.
33. Hu A, Mathey L, Tiesinga E, Danshita I, Williams CJ, and Clark CW. Phys. Rev. A, 84:041609, 2011.
34. White SR and Feiguin AE. Phys. Rev. Lett., 93:076401, 2004.
35. Daley AJ, Kollath C, Schollwöck U, and Vidal G. Journal of Statistical Mechanics: Theory and Experiment, page P04005, 2004.
36. Schollwöck U. J. Phys. Soc. Jpn., 74S:246, 2005.
37. Schollwöck U. Ann. Phys. (NY), 326:96, 2011.
38. Meyer H-D, Manthe U, and Cederbaum LS. Chem. Phys. Lett., 165:73, 1990.
39. Beck MH, Jäckle A, Worth GA, and Meyer H-D. Phys. Rep., 324:1, 2000.
40. Zöllner S, Meyer H-D, and Schmelcher P. Phys. Rev. Lett., 100:040401, 2008.
41. Wang H and Thoss M. J. Chem. Phys., 119:1289, 2003.
42. Manthe U. J. Chem. Phys., 128:164116, 2008.
43. Vendrell O and Meyer H-D. J. Chem. Phys., 134:044135, 2011.
44. Wang H and Shao J. J. Chem. Phys., 137:22A504, 2012.
45. Zanghellini J, Kitzler M, Fabian C, Brabec T, and Scrinzi A. Laser Phys., 13:1064, 2003.
46. Streltsov AI, Alon OE, and Cederbaum LS. Phys. Rev. Lett., 99:030402, 2007.
47. Alon OE, Streltsov AI, and Cederbaum LS. Phys. Rev. A, 77:033613, 2008.
48. Alon OE, Streltsov AI, Sakmann K, Lode AUJ, Grond J, and Cederbaum LS. Chem. Phys., 401:2, 2012.
49. Alon OE, Streltsov AI, and Cederbaum LS. Phys. Rev. A, 79:022503, 2009.
50. Wang H and Thoss M. J. Chem. Phys., 131:024114, 2009.
51. Huang K and Yang CN. Phys. Rev., 105:767, 1957.
52. McLachlan AD. Mol. Phys., 8:39, 1963.
53. Dirac PAM. Proc. Cambridge Philos. Soc., 26:376, 1930.
54. Frenkel J. Wave Mechanics. Clarendon Press, Oxford, 1934.
55. Broeckhove J, Lathouwers L, Kesteloot E, and Van Leuven P. Chem. Phys. Lett., 149:547, 1988.
56. Cao L, Krönke S, Vendrell O, and Schmelcher P. to be published.
57. Streltsov AI, Alon OE, and Cederbaum LS. Phys. Rev. A, 81:022124, 2010.
58. Kevrekidis PG, Frantzeskakis DJ, and Carretero-González R, editors. Emergent Nonlinear Phenomena in Bose-Einstein Condensates, volume 45 of Springer Series on Atomic, Optical, and Plasma Physics. Springer Berlin / Heidelberg, 2008.
59. Sakmann K, Streltsov AI, Alon OE, and Cederbaum LS. Phys. Rev. A, 78:023615, 2008.
60. Glauber RJ. Phys. Rev., 130:2529, 1963.
61. Streltsov AI, Alon OE, and Cederbaum LS. Phys. Rev. A, 81:022124, 2010.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters