A unified ab-initio approach to the correlated quantum dynamics of ultracold fermionic and bosonic mixtures

TL;DR

This paper introduces a unified computational method, ML-MCTDHX, for simulating the correlated quantum dynamics of ultracold mixtures of fermions and bosons, capturing intra- and inter-species correlations beyond mean-field approximations.

Contribution

The authors extend the ML-MCTDHB method to fermionic systems, creating a versatile approach for studying correlated dynamics in mixed quantum gases.

Findings

Demonstrated the method on Bose-Fermi and Fermi-Fermi collisions.

Showed correlations significantly influence transmission and reflection.

Identified correlation-driven effects beyond mean-field predictions.

Abstract

We extent the recently developed Multi-Layer Multi-Configuration Time-Dependent Hartree method for Bosons (ML-MCTDHB) for simulating the correlated quantum dynamics of bosonic mixtures to the fermionic sector and establish a unifying approach for the investigation of the correlated quantum dynamics of mixture of indistinguishable particles, be it fermions or bosons. Relying on a multi-layer wave-function expansion, the resulting Multi-Layer Multi-Configuration Time-Dependent Hartree method for Mixtures (ML-MCTDHX) can be adapted to efficiently resolve system-specific intra- and inter-species correlations. The versatility and efficiency of ML-MCTDHX is demonstrated by applying it to the problem of colliding few-atom mixtures of both Bose-Fermi and Fermi-Fermi type. Thereby, we elucidate the role of correlations in the transmission and reflection properties of the collisional events. In particular, we present examples where the reflection (transmission) at the other atomic species is a correlation-dominated effect, i.e.\ it is suppressed in the mean-field approximation.

Equations

$$\hat{H}=\sum _{\sigma =A,B}\left[{\hat{H}}_{\sigma }+{\hat{V}}_{\sigma }\right]+{\hat{W}}_{AB},$$

$${\hat{H}}_{\sigma }=\sum _{i=1}^{N_{\sigma }}\left[\frac{({\hat{p}}_i^{\sigma }{)}^2}{2{\mathcal{M}}_{\sigma }}+u_{\sigma }({\hat{x}}_i^{\sigma })\right]\equiv \sum _{i=1}^{N_{\sigma }}{\hat{h}}_i^{\sigma }.$$

$${\hat{V}}_{\sigma }=\sum _{1\le i<j\le N_{\sigma }}{v}_{\sigma }({\hat{x}}_i^{\sigma },{\hat{x}}_j^{\sigma }),$$

$${\hat{W}}_{AB}=\sum _{i=1}^{N_A}\sum _{j=1}^{N_B}{w}_{AB}({\hat{x}}_i^A,{\hat{x}}_j^B),$$

$$|\Psi (t)⟩=\sum _{i,j=1}^M{A}_{ij}(t)|{\psi }_i^A(t)⟩|{\psi }_j^B(t)⟩.$$

$$|\mathbf{n}{⟩}_t^{\sigma }=\frac{1}{\sqrt{n_1!n_2!...n_{m_{\sigma }}!}}({\hat{a}}_{1,\sigma }^{†}{)}^{n_1}({\hat{a}}_{2,\sigma }^{†}{)}^{n_2}...({\hat{a}}_{m_{\sigma },\sigma }^{†}{)}^{n_{m_{\sigma }}}|vac⟩,$$

$$\displaystyle\left[\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{k,\sigma},\hat{a}^{\dagger}_{q,\sigma}\right]_{\pm}$$

$${\left[{\hat{a}}_{k,\sigma }^{†},{\hat{a}}_{q,\sigma }^{†}\right]}_{\pm }$$

$$|{\psi }_i^{\sigma }(t)⟩=\sum _{\mathbf{n}|N_{\sigma }}{C}_{i;\mathbf{n}}^{\sigma }(t)|\mathbf{n}{⟩}_t^{\sigma },i=1,...,M$$

$$S=\int_{0}^{t}d\tau\,\langle\Psi(\tau)|\big{(}\hat{H}-i\partial_{\tau}\big{)}|\Psi(\tau)\rangle,$$

$$i{\partial }_t{A}_{ij}=\sum _{q,p=1}^M⟨{\psi }_i^A{\psi }_j^B|\hat{H}|{\psi }_q^A{\psi }_p^B⟩A_{qp}.$$

$$i\partial_{t}C^{\sigma}_{i;{\bm{n}}}={}^{\sigma}\langle{\bm{n}}|(1-\hat{P}^{\sigma}_{1})\left[\big{(}\hat{H}_{\sigma}+\hat{V}_{\sigma}\big{)}|\psi_{i}^{\sigma}\rangle+\sum_{l,j,k,p=1}^{M}[\eta^{-1}_{1,\sigma}]^{i}_{l}\,[\eta_{2,\sigma\bar{\sigma}}]^{lk}_{jp}\,[\hat{W}_{\sigma|\bar{\sigma}}]^{k}_{p}\,|\psi_{j}^{\sigma}\rangle\right],$$

$$\displaystyle i\partial_{t}|\phi_{i}^{\sigma}\rangle=(1-\hat{P}^{\sigma}_{2})\left[\hat{h}_{\sigma}|\phi_{i}^{\sigma}\rangle+\sum_{p,s=1}^{m_{\sigma}}[\rho^{-1}_{1,\sigma}]^{i}_{p}\Big{(}\sum_{q,l=1}^{m_{\sigma}}[\rho_{2,\sigma}]^{pq}_{sl}[\hat{v}_{\sigma}]^{q}_{l}+\sum_{q,l=1}^{m_{\bar{\sigma}}}[\rho_{2,\sigma\bar{\sigma}}]^{pq}_{sl}[\hat{w}_{\sigma|\bar{\sigma}}]^{q}_{l}\Big{)}|\phi_{s}^{\sigma}\rangle\right].$$

$${\hat{\eta }}_{1,\sigma }=\sum _{i=1}^M{\lambda }_i^{\sigma }|{\Psi }_i^{\sigma }⟩⟨{\Psi }_i^{\sigma }|,$$

$${\hat{\rho }}_{1,\sigma }=\sum _{i=1}^{m_{\sigma }}{n}_i^{\sigma }|{\Phi }_i^{\sigma }⟩⟨{\Phi }_i^{\sigma }|.$$

$$|\Psi ⟩=\sum _{j_1=1}^{M_1}...\sum _{j_S=1}^{M_S}{A}_{j_1...j_S}|{\psi }_{j_1}^1⟩\cdots |{\psi }_{j_S}^S⟩$$

$${\rho }_1^{\sigma }(x;t)=⟨\Psi (t)⟩{\hat{\Psi }}_{\sigma }^{†}(x){\hat{\Psi }}_{\sigma }(x)\Psi (t),$$

$${\Delta }_{C{C}^{\prime }}^{\sigma }(t)=\frac{1}{2N_{\sigma }}\int dx\left|{\rho }_{1,C}^{\sigma }(x;t)-{\rho }_{1,C^{\prime }}^{\sigma }(x;t)\right|\in [0,1],$$

$$\delta {\lambda }_i^{C{C}^{\prime }}(t)=\frac{\left|{\lambda }_i^C(t)-{\lambda }_i^{C^{\prime }}(t)\right|}{{\lambda }_i^{C^{\prime }}(t)},$$

$$|\Psi (t)⟩=\sum _{k=1}^M\sqrt{{\lambda }_k(t)}|{\Psi }_k^A(t)⟩|{\Psi }_k^B(t)⟩,$$

$${\rho }_1^{\sigma }(x;t)=\sum _{k=1}^M{\lambda }_k{\rho }_{1,k}^{\sigma }(x;t),$$

$$\displaystyle\langle\psi_{i_{\sigma}}^{\sigma}|\hat{a}^{\dagger}_{k,\sigma}\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{q,\sigma}|\psi_{j_{\sigma}}^{\sigma}\rangle$$

$$\displaystyle\langle\psi_{i_{\sigma}}^{\sigma}|\hat{a}^{\dagger}_{k,\sigma}\hat{a}^{\dagger}_{q,\sigma}\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{q^{\prime},\sigma}\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{k^{\prime},\sigma}|\psi_{j_{\sigma}}^{\sigma}\rangle$$

$$Q_{\mathbf{n}}(k,q)$$

$$P_{\mathbf{n}}(k,q)$$

$$Q_{\mathbf{n}}(k,q)$$

$$P_{\mathbf{n}}(k,q)$$

$$d_{\mathbf{n}}(k,q)=\sum _{a=\min (k,q)+1}^{\max (k,q)-1}{n}_a$$

$$\displaystyle\langle\psi_{i_{A}}^{A}\psi_{i_{B}}^{B}|\big{(}\hat{H}_{\sigma}+\hat{V}_{\sigma}\big{)}|\psi_{i_{A}^{\prime}}^{A}\psi_{i^{\prime}_{B}}^{B}\rangle$$

$$\displaystyle\phantom{=}+\frac{1}{2}\sum_{r,s,u,v=1}^{m_{\sigma}}[v_{\sigma}]_{uv}^{rs}\;\langle\psi_{i_{\sigma}}^{\sigma}|\hat{a}^{\dagger}_{r,\sigma}\hat{a}^{\dagger}_{s,\sigma}\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{v,\sigma}\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{u,\sigma}|\psi_{i_{\sigma}^{\prime}}^{\sigma}\rangle$$

$$⟨{\psi }_{i_A}^A{\psi }_{i_B}^B|{\hat{W}}_{AB}|{\psi }_{i_A^{\prime }}^A{\psi }_{i_B^{\prime }}^B⟩$$

$$[{\hat{W}}_{A|B}{]}_j^i$$

$$[{\hat{W}}_{B|A}{]}_j^i$$

$$[{\hat{v}}_{\sigma }{]}_l^q$$

Full text

I Introduction

Degenerate Bose and Fermi gases allow for an exquisite experimental control of the inter-atomic interactions as well as the motional and the internal quantum states which renders them ideal prototypes to investigate many-body quantum phenomena Bloch et al. (2008). This holds in particular for mixtures of different components which could be different chemical species or isotopes as well as different spin degrees of freedom represented by the hyperfine states. Depending on the statistics of the species, one can realize Bose-Bose (BB) Modugno et al. (2002), Fermi-Fermi (FF) Wille et al. (2008); Kohstall et al. (2012) and Bose-Fermi (BF) mixtures Schreck et al. (2001); Ospelkaus et al. (2006a); Heinze et al. (2011) with or without pseudo-spin degrees of freedom.

These ultracold mixtures reveal a great variety of physical phenomena depending among others on the particle statistics, the interaction strength and the external confinement. Examples for phenomena in BB systems are the relative phase evolution Anderson et al. (2009), composite fermionization and phase separation Garcia-March et al. (2013); Zöllner et al. (2008a), binary non-linear Kamchatnov et al. (2014) or collective Mertes et al. (2007) excitations, and the miscible to immiscible phase transition Ticknor (2013); Nicklas et al. (2015). For the corresponding spinor systems Kawaguchi and Ueda (2012); Stamper-Kurn and Ueda (2013) in optical lattices, their ferromagnetic and antiferromagnetic behavior Eisenberg and Lieb (2002); Sadler et al. (2006), the interplay of spin and charge excitations Zvonarev et al. (2007), governed by spin-charge separation in one dimension Kleine et al. (2008), and spin-exchange interactions Lee et al. (2009) have been investigated.

In nature, material properties are often determined by the underlying electronic structure and thus in particular fermionic degrees of freedom play a decisive role. By advanced cooling techniques Onofrio (2017), ultracold fermionic gases can nowadays be prepared with deterministic particle number and quantum state control Serwane et al. (2011) even for two-particle systems Murmann et al. (2015a) freeing the way to observe the onset of many-body physics by successively increasing the particle number Wenz et al. (2013). These systems exhibit fermion pairing Zürn et al. (2013) and Wigner molecule formation Brandt et al. (2015) and are ideally suited to investigate the physics of magnetism Murmann et al. (2015b); Yannouleas et al. (2016); Deuretzbacher et al. (2014); Volosniev et al. (2014), spin dynamics Koller et al. (2015); Volosniev et al. (2015, 2016), and collective excitations Koller et al. (2016). For two species mixtures of fermions a transition from a molecular bosonic condensate Strecker et al. (2003) to a Bardeen-Cooper-Schrieffer super-fluid of weakly bound Cooper pairs Stoof (1994); Zwierlein et al. (2004); Ries et al. (2015) has been observed. For strong inter-species interaction, even such distinguishable fermions can be fermionized Zürn et al. (2012).

In the same spirit, binary BF mixtures can be mapped to a FF system for strong inter-bosonic interactions Girardeau and Minguzzi (2007) being describable by the Gaudin-Yang model Gaudin (1967); Yang (1967). The low energy physics can be described by the Tomonaga-Luttinger liquid theory Cazalilla and Ho (2003) revealing a charge-density wave phase, a fermion pairing phase, and a phase separation regime Mathey et al. (2004) comparable to a liquid of polarons. Moreover, collective oscillations such as the monopole Maruyama et al. (2005), dipole Maruyama and Bertsch (2006), and scissors Roy et al. (2017) mode have been analyzed for BF mixtures. For attractive interactions, collapse Ospelkaus et al. (2006b) or the formation of composite particles, i.e. bound boson fermion pairs, becomes possible Varney et al. (2008); Burovski et al. (2009). Additionally, these systems are ideally suited to understand so-called induced interactions Kinnunen and Bruun (2015) and their impact on the super-fluid transition temperature Heiselberg et al. (2000).

In order to capture the above indicated variety of quantum phenomena, one needs a highly flexible theoretical description which is able to go beyond mean-field approximations such as the Gross-Pitaevskii (GP) Pitaevskii (1961); Gross (1963) or the Hartree-Fock (HF) Hartree (1928); Fock (1930) theory, and includes inter- and intra- species correlations. Since analytically solvable many-body problems are rare Lieb and Liniger (1963); Lieb (1963); Busch et al. (1998); Klaiman et al. (2017) and exact diagonalization is restricted to very few degrees of freedom, sophisticated methods to solve the time-dependent many-body Schrödinger equation of mixture systems from first principles are required. Various wave-function propagation methods face the exponential scaling of complexity with the number of particles by expanding the system state $|\Psi_{t}\rangle$ w.r.t. a dynamically optimized truncated basis. This family of methods may be divided into two classes with different underlying concepts of correlations. On the one hand, there are methods, in which correlations are regarded as deviations from a state that factorized w.r.t. to certain a-priori given modes, i.e. from a Gutzwiller mean-field state. Famous examples are the tensor network methods based on matrix-product or projected entangled-pair states Schollwöck (2005); Verstraete et al. (2008); Schollwöck (2011), which are mostly tailored to lattice problems and proved to be highly successful in particular for one-dimensional systems. On the other hand, correlations can be defined as deviations from a, possibly (anti-) symmetrized, product state w.r.t. particles. Here, the correlation-free reference state is a single Hartree product in the case of distinguishable particles or indistinguishable bosons and a single Slater determinant in the case of fermions, respectively. This is exactly the correlation concept underlying the Multi-Configuration Time-Dependent Hartree (MCTDH) family of methods where the total state $|\Psi_{t}\rangle$ is expanded w.r.t. a time-dependent many-body basis. Originally designed for distinguishable degrees of freedom Meyer et al. (1990); Beck et al. (2000), MCTDH became tailored to indistinguishable fermions and bosons in terms of the MCTDHF Zanghellini et al. (2003); Kato and Kono (2004); Nest et al. (2005) and MCTDHB method Streltsov et al. (2007); Alon et al. (2008), by using Slater determinants and permanents of variationally optimized time-dependent single-particle states as a many-body basis, respectively, which can formally be treated in a unifying manner Alon et al. (2007a, 2012). Recent developments cover the usage of Wannier basis states for lattice systems Cao et al. (2015); Lode and Bruder (2016), internal degrees of freedom Lode (2016), linear-response theory Grond et al. (2012); Alon et al. (2014) and restricting the active configuration space Miyagi and Madsen (2013); Lévêque and Madsen (2016). A direct extension of MCTDHB and MCTDHF to bosonic, fermionic or BF mixtures Alon et al. (2007b, 2012), however, is impractical except for very small system sizes due to the exponential wall of complexity. To pursue, the concept of the so-called Multi-Layer MCTDH method Wang and Thoss (2003); Manthe (2008); Vendrell and Meyer (2011), namely the coarse-graining of degrees of freedom, has been applied to bosonic mixtures in the recently developed Multi-Layer MCTDH method for Bosons (ML-MCTDHB) Krönke et al. (2013); Cao et al. (2013). Here the central idea is to expand the state of the bosonic mixtures w.r.t. dynamically optimized truncated basis states for the constituting species, which are in turn efficiently represented as MCTDHB wave functions. Thereby, the variational wave-function ansatz can be tailored to system-specific intra- and inter-species correlations. In passing, we note that the multi-layering idea has also been applied directly to the structure of Fock space in the so-called Multi-Layer MCTDH method in (optimized) Second Quantization Representation Wang and Thoss (2009); Manthe and Weike (2017), which is based on the concept of correlations between (optimized) modes.

The above variational approaches have already been applied successfully to ultracold systems. Beginning from one-dimensional single species systems, correlated processes in e.g. the tunneling Zöllner et al. (2008b); Sakmann et al. (2009); Lode et al. (2012); Schurer et al. (2016), breathing Schmitz et al. (2013), soliton Streltsov et al. (2011); Krönke and Schmelcher (2015), and quench induced dynamics in lattice systems Mistakidis et al. (2015); Beinke et al. (2016); Koutentakis et al. (2017) including periodic driving Mistakidis and Schmelcher (2017), hybrid atom-ion systems Schurer et al. (2015), and dipolar systems Cao et al. (2017); Streltsov (2013); Fischer et al. (2015) have been explored. In the realm of binary mixture systems, energy transfer Krönke et al. (2015), mixture tunneling Pflanzer et al. (2010); Cao et al. (2012), mesoscopic charged molecules Schurer et al. (2017), and dark-bright solitons Katsimiga et al. (2016) have been investigated. Meanwhile, MCTDH type methods have also been applied to two- and three-dimensional ultracold atomic ensembles Klaiman et al. (2014); Cao et al. (2015); Bolsinger et al. (2017); Wells et al. (2015); Weiner et al. (2017) for analyzing e.g. the dimensional cross-over or vortex dynamics.

In this work, we pursue the path by the ML-MCTDHB approach further in order to derive and apply an efficient method for simulating the correlated quantum dynamics of arbitrary ultracold atomic mixtures, let them be of bosonic, fermionic or Bose-Fermi type. For this purpose, the ML-MCTDHB wave-function ansatz is generalized to also include fermionic species and a unifying approach is applied, which does not distinguish between the particle-statistics. Using this ML-MCTDHX wave-function representation, we derive the equations of motion (EOM) for the expansion coefficients as well as the variational optimized basis states. Hereupon, we apply the ML-MCTDHX method to the quantum dynamics of colliding small binary mixtures of BF and FF type in one spatial dimension demonstrating the applicability of our ML-MCTDHX approach. Collisions of ultracold atomic clouds, being investigated theoretically Zin and Wasak (2017); Ögren and Kheruntsyan (2009); Deuar and Drummond (2007); Chwedeńczuk et al. (2006); Ziń et al. (2005); Bach et al. (2002) as well as experimentally Perrin et al. (2007); Bongs and Sengstock (2004); Vogels et al. (2002); Kasevich (2002), can be strongly influenced by correlations becoming observable in the transmission and reflection properties Deuar and Drummond (2007); Norrie et al. (2005, 2006). We reveal such correlation-induced transmission and reflection in a BF and a FF mixture, respectively. Moreover, we show that analyzing the many-body wave-function in terms of the natural orbitals of the species and the one-body density matrices allows for dedicated insights into the beyond mean-field dynamics.

This work is structured as follows. In section II, we derive the ML-MCTDHX method in detail. Starting with the wave-function ansatz, we compute and discuss the variational EOM in detail. Section III is devoted to the collisional dynamics of a BF and a FF mixture. After an exemplary detailed convergence study, we analyze the temporal evolution of the mixtures especially focusing on the impact of correlations. Finally, we provide our conclusion in section IV.

II Theory

In this section, we explicate the working principles of the ML-MCTDHX method for multi-component quantum gases by considering a binary mixture of indistinguishable particles. The various possibilities of having a BB, FF or BF mixture are treated within the same formalism and the particle-statistics turns out to enter only in the evaluation of certain entities. Discussing first the physical systems we aim to address, we then state the ML-MCTDHX wave-function ansatz and derive the corresponding equations of motion. Thereafter, we comment on limiting cases, convergence and the generalization of the binary mixture case to systems with more than two species. Finally, it is briefly addressed how we can treat atoms in more than one spatial dimensions and how we can incorporate internal degrees of freedom such as (pseudo-)spin into our approach.

II.1 Physical systems

In order to demonstrate the main idea of ML-MCTDHX, we consider a mixture of two species labeled with $A$ and $B$, respectively. Each species $\sigma=A,B$ shall be constituted by a fixed number of atoms $N_{\sigma}$ and can be either fermionic or bosonic. For simplicity, we first assume that the particles may only move in one-dimensional space and that their internal degrees of freedom, if existing, are frozen. The Hamiltonian of the binary mixture

[TABLE]

is decomposed into the single-particle and intra-species interaction terms of the $\sigma$ species, $\hat{H}_{\sigma}$ and $\hat{V}_{\sigma}$, respectively, as well as the inter-species interaction $\hat{W}_{AB}$. For being concrete, we consider a model with atoms in continuous space (treating discrete lattice systems is straightforward, see e.g. Cao et al. (2015); Lode and Bruder (2016)), where the single-particle Hamiltonian typically reads

[TABLE]

Here $\mathcal{M}_{\sigma}$ denotes the mass of a $\sigma$-species particle and $\hat{x}_{i}^{\sigma}$, $\hat{p}_{i}^{\sigma}$ refer to the position, momentum operator of the $i$-th $\sigma$-species particle, respectively. $u_{\sigma}(\hat{x})$ is an arbitrary, in general, species-selective external potential. In this paper, we are solely concerned with two-body interactions 111The extension to e.g. three-body terms is conceptually straightforward:

[TABLE]

where in the context of ultracold atoms the interaction potentials $v_{\sigma}(x_{1},x_{2})$, $w_{AB}(x_{1},x_{2})$ are usually of contact or dipolar type Stringari and Pitaevskii (2003); Pethick and Smith (2008). As ML-MCTDHX constitutes a wave-function propagation method, all terms of the Hamiltonian may be time-dependent.

As commonly done, we discretize the model at hand by assigning a typically large but finite number of $G_{\sigma}$ single-particle basis states to each particle of species $\sigma$, for instance a discrete-variable representation (DVR) Light and Carrington (2000); Beck et al. (2000) or, for lattice systems in the tight-binding approximation, Wannier states. These states span the single-particle Hilbert space of the $\sigma$-species $\mathfrak{h}^{\sigma}$. In the following, we assume that this discretization is fine enough to resolve the relevant physics such that we may refer to $\mathscr{H}=(\hat{S}_{A}\otimes_{i=1}^{N_{A}}\mathfrak{h}^{A}_{i})\bigotimes\,(\hat{S}_{B}\otimes_{j=1}^{N_{B}}\mathfrak{h}^{B}_{j})$ as the complete Hilbert space of the binary mixture. Here, $\mathfrak{h}^{\sigma}_{i}\equiv\mathfrak{h}^{\sigma}$ and $\hat{S}_{\sigma}$ denotes the (anti-)symmetrization operator if the species $\sigma$ is bosonic (fermionic). The complete Hilbert space $\mathscr{H}$ is typically far too high-dimensional for simulating quantum dynamics. In order to face this problem, the state of the total system is expanded with respect to a dynamically optimized, truncated many-body basis in ML-MCTDHX, whose construction is the subject of the subsequent section.

II.2 Construction of the many-body basis

In the ML-MCTDHX method, the state of the binary mixture $|\Psi(t)\rangle$ is firstly expanded as

[TABLE]

In the above expansion, most importantly, not only the coefficients $A_{ij}$ are time-dependent but also the states $|\psi_{i}^{\sigma}(t)\rangle$ themselves. The states $\{|\psi_{i}^{\sigma}(t)\rangle\}_{i=1}^{M}$ with $M\leq\dim{(\hat{S}_{\sigma}\otimes_{i=1}^{N_{\sigma}}\mathfrak{h}^{\sigma}_{i})}$ form a set of orthonormal functions for the $\sigma$ species, which can be viewed as a (truncated) orthonormal species basis 222To avoid redundant terms, the same upper summation bound $M$ can be chosen for both species (see Sec. II.5).. It is important to notice that these species basis states (SBSs) do not form a complete basis in the full Hilbert space of the $\sigma$ species $\hat{S}_{\sigma}\otimes_{i=1}^{N_{\sigma}}\mathfrak{h}^{\sigma}_{i}$ in general but only represent a basis in the time-dependent truncated sub-space being relevant for the expansion (5).

The SBSs are now constructed in the following way. First, we assign an orthonormal set of $m_{\sigma}$ time-dependent single-particle functions (SPFs) $\{|\phi_{n}^{\sigma}(t)\rangle\}_{n=1}^{m_{\sigma}}$ to each species, thereby truncating the single-particle Hilbert space $\mathfrak{h}^{\sigma}$ to a lower-dimensional, time-dependent sub-space. Second, given this time-dependent truncated single-particle space ${\rm span}\{|\phi_{n}^{\sigma}(t)\rangle\}_{n=1}^{m_{\sigma}}$, the correspondingly available truncated sub-space for the $N_{\sigma}$ particles of kind $\sigma$ is spanned by the occupation number states

[TABLE]

where $|vac\rangle$ is the vacuum state and $\hat{a}^{\dagger}_{i,\sigma}$ (\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{i,\sigma}) creates 333Note that the order of creation operators is important for fermions and fixed by definition (6). (annihilates) one particle in state $|\phi_{i}^{\sigma}(t)\rangle$. The second quantization operators obey the (anti-)commutation relations

[TABLE]

with $[\hat{A},\hat{B}]_{\pm}=\hat{A}\hat{B}\pm\hat{B}\hat{A}$ where the $+$ and the $-$ correspond to the case of the $\sigma$ particles being fermions or bosons, respectively. The vector ${\bm{n}}=(n_{1},n_{2},\cdots,n_{m_{\sigma}})$ with integers $n_{r}$ (for a fermionic species $n_{r}\in\{0,1\}$) determines the number of particles in a certain single-particle state with the restriction $\sum_{r=1}^{m_{\sigma}}n_{r}=N_{\sigma}$ such that $K_{\sigma}=\binom{N_{\sigma}+m_{\sigma}-1}{m_{\sigma}-1}$ [$K_{\sigma}=\binom{m_{\sigma}}{N_{\sigma}}$] different number states for $N_{\sigma}$ bosons (fermions) exist.

Finally, we truncate the available space ${\rm span}\{|{\bm{n}}\rangle^{\sigma}_{t}\}$ by considering only a reduced number $M\leq\min\{K_{A},K_{B}\}$ of optimally chosen time-dependent, pairwise orthonormal superpositions of $|{\bm{n}}\rangle^{\sigma}_{t}$

[TABLE]

where symbol ’${\bm{n}}|N_{\sigma}$’ indicates that the summation runs over all $K_{\sigma}$ possible number states for an $N_{\sigma}$-particle system. These states are the aforementioned SBSs, on which the expansion (5) is based.

In summary, the state of the system $|\Psi(t)\rangle$ is expanded with respect to the SBSs $\{|\psi_{i}^{\sigma}(t)\rangle\}_{i=1}^{M}$ and the SPFs $\{|\phi_{i}^{\sigma}(t)\rangle\}_{i=1}^{m_{\sigma}}$, which corresponds to a two-layer Hilbert space truncation: On the so-called particle layer, the assignment of a set of SPFs to each species truncates the single-particle Hilbert space $\mathfrak{h}_{n}^{\sigma}$ to a smaller sub-space of dimension $m_{\sigma}<G_{\sigma}$. On the so-called species layer, one takes only a reduced number of $M<\min\{K_{A},K_{B}\}$ time-dependent SBSs for each species into account instead of expanding the state of the binary mixture w.r.t. all conceivable number-state configurations of the two species $|{\bm{n_{A}}}\rangle^{A}_{t}|{\bm{n_{B}}}\rangle^{B}_{t}$ with the SPFs as the underlying single-particle states. Thereby, $m_{\sigma}$ and $M$ constitute the main numerical control parameters (see the discussion in Sect. II.4). In the following section, we explain how the optimal time-dependent basis states are found. For simplicity, we drop the time arguments in the notation in the remainder of this work.

II.3 Equations of motion

In order to derive equations of motion (EOM) for the expansion coefficients $A_{ij}$ and for the time-dependent basis states $|\phi_{i}^{\sigma}\rangle$, $|\psi_{i}^{\sigma}\rangle$, which ensure that these basis states move in an optimal manner, we employ the Langrangian variational principle, which turns out to be technically most convenient Alon et al. (2008); Broeckhove et al. (1988). Thus, we aim at finding the stationary point of the action functional (setting $\hbar=1$)

[TABLE]

under the constraints of conserving the normalization of the total state as well as the orthonormality of the SPFs and SBSs. Enforcing these constraints by Lagrange multipliers, however, does not give unique EOM, since expansions w.r.t. time-dependent basis states such as (5) are not unique: Time-dependence can be shuffled between the expansion coefficients and the time-dependent basis states by means of appropriate (time-dependent) unitary transformations of the expansion coefficients and the basis states leaving the total state invariant. This “gauge” degree of freedom can be fixed by setting the so-called constraint operators $\langle\phi_{q}^{\sigma}|i\partial_{t}|\phi_{k}^{\sigma}\rangle$ and $\langle\psi_{q}^{\sigma}|i\partial_{t}|\psi_{k}^{\sigma}\rangle$ to some hermitian matrices Beck et al. (2000), which we take to be zero in this work without loss of generality and without affecting the accuracy of the method. Now let us summarize the resulting EOM, while presenting the explicit formulas for the ingredients in Appendix A - D.

When performing the variation of $A^{*}_{ij}$, we obtain the following EOM

[TABLE]

As expected, the expansion coefficients obey a linear Schrödinger equation with a, by virtue of the SBSs, time-dependent Hamiltonian matrix. Varying correspondingly the coefficients $(C^{\sigma}_{i;{\bm{n}}})^{*}$ gives

[TABLE]

where $\bar{\sigma}=A(B)$ if $\sigma=B(A)$. The SBSs dynamics is driven by both the intra-species Hamiltonian $\hat{H}_{\sigma}+\hat{V}_{\sigma}$ and the coupling to the other species: the mean-field operator matrix $[\hat{W}_{\sigma|\bar{\sigma}}]^{k}_{p}=\langle\psi_{k}^{\bar{\sigma}}|\hat{W}_{AB}|\psi_{p}^{\bar{\sigma}}\rangle$, being an effective one-body operator, describes the action of the inter-species interaction on the $\sigma$ species conditioned on the transition of the $\bar{\sigma}$-species from the $p$-th to the $k$-th SBS. This operator is weighted with the corresponding two-species density matrix $[\eta_{2,\sigma\bar{\sigma}}]^{lk}_{jp}$ in the SBS representation and contracted with the inverse of the reduced density matrix 444As one can see in the spectral representation, the inverse of the reduced density matrix makes SPFs of low population rotate fast in oder to become dynamically optimized. Note that (near) singular density matrices have to be regularized (see Beck et al. (2000) and the recently proposed alternative to the regularization Manthe (2015)). of the $\sigma$ species in the SBS representation, $[\eta_{1,\sigma}]^{i}_{l}$. Summarizing, the inter-species interaction becomes manifest in a two-fold manner: it acts as an effective one-body operator on the $\sigma$ species but also leads to a coupling of $i\partial_{t}C^{\sigma}_{i;{\bm{n}}}$ to the expansion coefficients of all the other SBSs. The operator $1-\hat{P}^{\sigma}_{1}\equiv 1-\sum_{r=1}^{M}|\psi_{r}^{\sigma}\rangle\!\langle\psi_{r}^{\sigma}|$ denotes the projector onto the orthogonal complement of the space spanned by the SBSs w.r.t. ${\rm span}\{|{\bm{n}}\rangle^{\sigma}_{t}\}$ naturally preventing the SBSs from moving within the already spanned subspace ${\rm span}\{|\psi_{r}^{\sigma}\rangle\}_{r=1}^{M}$.

Finally varying the SPFs, we end up with

[TABLE]

Here again, the projector $1-\hat{P}^{\sigma}_{2}\equiv 1-\sum_{i=1}^{m_{\sigma}}|\phi_{i}^{\sigma}\rangle\!\langle\phi_{i}^{\sigma}|$ enforces that the SPFs may only rotate out of the subspace spanned by the instantaneous SPFs. Such a rotation can be induced by the single-particle Hamiltonian $\hat{h}_{\sigma}\equiv\hat{h}^{\sigma}_{1}$ and by the interaction with particles of the same (other) species via the mean-field operator matrices $[\hat{v}_{\sigma}]^{q}_{l}$ ($[\hat{w}_{\sigma|\bar{\sigma}}]^{q}_{l}$). These operators have potential-term character 555Given that the particles are interacting via interaction potentials $v_{\sigma}(\hat{x}_{1}^{\sigma},\hat{x}_{2}^{\sigma})$, $w_{AB}(\hat{x}_{1}^{A},\hat{x}_{1}^{B})$, i.e. local operators, as it is natural to assume. and are given by the interaction-potential operators conditioned on the interaction partner undergoing a transition from the $l$-th to the $q$-th SPF: $[\hat{v}_{\sigma}]^{q}_{l}=\langle\phi_{q}^{\sigma}|v_{\sigma}(\hat{x}_{1}^{\sigma},\hat{x}_{2}^{\sigma})|\phi_{l}^{\sigma}\rangle$ and $[\hat{w}_{\sigma|\bar{\sigma}}]^{q}_{l}=\langle\phi_{q}^{\bar{\sigma}}|w_{AB}(\hat{x}_{1}^{A},\hat{x}_{1}^{B})|\phi_{l}^{\bar{\sigma}}\rangle$. These transitions are weighted with the corresponding reduced two-body density matrix elements of two $\sigma$ particles (one $\sigma$ and one $\bar{\sigma}$ particle) in the SPF representation, namely $[\rho_{2,\sigma}]^{pq}_{sl}$ ($[\rho_{2,\sigma\bar{\sigma}}]^{pq}_{sl}$). Last, these terms are contracted with the inverse of the reduced one-body density matrix Note (4) of a $\sigma$ particle, $[\rho^{-1}_{1,\sigma}]^{i}_{p}$.

The coupled EOM (11), (12), and (13), which have to be solved simultaneously, formally appear alike the EOM of the ML-MCTDHB method for bosonic mixtures Krönke et al. (2013); Cao et al. (2013), which has been achieved by the employed unifying notation. As a matter of fact, the particle-statistics enter the ML-MCTDHX EOM only via the second quantization algebra when applying $\hat{H}_{\sigma}$, $\hat{V}_{\sigma}$ and $[\hat{W}_{\sigma|\bar{\sigma}}]^{k}_{p}$ to a SBS in (12) as well as in the calculation of the so-called reduced one- and two-body transition matrices Cioslowski (2000); Mazziotti (2007) \langle\psi_{k}^{\sigma}|\hat{a}^{\dagger}_{i,\sigma}\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{j,\sigma}|\psi_{l}^{\sigma}\rangle and \langle\psi_{k}^{\sigma}|\hat{a}^{\dagger}_{i,\sigma}\hat{a}^{\dagger}_{j,\sigma}\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{q,\sigma}\hat{a}^{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\dagger}}_{p,\sigma}|\psi_{l}^{\sigma}\rangle (see Appendix A). The latter entities measure the overlap between the many-body states $|\psi_{l}^{\sigma}\rangle$, $|\psi_{k}^{\sigma}\rangle$ after removing one (two) particles, i.e. probe one (two) body properties in the transitions between many-body states, and constitute essential building blocks for constructing the Hamiltonian matrix $\langle\psi_{i}^{A}\psi_{j}^{B}|\hat{H}|\psi_{q}^{A}\psi_{p}^{B}\rangle$, the mean-field operator matrices $[\hat{W}_{\sigma|\bar{\sigma}}]^{k}_{p}$ and the reduced density matrices $[\rho_{1,\sigma}]^{i}_{p}$, $[\rho_{2,\sigma}]^{pq}_{sl}$ and $[\rho_{2,\sigma\bar{\sigma}}]^{pq}_{sl}$ (see Appendix B - D). In this way, the reduced transition matrices are reminiscent of the reduced density operators in the unified formulation Alon et al. (2007b, 2012) of the MCTDH theory for mixtures of indistinguishable bosons and / or fermions, which, however, is limited to small particle numbers since the additional truncation to $M$ SBSs is not performed.

II.4 Limiting cases and convergence

Whenever the basis set on one layer is not truncated, the projector $\hat{P}^{\sigma}_{1}$ or $\hat{P}^{\sigma}_{2}$ in the corresponding EOM becomes the identity such that the corresponding basis states do not move (in the chosen constraint operator gauge). In particular, the ML-MCTDHX EOM reduce to the MCTDH theory for BB, FF and BF mixtures in its unifying formulation Alon et al. (2007b, 2012) if $M=\min\{K_{A},K_{B}\}$. When choosing $m_{\sigma}=G_{\sigma}$ in addition, our wave-function ansatz is of full configuration-interaction type, covering the complete Hilbert space $\mathscr{H}$. In the other extreme case of neglecting all correlations by setting $M=1$ and $m_{\sigma}=1$ for bosons ($m_{\sigma}=N_{\sigma}$ for fermions), the ML-MCTDHX EOM reduce to coupled GP (HF) EOM, depending on the statistics of the components. In between these extreme cases for not too strong inter-species correlations, one expects that it is sufficient to take only $M\ll\min\{K_{A},K_{B}\}$ variationally optimized SBSs into account, which leads to higher efficiency in the representation of the many-body state compared to the direct extension of the MCTDHF and MCTDHB methods for treating mixtures of indistinguishable degrees of freedom Alon et al. (2007b, 2012). Thus when using ML-MCTDHX, $m_{\sigma}$ and $M$ need to be chosen according to the dominant correlations in the system such that the expansion (5) is efficient, i.e. involves as few coefficients as possible, while capturing all relevant correlations. We remark that the EOM can be shown to conserve energy and norm Beck et al. (2000) as well as single-particle symmetries (see Cao et al. (2013) for details) independently of the chosen number of SPFs and SBSs.

In order to judge the convergence of a simulation, one has to systematically increase $m_{\sigma}$ and $M$, i.e. enlarge the subspace within which ML-MCTDHX finds the variationally optimal solution, and compare the results for observables of interest, e.g. the density distribution. Furthermore, one can inspect the spectral decomposition of the reduced density operator of the whole species $\sigma$

[TABLE]

and the reduced one-body density operator of the species $\sigma$

[TABLE]

The eigenstates $|\Psi^{\sigma}_{i}\rangle$ [$\,|\Phi^{\sigma}_{i}\rangle\,$] are called natural species functions (NSFs) [natural orbitals (NOs)] and the eigenvalues $\lambda_{i}^{\sigma}$ [$\,n_{i}^{\sigma}\,$], which sum up to unity, are denoted as natural populations (NPs). We label the NSFs and NOs in decreasing sequence w.r.t. their natural populations in the following and note that $\lambda_{i}^{A}=\lambda_{i}^{B}\equiv\lambda_{i}$ holds for binary mixtures (see Sect. II.5). Convergence can now be judged on the basis of the convergence of the individual natural populations. Thereby, a negligibly small lowest natural population $\lambda_{M}^{\sigma}$ ($n^{\sigma}_{m_{\sigma}}$) means that the corresponding NSF (NO) essentially does not contribute to the numerically obtained $|\Psi(t)\rangle$. The latter provides in practice a good indicator for the considered basis being sufficiently large, even though it is not a strict convergence criterion Cosme et al. (2016). We note that the convergence behavior and in particular the degree of convergence in general depend on the character of the observable of interest, of course. In Sect. III, we will illustrate how to control convergence of ML-MCTDHX simulations and, moreover, employ the NSFs and NOs as well as their populations as an analysis tool for a descriptive unraveling of correlated many-body states.

II.5 Generalizations

Having outlined the ML-MCTDHX wave-function ansatz and EOM for a one-dimensional binary mixture, we briefly comment on the possible generalizations, which we can also treat with our current implementation. First, more than two species can straightforwardly be taken into account by employing the ansatz

[TABLE]

with $S$ denoting the number of species. As a consequence, the calculation of the one- and two-species density matrix has to be adapted (see Appendix D) and the EOMs for the SBSs and SPFs of the $\sigma$ species feature $S-1$ inter-species mean-field operator matrix terms. While for $S>2$ the numbers of SBSs $M_{\sigma}$ for the different species do not have to coincide, we remark that taking unequal numbers $M_{A}$, $M_{B}$ for a binary mixture does not improve the accuracy of the wave-function ansatz. This is a consequence of the Schmidt decomposition for bipartite systems Nielsen and Chuang (2000); Horodecki et al. (2009), which states that the maximal strength of interspecies correlations, which can be resolved, only depends on $\min\{M_{A},M_{B}\}$ such that we may set $M_{A}=M_{B}=M$ without loss of generality. This feature is also reflected in the aforementioned coinciding natural populations of the NSFs, $\lambda_{i}^{A}=\lambda_{i}^{B}\equiv\lambda_{i}$.

Second, the concept of multi-layering can be applied to simulate also three-dimensional systems. Here, the idea is to expand the time-dependent three-dimensional SPFs $|\phi_{i}^{\sigma}\rangle$ w.r.t. Hartree products of time-dependent one-dimensional basis states for the different spatial directions. By adjusting the numbers of these one-dimensional basis states, correlations between the spatial directions can be taken into account. If the latter are not too strong, as in the physically interesting cross-over regime between one and three dimensions, this approach allows for efficient simulations on extremely large spatial grids Cao et al. (2013); Bolsinger et al. (2017, ). In a completely similar manner, we can account for internal degrees of freedom in ML-MCTDHX, allowing for simulating in particular interacting ultracold fermions in different pseudo-spin states including also spin-changing collisions Krauser et al. (2014). For this purpose, the SPFs $|\phi_{i}^{\sigma}\rangle$ are expanded w.r.t. Hartree products of time-dependent spatial wave-functions and time-independent spin basis states.

III Applications: Collisional dynamics of binary mixtures

In this section, we apply the ML-MCTDHX method for simulating the correlated collision dynamics of various binary mixtures in order to demonstrate the versatility and certain features of the approach. After discussing the general setup in Sect. III.1, we focus on the convergence behavior by considering a simple scenario involving a few-atom BF mixture (Sect. III.2). Thereafter, we proceed to more complex scenarios involving larger atom numbers and illustrate how manifestations of correlations can be unraveled by means of tailored analysis tools. Here, we consider the cases of a particle-number imbalanced BF (Sect. III.3) and a balanced FF mixture (Sect. III.4). We note that exploring the observed physical effects in depth clearly goes beyond the scope of this work.

III.1 Setup

We consider a binary mixture consisting of $A$ and $B$ species of the same mass 666For simplicity, we consider the constituents to possess the same mass also for the case of a BF mixture., being confined in harmonic-oscillator potentials of same frequencies. Working in harmonic oscillator units, the single-particle Hamiltonian takes the form $h_{i}^{\sigma}=-\frac{1}{2}\partial^{2}_{x_{i}^{\sigma}}+\frac{1}{2}(x_{i}^{\sigma}-x_{0}^{\sigma})^{2}$ where $x_{0}^{\sigma}$ (${\sigma}\in\{A,B\}$ is the species index) denotes the spatial offset added initially to the harmonic confinement for each of the species. The two-body intra-species and inter-species interaction potential is given by $v_{\sigma}(x_{i}^{\sigma},x_{j}^{\sigma})=g_{\sigma}\delta(x_{i}^{\sigma}-x_{j}^{\sigma})$ and $w_{AB}(x_{i}^{A},x_{j}^{B})=g_{AB}\delta(x_{i}^{A}-x_{j}^{B})$ respectively.

Below, we explore the dynamics of both a BF and a FF mixture with each of the individual species being considered as spinless 777We remark here that even one of the species would possesses spin degrees of freedom they are considered to be frozen as is the case of a spin polarized ensemble.. Each system is initialized in the ground state of the corresponding species with an initial finite spatial offset $x_{0}^{A}=-x_{0}^{B}$. To induce the dynamics, we quench this spatial offset to zero, i.e. $x_{0}^{A}=-x_{0}^{B}=0$, letting the initially separated atomic clouds collide.

III.2 Few-atom Bose-Fermi mixture: Convergence study

In the present section, we consider the dynamics of a few-body BF mixture and show in detail the convergence behavior of the ML-MCTDHX method. In particular, the mixture consists of $N_{B}=2$ bosons and $N_{F}=2$ fermions with weak bosonic intra-species repulsive interaction $g_{BB}=0.05$ and strong repulsive inter-species interaction $g_{BF}=1.0$, while the initial spatial offset is $x_{0}^{F}=-x_{0}^{B}=2$.

To study the dynamics we employ the evolution of the reduced one-body density for the ${\sigma}$-th species

[TABLE]

where $\hat{\Psi}_{\sigma}^{\dagger}(x)$ ($\hat{\Psi}_{\sigma}(x)$) refers to the field operator that creates (annihilates) a $\sigma$ species particle at position $x$. Figs. 1 ($a$)-($d$) present the dynamics from a single particle perspective via the one-body density for both the fermionic and the bosonic species. As shown $\rho^{\sigma}_{1}(x;t)$ exhibits oscillatory patterns possessing the frequency $\omega=1$ of the harmonic trap such that at regular time intervals, $\Delta t=\frac{(2n+1)T}{4}$ ($T=2\pi$ being a corresponding period) the two species collide. In the mean-field case, see Figs. 1 ($a$), ($b$), the only effect observed is the accumulation of the single particle density before a collision event for both the fermionic and bosonic species. However in the correlated case, see Figs. 1 ($c$), ($d$), the out-of-equilibrium dynamics differs significantly. As it can be seen, the density distribution $\rho^{F}_{1}(x;t)$ (as well as $\rho^{B}_{1}(x;t)$) oscillates as a whole up to the first collision event (see $t=\frac{1}{2}T$) and then splits into counter-propagating fragments. In addition, during each collision event the density shows a characteristic peak, which is more pronounced than in the mean-field case. The observed splitting process can be seen as a partial reflection of the particles during a collision. Furthermore note that the reflected part of the one-body density lies within the spatial region for which the density of the other species is dominant and inherits its density structure.

Since beyond mean-field effects apparently play a role already on the level of the reduced one-body density, we first study the convergence behavior of $\rho^{\sigma}_{1}(x;t)$ upon increasing the basis size in our variational wave-function ansatz [see Eqs. (5), (9)]. Figs. 1 ($e$)-($j$) show the profiles of the one-body density for both the fermions and the bosons using different $M-(m_{F},m_{B})$ approximations at various time instants during the evolution. For increasing number of SBSs $M$ (keeping the number of SPFs $m_{B}$, $m_{F}$ fixed) the corresponding profiles are almost indistinguishable (see solid red and light blue lines respectively). For instance, the corresponding one-body densities obtained e.g. within the approximations 3-(8,8) and 8-(8,8) respectively, deviate at most by $10\%$. However, the difference between the 7-(8,8) and 8-(8,8) approximations becomes negligible, showing a maximal deviation of the order of $1.2\%$ (results not shown here). The same behavior can be observed for increasing number of SPFs, namely by increasing $m_{B}$ and/or $m_{F}$ for a fixed number of NSFs $M$ (see light blue and black dashed lines). Here the maximum deviation between the approximations 8-(8,8) and 8-(10,10) is of the order of $2\%$. However, the mean-field approximation 1-(2,1) differs significantly from the correlated approach 8-(10,10) possessing completely different density profiles (see solid red and dashed black lines).

In order to quantitatively judge the convergence behavior of $\rho^{\sigma}_{1}(x;t)$ over the whole evolution time, we inspect the spatially integrated density difference

[TABLE]

where $N_{\sigma}$ refers to the total number of particles for the species $\sigma$, $\rho^{\sigma}_{1,C}(x;t)$ ($\rho^{\sigma}_{1,C^{\prime}}(x;t)$) denotes the one-body density of species $\sigma$ calculated within the configuration $C=M-(m_{F},m_{B})$ [$C^{\prime}=M^{\prime}-(m^{{}^{\prime}F},m^{{}^{\prime}B})$]. To judge about convergence on the single-particle level we calculate $\Delta_{CC^{\prime}}^{\sigma}(t)$ by successively increasing $M$ and fixed $m_{F}$, $m_{B}$ or with successively increasing the $m_{B}$ or $m_{F}$ keeping $M$ fixed.

Figs. 1 ($k$), ($l$) show in a transparent way an adequate convergence of our results both for fermions [see Fig. 1 ($k$)] and bosons [see Fig. 1 ($l$)] for increasing number of SBSs and in particular for $M>7$. Indeed, $\Delta_{CC^{\prime}}^{\sigma}(t)$ testifies large deviations for small number of SBSs (e.g. incrementing from $M=1$ to $M=2$), while for increasing $M$ (here $M=7$ to $M=8$) $max[\Delta_{CC^{\prime}}^{\sigma}(t)]\approx 1\%$. The latter indicates convergence of the reduced one-body densities. The same observations hold for increasing number of SPFs in the fermionic and/or the bosonic component of the binary mixture. As shown in Figs. 1 ($m$), ($n$) the one-body densities for both the fermions and the bosons respectively are more sensitive to the addition of more orbitals $m_{F}$ in the fermionic component. This is an expected result as spinless fermions, in contrast to bosons, according to the Pauli principle cannot tend to condense to a particular state. An adequate convergence for increasing number of fermionic orbitals between the configurations $C=8-(8,10)$, $C^{\prime}=8-(9,10)$ (see solid red line) and $C=8-(9,10)$, $C^{\prime}=8-(10,10)$ (see dashed orange line) is observed. In particular, $\max[\Delta_{CC^{\prime}}^{\sigma}(t)]\approx 0.5\%$ (see solid light blue and blue dashed lines) for $m_{F},m_{B}>9$. Summarizing, we observe that the accuracy of the calculation is more sensitive upon increasing the number of SBSs than incrementing the number of SPFs. Also, to achieve convergence we need to add more SPFs for the fermions than for the bosons which can be explained by the fact that bosons are allowed to condense in the same state but fermions not.

Having investigated the convergence behavior of local one-body observables, we now turn to more involved quantities, namely the populations of NOs and NSFs, which reflect certain aspects of the structure of the many-body state (see Sec. II D). The occupations of the NSFs and NOs are important for both the interpretation of the dynamics as well as to check the convergence of the method. Fig. 2 ($a_{1}$) presents the evolution of the populations of the NSFs for the system examined in Fig. 1. As it can be seen, $\lambda_{1}(t=0)=1$, implying the absence of inter-species correlations initially, while as time evolves $\lambda_{1}(t)$ deviates significantly from unity and $\lambda_{i}(t)$, $i\neq 1$, acquire non-negligible populations. Indeed all $\lambda_{i}(t)$, $i=1,2,...$, change significantly during a collision event but stay approximately constant in between. Figs. 2 ($a_{2}$), ($a_{3}$) illustrate in a transparent way, the convergence regarding the first and last NSF population respectively for varying number of used SBSs, namely $M=6,7,8$. Concerning the occupation of the least contributing NSF for the case of $M=8$ we observe that it is always very small, i.e. less than $10^{-3}$. The latter indicates that the used number of SBSs is sufficient for the convergence of the observables of interest (see Sect. II.4). The occupations of the first NSF during the evolution for different orbital configurations are almost indistinguishable. To quantitatively demonstrate convergence we define the relative difference of the population of the $i$-th NSF

[TABLE]

where $\lambda_{i}^{C}(t)$ [$\lambda_{i}^{C^{\prime}}(t)$] denotes the population of the $i$-th NSF calculated within the $C=M-(m_{F},m_{B})$ [$C^{\prime}=M^{\prime}-(m^{{}^{\prime}F},m^{{}^{\prime}B})$] orbital configuration. Figs. 2 ($a_{4}$), ($a_{5}$) present the relative difference for the first $\delta\lambda_{1}^{CC^{\prime}}(t)$ and the second $\delta\lambda_{2}^{CC^{\prime}}(t)$ population of the NSFs. A decreasing relative difference for increasing $M$ is clearly visible and the relative deviations are very small anyway, i.e. less than $4\%$ when comparing the $C=7-(10,10)$ calculation with $C^{\prime}=8-(10,10)$, for example, which further indicates convergence.

Moreover, Fig. 2 ($b_{1}$) presents the natural populations (NPs) of the fermionic natural orbitals during the dynamics. As $N_{F}=2$, the fermionic ensemble initially occupies a HF-like state characterized by the occupation of the first two orbitals, i.e. $n_{1}^{F}(0)=n_{2}^{F}(0)\approx 0.5$. However, following the first collision, other fermionic orbitals acquire significant populations. To judge about the obtained convergence Figs. 2 ($b_{2}$), ($b_{3}$) show the evolution of the first two and of the last fermionic NPs respectively using different number of fermionic orbitals $m_{F}$. Concerning the dependence of $n_{1}^{F}(t)$, $n_{2}^{F}(t)$ on the numerical configuration, we observe for the presented orbital numbers that they are almost indistinguishable during the evolution, see Fig. 2 ($b_{2}$). The latter is also confirmed by their corresponding relative differences $\delta n_{1}^{F,CC^{\prime}}(t)$ and $\delta n_{2}^{F,CC^{\prime}}(t)$ [defined in the same manner as the relative difference of the population of the NSFs, see Eq.(19)], shown in Figs. 2 ($b_{4}$), ($b_{5}$) respectively. A decaying relative difference for increasing number of fermionic orbitals [see $\delta n_{1}^{F,CC^{\prime}}(t)$ for $C=8-(8,10)$ and $C^{\prime}=8-(9,10)$ and $\delta n_{1}^{F,CC^{\prime}}(t)$ between $C=8-(9,10)$ and $C^{\prime}=8-(10,10)$] is observed with a corresponding deviation smaller than $0.6\%$, thereby indicating convergence. Finally, the smallest natural population, see Fig. 2 ($b_{3}$), is less than $10^{-3}$, indicating that there is no necessity of adding more orbitals to the fermionic component. Next we turn our attention to the evolution of the bosonic NPs, shown in Fig. 2 ($c_{1}$). The same overall behavior as for the fermionic NPs is observed. To illustrate convergence we explicitly show in Figs. 2 ($c_{2}$), ($c_{3}$) the evolution of the first and the last NP for varying number of bosonic NOs $m_{B}$. The time-evolution of the first NP, see Fig. 2 ($c_{2}$), remains almost invariant upon increasing the orbital numbers, while the last one, see Fig. 2 ($c_{3}$), is always below $10^{-3}$. Finally, the corresponding relative differences of the first [$\delta n_{1}^{B,CC^{\prime}}(t)$] and the second [$\delta n_{2}^{B,CC^{\prime}}(t)$] bosonic NPs using different number of $m_{B}$, are presented in Figs. 2 ($c_{4}$), ($c_{5}$). As expected, both $\delta n_{1}^{B,CC^{\prime}}(t)$ and $\delta n_{2}^{B,CC^{\prime}}(t)$ decrease for increasing number of bosonic NOs $m_{B}$ (see purple and green lines).

Summarizing, we have studied the convergence behavior of the reduced one-body density and the NOs and NSF populations, which are sensitive to both intra- and inter-species correlations. For the considered configurations, we have in particular seen that the accuracy is much more affected by the truncation on the species layer than on the particle layer. We have made similar observations in the other collision scenarios, too.

III.3 Particle-number imbalanced Bose-Fermi mixture

In the following, we investigate a complex collision scenario where again, for strong enough inter-species interactions, correlations qualitatively alter the transmission properties of the colliding species. Here, our main aim is to show that the NSFs, which can be easily extracted from the ML-MCTDHX wave-function ansatz, allow for linking the structure of the correlated many-body wave function to beyond mean-field effects in the reduced one-body density. Specifically, we consider a BF mixture consisting of $N_{B}=10$ bosons and $N_{F}=3$ fermions, while the initial positions of the species-dependent harmonic trap are $x_{0}^{F}=-x_{0}^{B}=2$. The bosonic intra-species interactions is kept always weak ($g_{BB}=0.05$) and the cases of weak and stronger inter-species interactions, namely $g_{BF}=0.5$ and $g_{BF}=2.0$ respectively, are examined.

For weak inter-species interactions, an inspection of the dynamics of the one body density reveals qualitatively similar characteristics between the mean-field, [see Figs. 3 ($a$), ($b$)] and the correlated [see Figs. 3 ($c$), ($d$)] description. In both cases the two initially separated clouds collide after the quench and perform oscillations with the frequency of the trap. The collisional events show in the one-body density evolution density peaks. These density peaks appear to be more pronounced in the mean-field case. Additionally, as it can be seen after the first collision a small fragment of the fermionic cloud is swept away from the majority cloud as a consequence of the collision with the bosonic ensemble. For times following the first collision, the bosonic density lies in between the above-mentioned majority fermionic cloud and the small reflected fermionic fragment. Also, the sweeping of the minority fermionic fragment is intensified in the course of the dynamics especially following each collision event. The latter is stronger in the correlated dynamics as compared to the mean-field behaviour.

Figs. 3 ($e$)-($h$) present the dynamics for stronger interactions, namely $g_{BF}=2.0$. In this case the evolution of the system within the mean-field approach is qualitatively different from the full many-body calculation. Indeed, in the mean-field case, the two species are mainly reflected, see Figs. 3 ($e$), ($f$). Only a tiny portion of the fermionic cloud transmits through the bosonic ensemble after each collision [not visible in Fig. 3 ($e$)]. Taking correlations into accout, however, one part of the fermionic density goes through the bosonic species after the first and even more after the second collision, see Figs. 3($g$), ($h$), and acquires more population during the dynamics.

In order to trace back the transmitted fermionic density to beyond HF behavior, one can evaluate the NOs densities and their populations. Being strongly averaged quantities however, the NOs do not allow for distinguishing intra- and inter-species correlation effects. As we shall see here, the latter can be accomplished for a binary mixture by the Schmidt decomposition Horodecki et al. (2009); Krönke et al. (2015)

[TABLE]

with $\{|\Psi_{k}^{\sigma}(t)\rangle,\;k=1,...,M\}$ being orthonormal. This decomposition can be easily obtained by diagonalizing the single-species reduced density matrices $[\eta_{1,\sigma}]^{i}_{l}$, which are typically compact objects in ML-MCTDHX calculations. Then, the one-body reduced density $\rho^{\sigma}_{1}(x;t)$ of the $\sigma$ species can be decomposed as

[TABLE]

where $\rho^{\sigma}_{1,k}(x;t)=\braket{\Psi^{\sigma}_{k}(t)}{\hat{\Psi}_{\sigma}^{\dagger}(x)\hat{\Psi}_{\sigma}(x)}{\Psi^{\sigma}_{k}(t)}$ denotes the one-body density matrix for the $k$-th NSF. Indeed, the decomposition of the one-body density into $\lambda_{k}(t)\rho^{\sigma}_{1,k}(x;t)$ modes [see Eq. (21)] sheds light on the correlated dynamics and the corresponding differences to the mean-field case. The leading contribution $\lambda_{1}(t)\rho^{F}_{1,1}(x;t)$, see Figs. 4 ($a$), ($b$), is similar to the mean-field one-body density $\rho^{F}_{1,{\rm MF}}(x;t)$ as the fermionic density is shown to be completely reflected after each of the collisions. On the contrary, $\lambda_{2}(t)\rho^{F}_{1,2}(x;t)$ clearly shows that a small fragment of the fermions is transmitted through the bosonic cloud and the remaining fermionic cloud is reflected, see Figs. 4 ($c$), ($d$). Furthermore, in $\sum_{k>2}\lambda_{k}(t)\rho^{F}_{1,k}(x;t)$ we also observe that a fermionic fragment passes through the bosonic cloud, see Figs. 4 ($e$), ($f$). For the bosonic component, we do not find striking structural differences between the NSF densities but only the tendency that less populated NSFs are more smeared out in space, cf. $\sum_{k>2}\lambda_{k}(t)\rho^{B}_{1,k}(x;t)$ with $\lambda_{1}(t)\rho^{B}_{1,1}(x;t)$ and $\lambda_{2}(t)\rho^{B}_{1,2}(x;t)$.

Concluding, we have linked the transmitted fermionic density to inter-species correlations by diagonalizing the density operator of the fermionic species. Because of inter-species interactions, the state of the $N_{F}$ fermions becomes mixed. Its dominant NSF density lacks a transmitted fraction and resembles the behavior of the corresponding mean-field calculation, whereas the sub-dominant NSFs feature a transmitted fraction.

III.4 Fermi-Fermi mixture

Finally, we study the collision dynamics of a FF mixture for a two-fold purpose. We can relate - this time - a reflected fraction of the density to inter-species correlations by applying the already introduced NSF analysis. Furthermore, a NO analysis reveals that the occupied single particle states can be classified into certain categories, manifesting the emerging NO structure of the beyond HF dynamics.

To illustrate the emergence of inter- and intra-species correlations in fermionic ensembles we consider a binary mixture consisting of $N_{A}=N_{B}=6$ fermions and $x_{0}^{A}=-x_{0}^{B}=4.5$ in order to avoid spatial overlap between the species initially. We emphasize that the considered inter-species interactions are relatively weak, namely $g_{AB}=0.4$, but still give rise to intriguing correlation effects.

The effect of the inter-species interaction is manifested after a collision by a small portion of the fermionic density which oscillates in the opposite direction with respect to the majority of the fermionic particles, see Fig. 5 ($b$). This is in contrast to the mean-field HF case, see Fig. 5 ($a$), where no such effect is observed and only coherent oscillations of the fermions are present. To shed light onto the collisional fermionic dynamics we employ the corresponding Schmidt decomposition. The density evolution of the dominant NSF [see Fig. 5 ($c$)], shows that all the fermions oscillate coherently in the same direction, in agreement to the HF case. The second dominant NSF [see the density evolution in Fig. 5 ($d$)] consists of the majority of the fermionic density oscillating in the one direction and a small fragment in the opposite direction. Higher order NSF contributions correspond to states with a somewhat more extended density oscillating in the opposite direction [see Fig. 5 ($e$)].

The populations of the above-mentioned modes, presented in Fig. 5 ($f$), change over time in particular after each collision event (taking place at $t_{n}=\frac{(2n+1)\pi}{2}$, $n\in\mathbb{N}$) and exhibit an oscillatory stepwise pattern. To compare the depletion rate for the first NSF, i.e. the mode where all the particles oscillate around the center of the trap in the same direction, we show its occupation for different inter-species interactions $g_{AB}$, see Fig. 5 ($g$). It is manifest that stronger interactions lead to a faster built-up of inter-species correlations. Since reflected density fractions are exclusive signatures of higher order NSFs [compare Figs. 5 ($b$), ($d$), ($f$)], measuring the density profile can be used to probe inter-species correlations. Moreover, the interaction dependence observed in Fig. 5 ($g$) can be explained by the fact that the reflection rate for the scattering between each pair of particles gets higher for increasing interaction and therefore the probability of all the particles to be transmitted decreases. Another interesting observation is that the depletion ($1-\lambda_{i}$) reaches a maximum value and subsequently reduces, e.g. see Figs. 5 ($d$), ($g$). This is due to the collision process between the two species which also allows for population transfer from the second NSF to the first.

Similarly to the bosonic case, where the populations of the NOs quantify the degree of fragmentation Penrose and Onsager (1956); Mueller et al. (2006), NOs and their populations provide also valuable insights into the structure of the many-body state of interacting fermions. Due to the fermionic statistics an amount of the initially occupied orbitals are almost unperturbed (core orbitals) during the dynamics, while the remaining orbitals (valence orbitals) exhibit depletion towards higher lying excited orbitals not populated in the initial state. In the present case, see Fig. 5 ($h$), we identify (from the corresponding populations) the first three $\Phi_{1}$ to $\Phi_{3}$ orbitals as the core orbitals, the orbitals $\Phi_{4}$ to $\Phi_{6}$ as the valence orbitals, and finally the orbitals $\Phi_{7}$ to $\Phi_{9}$ as the significantly occupied excited orbitals. The core orbitals are characterized by low depletion ($1/N_{\sigma}-n_{i}^{\sigma}$), as shown in the inset of Fig. 5 ($h$), while the depletion of the valence orbitals is related to the occupation of the excited orbitals. Figs. 6 ($a$)-($c$) present $\left|\Phi_{i}(x,t_{1})\right|^{2}$ with $i\in\{1,\dots,9\}$ at $t_{1}=9.8\approx\frac{3\pi}{2}$ where the two species have maximal separation. The core orbitals possess a similar appearance as the three lowest-energy eigenstates of the harmonic oscillator [visible in the nodal structure in Fig. 6 ($a$)]. On the contrary, the three valence orbitals [see Fig. 6 ($b$)] of each species are closer to each other than the core orbitals. Finally, the excited orbitals [see Fig. 6 ($c$)] of one species are located in the spatial region where the density of the other species is dominant. Hence, the valence orbitals of a given species overlap with the excited orbitals of the other species. The latter reveals that the valence orbitals become more affected than the core orbitals by the inter-species interaction. It is important to remark that the same kind of orbital structure persists during the evolution, see for instance Figs. 6 ($d$)-($f$) where $\left|\Phi_{i}(x,t_{1})\right|^{2}$ with $i\in\{1,\dots,9\}$ at $t_{2}=46.9$ is presented. Note that the latter is again a point of maximal separation. Our analysis for the NOs implies the following intuitive picture for the beyond HF fermionic orbitals. The energetically lowest initially occupied orbitals cannot be excited to unoccupied excited states by interactions with another fermionic species for low interaction strengths, in contrast to the energetically higher initially occupied orbitals being more susceptible to such excitations.

In conclusion, we have studied a binary FF ensemble and showed that beyond HF effects can be observed even in the weakly interacting regime. The behavior of the ensemble on the level of NSFs is similar to the aforementioned BF case where the first NSF resembles the mean-field dynamics while the higher order contributions give rise to the observed beyond mean-field effects. Additionally, the structure of the NOs reveals that only higher-lying initially occupied fermionic states can be excited due to the presence of inter-species interactions.

IV Conclusions

In this work, we present the ML-MCTDHX wave-function propagation method, which is designed for efficient ab-initio simulations of arbitrary multi-component ultracold atomic ensembles out of equilibrium. Employing both the concept of a variationally optimally moving many-body basis and a multi-layer ansatz for the wave-function expansion allows for adapting ML-MCTDHX to system-specific intra- and inter-species correlations, which greatly improves its efficiency. By virtue of the so-called reduced transition matrices, we can formulate the ML-MCTDHX EOM in a unified manner, treating bosonic and fermionic components formally on an equal footing. Moreover, using the reduced transition matrices as central building blocks for evaluating the EOM proves to be computationally beneficial.

We demonstrate the power of the ML-MCTDHX method by applying it to the collision dynamics of few-atom binary mixtures of Bose-Fermi and Fermi-Fermi type in one dimension. Beginning with a detailed analysis of the convergence behavior of the densities, natural orbitals, and natural species functions by successively enlarging the employed basis sets, we show how the intra- and inter-species correlations are described with increasing accuracy. The truncation on the species layer is found to affect dominantly the accuracy of the computation, while, however, reducing the required number of coefficients for a converged simulation significantly: For e.g. the Fermi-Fermi collision scenario, a wave-function expansion w.r.t. to all numberstates $|{\bm{n_{A}}}\rangle^{A}_{t}|{\bm{n_{B}}}\rangle^{B}_{t}$ would involve $K_{A}K_{B}\approx 25\cdot 10^{6}$ configurations, whereas the species-layer truncation reduces the corresponding number of coefficients by a factor of about $500$ to $M^{2}+M(K_{A}+K_{B})\approx 5\cdot 10^{4}$. Thereafter, we reveal the impact of correlations on the collision dynamics induced by the interspecies interaction. For a colliding BF mixture, we in particular observe a correlation-induced transmitted fraction, which is highly suppressed in corresponding mean-field calculations. This process takes place only in the sub-dominant natural species functions leading to a mixed state for the individual species. For a colliding FF mixture, we similarly find that, even at relatively weak interactions, correlations lead to a partial reflection via the sub-dominant natural species function. On the single-particle level, we classify the natural orbitals as core, valence and excited orbitals, where the population of the latter is a further signature of the beyond HF physics observed.

In the future, ML-MCTDHX will give access to the quantum dynamics of a multitude of physical mixture systems such that correlation-induced processes can be analyzed in detail. Here, we have in particular induced interactions, dephasing, self-localization, entangled state generation, as well as inter-species energy and momentum transfer in mind with a particular focus on the question how the distinguishability and particle statistics affects the physical behavior. Finally, we would like to stress that in order to unravel the correlated processes, further development of analysis tools for many-body systems needs to be put forward such as the here employed natural-species-function analysis shown to separate mean-field and correlated behavior.

Acknowledgements.

The authors acknowledge inspiring discussions with H.-D. Meyer. This work has been financially supported by the excellence cluster ’The Hamburg Centre for Ultrafast Imaging - Structure, Dynamics and Control of Matter at the Atomic Scale’ of the Deutsche Forschungsgemeinschaft. Financial support by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the SFB 925 ”Light induced dynamics and control of correlated quantum systems” is gratefully acknowledged by L.C., S.I.M. and P.S.

Appendix A Reduced few-body transition matrices

The reduced one and two-body density matrices, the mean-field operators and the Hamiltonian matrix are key ingredients in the ML-MCTDHX EOM (11), (12) and (13). The calculation of these quantities shares a common step, which is to calculate the reduced one- and two-body transition matrices:

[TABLE]

Here, $\hat{k}$ refers to an occupation number vector with zero occupations except for the occupation of the $k$-th SPF being set to unity. The particle-exchange symmetry enters the calculation only via the occupation number dependent factors $Q_{{\bm{n}}}(k,q)$ and $P_{{\bm{n}}}(k,q)$, which read for bosonic species

[TABLE]

and for fermionic species

[TABLE]

where we have defined

[TABLE]

and $\theta(k,q)=1$ if $k>q$ and zero otherwise.

Appendix B Hamiltonian matrix

In order to calculate the Hamiltonian matrix $\langle\psi_{i_{A}}^{A}\psi_{i_{B}}^{B}|\hat{H}|\psi_{i_{A}^{\prime}}^{A}\psi_{i^{\prime}_{B}}^{B}\rangle$ in SBS representation, as needed for the top-layer EOM (11), one may use the second-quantization representation of $\hat{H}$. Bearing in mind that the SBSs consist only of number states in which SPFs are occupied 888But no other states of the full single-particle Hilbert space $\mathfrak{h}^{\sigma}$., only the creation and annihilation operators of SPFs provide non-vanishing contributions. Consequently, we find $\langle\psi_{i_{A}}^{A}\psi_{i_{B}}^{B}|\hat{H}|\psi_{i_{A}^{\prime}}^{A}\psi_{i^{\prime}_{B}}^{B}\rangle$ to be the sum of the following terms

[TABLE]

with $[h_{\sigma}]_{s}^{r}=\langle\phi_{r}^{\sigma}|\hat{h}_{\sigma}|\phi_{s}^{\sigma}\rangle$ as well as $[v_{\sigma}]_{uv}^{rs}=\langle\phi_{r}^{\sigma}\phi_{s}^{\sigma}|\hat{v}_{\sigma}|\phi_{u}^{\sigma}\phi_{v}^{\sigma}\rangle$, and

[TABLE]

with $[w_{AB}]_{sv}^{ru}=\langle\phi_{r}^{A}\phi_{u}^{B}|\hat{w}_{AB}|\phi_{s}^{A}\phi_{v}^{B}\rangle$ and the transition matrices given in Appendix A.

Appendix C Mean-field operator matrices

The mean-field operator matrices on the species-layer $[\hat{W}_{\sigma|\bar{\sigma}}]^{i}_{j}=\langle\psi_{i}^{\bar{\sigma}}|\hat{W}_{AB}|\psi_{j}^{\bar{\sigma}}\rangle$ [see Eq. (12)] can be expressed in terms of the reduced one-body transition matrix as well

[TABLE]

On the particle layer, the corresponding intra- and inter-species mean-field operator matrices read

[TABLE]

Appendix D Density matrices

The reduced density matrices, entering the EOM for the time-dependent basis states (12) and (13), can be calculated as follows. The reduced density matrix of the $\sigma$ species in SBS representation is given by

[TABLE]

where the sum over the multi-index $J^{\sigma}$ denotes tracing out all the species but $\sigma$ and the multi-index $J^{\sigma}_{i}$ refers to the corresponding configuration of all species with the index for the $\sigma$ species being set to $i$. For a binary mixture, the reduced two-species density matrix reads [\eta_{2,AB}]^{i_{A}i_{B}}_{j_{A}j_{B}}=[\eta_{2,BA}]^{i_{B}i_{A}}_{j_{B}j_{A}}={\big{(}A_{i_{A}i_{B}}\big{)}^{*}A_{j_{A}j_{B}}}, while for $S>2$ species we have

[TABLE]

with the sum running over all SBSs of species other than $\sigma$, $\kappa$ and $J^{\sigma\;\kappa}_{i_{\sigma}i_{\kappa}}$ denoting the corresponding multi-index with the $\sigma$ ($\kappa$) species being in the $i_{\sigma}$-th ($i_{\kappa}$-th) SBS.

The reduced one-body density matrix of the $\sigma$ species can now be expressed as

[TABLE]

Correspondingly, we find for the reduced two-body intra-species density matrix

[TABLE]

Finally, the reduced two-body inter-species density matrix can be calculated as follows

[TABLE]