Non-equilibrium quantum dynamics of ultra-cold atomic mixtures: the multi-layer multi-configuration time-dependent Hartree method for bosons

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 to describe intra- and inter-species correlations accurately and efficiently. To demonstrate the beneficial scaling and efficiency of the method, we explored the correlated tunneling dynamics of two species with repulsive intra- and inter-species interactions, to which a third species with vanishing intra-species interaction was 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 second species exhibit a bunching tendency, whose strength can be influenced by their coupling to the third species.


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 growing interest in the 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 [8,[10][11][12] and quantum-phase transitions [13]. As the mean-field theory becomes exact for weak interactions and large particle numbers [3,4,14], 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 2-60 atoms [13]. Therefore, the transition from few-to many-body behavior is, in particular for strongly correlated quantum dynamics, a subject of immediate interest.
Moreover it is meanwhile also experimentally routinely achievable to trap and manipulate different components or species 4 , 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], spincharge 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 this work is to develop a broadly applicable and efficient ab initio method for the quantum dynamics of such mixtures to explore the fundamental dynamical processes in trapped ultra-cold multi-species setups and 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 exponential 5 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, which allows this scaling to be softened, 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. [28,40]). Later MCTDH theory was 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 (MCTDH for bosons (MCTDHB)) [46,47]. Furthermore, 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 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 manybody system into account. We call this method ML-MCTDHB. The multi-layer structure of our many-body wave function ansatz allows us to adapt our many-body basis to systemspecific inter-and intra-species correlations, which leads to a beneficial scaling. Moreover, the bosonic exchange symmetry is directly employed for 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 interspecies interaction strengths. In particular, bosons of different kinds show a bunching tendency and the inter-species interaction strengths allow 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 behavior 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.

The multi-layer multi-configuration time-dependent Hartree boson method
Let us consider an ensemble of S bosonic species. In the ultra-cold regime, the interaction between neutral atoms can be modeled by a contact interaction [3,4,51]. 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 N σ bosons of species σ = 1, . . . , S readŝ Here,Ĥ σ denotes the one-body Hamiltonian of the species σ containing a general speciesdependent trapping potential U σ , andV σ andŴ σ σ refer to the intra-species interaction of species σ and to the inter-species interaction between σ and σ bosons, respectively: Note that the intra-and inter-species interaction strengths g σ , g σ σ have to be properly renormalized with respect to their values in three-dimensional (3D) space as a consequence of dimensional reduction [7]. We remark that the Hamiltonian may be explicitly time dependent for studying driven systems.

Wave function 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 | (t) , we employ a time-dependent, with the system comoving basis and restrict ourselves to the following class of ansatze. For each species σ , we take M σ timedependent orthonormal species states |ψ (σ ) i (t) (i = 1, . . . , M σ ), i.e. states of all the N σ 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: Each species state |ψ (σ ) i (t) refers to a system of N σ indistinguishable bosons and should therefore be expanded in terms of bosonic number states | n σ t : where we allow each σ boson to occupy m σ time-dependent single-particle functions (SPFs) |φ (σ ) j (t) , indicated by the time dependence of bosonic number states | n σ t . The integer vector n = (n 1 , . . . , n m σ ) contains the occupation number n j of the jth SPF such that all n j 's sum up to N σ , indicated by the symbol ' n|N σ ' in the summation.
Summarizing, our wave function ansatz consists of three layers. The expansion coefficients A i 1 ,...,i S (t) form the top layer. Then we have the C σ i; n (t)'s on the species layer, which allow the species states to move with the system, and finally, on the particle layer, the SPFs |φ (σ ) j (t) 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 |ψ (σ ) i (t) and the SPFs |φ (σ ) j (t) . This two-fold time dependence allows for significantly reducing the number of basis states, leading to a very efficient algorithm. Also note that our ML-MCTDHB approach to mixtures conceptually differs from ML-MCTDH in second-quantization representation [50] by the fact that we only employ three layers, one for the A-coefficients, one for the whole species and one for the single bosons, but allow for a time-dependent singleparticle basis.
Keeping the number of grid points for representing the SPFs fixed, the numbers of species states M σ and SPFs m σ serve as numerical control parameters: taking m σ to be equal to the number of grid points and M σ equal to the number of number state configurations, i.e. (5) and (6) proves to be numerically exact. Opposite to this full configuration interaction (CI) limit, the choice m σ = M σ = 1 leads to the mean-field or Gross-Pitaevskii approximation [3,4]. In between these two limiting cases, any choice with m σ equal to or smaller than the number of grid points and M σ ( 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 M σ = 1 but m σ > 1 might be sufficient.

Equations of motion
Our final task is to find appropriate equations of motion for the ansatz constituents A i 1 ,...,i S (t), C σ i; n (t) and |φ (σ ) j (t) , whose time dependence we will omit in the notation from now on. To find the variationally optimal wave function | (t) within our class of ansatzes for given M σ , m σ , 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 δ |(i ∂ t − H )| = 0 with |δ being a variation within our ansatz class (h ≡ 1) [53,54], which turns out to be equivalent to McLachlan's variational principle on our manifold of wave function ansatze [55].
The variation of the top layer coefficients A i 1 ,...,i S gives us the usual linear equation of motion known from matrix mechanics: where the Hamiltonian matrix with respect to Hartree products of species states becomes time dependent due to the coupling to the C σ i; n coefficients and to the SPFs. Its explicit form is given in appendix A.
Varying the species state expansion coefficients C σ i; n , we obtain the following equations of motion on the species layer: 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 [w σ σ ] coupling both SPFs and species states. The reduced density matrix of the species σ and the subsystem constituted by the species σ and σ (σ = σ ) enter (8) as [η 1,σ ] and [η 2,σ σ ], respectively (cf (B.2) and (B.3)). The orthonormality of the species states is ensured by the projectorP 1 Formulae for the above ingredients are given in appendices A and B, 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 nonlinear integro-differential equations: Here, [ρ 1,σ ] denotes the reduced density matrix of a σ boson and [ρ 2,σ σ ], [ρ 2,σ σ ] (σ = σ ) refer to the reduced two-body density matrix of two σ bosons, a σ and a σ boson, respectively (cf (C.3)-(C.5)).ĥ σ corresponds to the one-body Hamiltonianˆp 2 σ 2m σ + U σ (x σ ) and the intra-and inter-species interactions enter these equations of motion in the form of the mean-field operator matrices [v σ ] and [ŵ σ σ ], respectively. All these ingredients are explicated in appendix C. The s | again ensures the orthonormality of the SPFs. So we have arrived at a set of highly coupled evolution equations (7)- (9), whose general properties we analyze in the following section.

Properties of the multi-layer multi-configuration time-dependent Hartree boson
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 (single-particle) symmetry, ML-MCTDHB respects both the symmetry of the SPFs and 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. m σ equal to the number of grid points and M σ = (N σ + m σ − 1)!/[N σ !(m σ − 1)!], the projectors in (8) and (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 A-coefficients. The full CI limit, however, is numerically only manageable for extremely small particle numbers, whereas ML-MCTDHB being based on a smaller but optimally comoving with the system basis can treat much larger ensembles. In the opposite mean-field limit m σ = M σ = 1, the time dependence of the Aand the C-coefficients is given by trivial phase factors. With all the various reduced density matrices being equal to the c-number one, 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 projectorP 2;σ .
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 [ρ 1,σ ], [ρ 2,σ σ ] and [ρ 2,σ σ ]. The entropy of a species as well as correlations between two species can be deduced from [η 1,σ ] and [η 2,σ σ ], 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 M σ = (N σ + m σ − 1)!/[N σ !(m σ − 1)!], the ML-MCTDHB theory becomes equivalent to MCTDHB [46,47] and its generalization to mixtures [48], respectively. If, however, fewer species states are sufficient for a converged simulation, ML-MCTDHB proves to have a better scaling. With n being the number of grid points, one has to pay: complex coefficients for storing an ML-MCTDHB wave function, which should be compared with the costs for a corresponding MCTDHB expansion: For a detailed scaling comparison of the MCTDH-type methods, we refer to [56].

Application to correlated tunneling dynamics
Let us now explore the tunneling dynamics of three bosonic species, referred 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 m σ for all the bosons. Furthermore, each species shall consist of N σ = 6 bosons and shall experience the very same trapping potential made of a harmonic trap superimposed with a Gaussian at the trap center, i.e. U σ (x) ≡ U (x) = x 2 /2 + h/ √ 2πs 2 exp(−x 2 /2s 2 ) in harmonic oscillator unitsh = m σ = ω = 1. We choose h = 3 and s = 0.2 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 1.63 from the first excited band and its level spacing amounts to E ≈ 0.23, leading to a tunneling period of T = 2π/ E ≈ 27 for non-interacting particles. For the contact interaction strengths, we take g A (N A − 1) = 0.2, g B = 0.75 g A and g AB = 0.05 √ g A g B . Furthermore, the C bosons are assumed to have no intra-species interaction, i.e. g C = 0, but an attractive, vanishing or repulsive coupling to the bosons of species A and B: g XC ≡ g AC = g BC ∈ {−0.5 g AB , 0.0, 0.5 g AB }. Anticipating the results, we will show that this very weak interaction of strength g XC has a significant impact on the correlations 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, M σ ≡ M, and SPFs 6 , m σ ≡ m. 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 twobody 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.

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. g XC < 0, up to time t = 100. 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. (ii) g XC = 0 (green), (iii) g XC = 0.5 g AB (red), (iv) g XC = −0.5 g AB (blue).
A further phenomenon not captured in the mean-field picture is explained in figure 2: the probabilities of finding two bosons of the same species in the same well oscillate between 0.5 and 1.0 with the frequency 2/T in the mean-field calculation. In the many-body calculation, however, the probability of 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 remains oscillating between 0.5 and 1.0.
For discussing the convergence of the simulation, the ML-MCTDHB calculation for m = M = 3 is compared with the results for m = M = 4 in figures 1 and 2. The single-particle probabilities show 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(a), 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 ∼70, a third species state gains a weight more than 1%. Thus, much less species states than M = (N + m − 1)!/[N !(m − 1)!], i.e. the full CI limit on the species layer, are enough for a fair representation of the total wave function. Figure 3(b) 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 m = 3 to 4 just leads to a reshuffling of the third-highest natural population without affecting the results. The natural populations corresponding to a B boson show similar behavior 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 1.2% in the long-time propagation up to t = 300 (not shown). Note that the extra species layer is crucial for this convergence check. Our m = M = 4 simulation lasted roughly a week 7 , while a corresponding MCTDHB calculation would require 146 times more coefficients. 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 g AB > 0 with g XC = 0, g XC > 0 or g XC < 0 and g AB = g XC = 0. Although all the inter-species interaction strengths are much smaller than g A and g B , their concrete values have a strong influence on the tunneling dynamics: for no inter-species interactions is there a partial revival of the tunneling oscillation after a temporal equilibration. In the case of g AB > 0, 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 C and the other species does not result in a complete temporal equilibration to a probability of 0.5 but rather leads to a reduction of the amplitude of the oscillations around the equilibration value to 0.12. 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).

Long-time propagation and build-up of correlations
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 P(σ, L; σ , L) denote the probability of finding a σ and a σ boson in the left well and let P(σ, L) (P(σ , L)) be the probability for finding a σ (σ ) boson in the left well. Then the abovementioned correlation measure reads f L L (σ, σ ) := P(σ, L; σ , L)/[P(σ, L)P(σ , L)] and we similarly define f R R (σ, σ ) for the right well. Note that ( f L L , f R R ) is a straightforward extension of the diagonal elements of the g 2 coherence/correlation measure [60] to spatially discrete systems with distinguishable components. The dynamics of the centre of mass positions of the σ and the σ species have an impact on f L L and f R R , of course. To diminish this impact, we finally construct our correlation measure for finding a σ and a σ boson in the same well as f (σ, σ ) := ([ f L L (σ, σ ) 2 + f R R (σ, σ ) 2 ]/2) 1/2 . If the σ and the σ bosons tunnel independently, f (σ, σ ) will be unity. A value of f (σ, σ ) 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 f (A, B) of up to 40% above unity. For 75 t 225, 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 interspecies interactions with the C species, i.e. if g AB > 0 is the only non-vanishing inter-species interaction strength, the correlation measure f (A, B) lies mostly in between 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 back action 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 [η 1,C ], which monotonically increases up to 6.2% (g XC < 0) and 14% (g XC > 0), respectively.

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 work. Being based on an expansion in terms of permanents and on a multi-layer ansatz, our ML-MCTDHB method optimally and efficiently takes 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 SPFs 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 SPFs such that one obtains a variationally optimal representation of the many-body wave function at any instant in time. This allows us to achieve convergence with a much smaller basis than methods based on a time-independent basis. Moreover, if the inter-species interactions are not too strong, i.e. do not require consideration of as many species states as there are number state configurations for a given number of SPFs, 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 singlespecies state and one SPF, ML-MCTDHB reduces to coupled Gross-Pitaevskii equations.
Employing ML-MCTDHB for a tunneling scenario of three species, we have entered a parameter regime that is hardly accessible by other methods in such 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 twobody 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-onedimensional 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 SPFs, which allows for optimally describing bosons in quasi-one-or -two-dimensional traps embedded in 3D 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.