Microscopic scattering theory for interacting bosons in weak random potentials

We develop a diagrammatic scattering theory for interacting bosons in a three-dimensional, weakly disordered potential. Based on a microscopic N-body scattering theory, we identify the relevant diagrams including elastic and inelastic collision processes that are sufficient to describe diffusive quantum transport. By taking advantage of the statistical properties of the weak disorder potential, we demonstrate how the N-body dynamics can be reduced to a nonlinear integral equation of Boltzmann type for the single-particle diffusive flux. Our theory reduces to the Gross-Pitaevskii mean field description in the limit where only elastic collisions are taken into account. However, even at weak interaction strength, inelastic collisions lead to energy redistribution between the bosons - initially prepared all at the same single-particle energy - and thereby induce thermalization of the single-particle current. In addition, we include also weak localization effects and determine the coherent corrections to the incoherent transport in terms of the coherent backscattering signal. We find that inelastic collisions lead to an enhancement of the backscattered cone in a narrow spectral window for increasing interaction strength.


Introduction
In recent years, increasing interest has been devoted to the behaviour of ultracold atoms in disordered potentials. Whereas the first experiments [1,2,3] concentrated on the realization of Anderson localization [4] in one dimension, this intriguing disorder effect -which leads to complete suppression of diffusive transport due to destructive interference -has now also been observed in three dimensions [5,6]. The 3D case is especially interesting since it exhibits a transition from extended to localized singleparticle eigenstates: In the absence of interactions, particles with low energy are localized, whereas those with higher energy (in comparison with the strength of the disorder potential) propagate diffusively in the random potential. Also in the latter case -on which we concentrate in the present paper -wave interference effects are relevant, though less pronounced: They lead to weak localization [7] (i.e. reduction of the diffusion constant instead of complete suppression of diffusion) and, associated with that, coherent backscattering [8,9,10] (i.e. enhancement of backscattering), which has recently been observed also with atomic matter waves [11,12,13]. Beyond the scope of [11,12] lies the investigation of the interplay between disorder and interactions, where it is not well understood, especially in the higher-dimensional case, to what extent interaction leads to a loss of coherence, i.e. to a breakdown of localization effects [14,15]. Most theoretical works, e.g. [16,17,18,19,20], focus on the regime of -or close tothermal equilibrium and examine, e.g., the effect of disorder on the condensate fraction, superfluid fraction or the sound velocity [17,18,20]. In this case, weak interactions can usually be treated perturbatively, e.g., by introducing Bogoliubov quasiparticles [21].
In contrast, the present paper investigates a stationary scattering setup far from thermal equilibrium. Here, bosonic atoms are continuously emitted from a coherent source ('atom laser' [22,23]) and guided into the random potential until a stationary scattering state is reached. Theoretical studies of this scattering scenario so far either neglect the interparticle interaction [24,25], treat it on the mean-field level [26,27,28,29], or apply a Hartree-Fock-Bogoliubov approach [30], which is appropriate in the case of a large condensate fraction. If all atoms enter the scattering region at fixed initial energy, the Gross-Pitaevski equation obtained within the mean-field approach predicts either a stationary regime with the same final energy for all scattered atoms, or a non-stationary, time-dependent behavior [26,30]. In contrast, according to the microscopic scattering theory developed in the present paper, atoms exchange energy with each other due to mutual collision events, leading to strong depletion of the condensate already for small interactions. As shown in [31], this finally leads to a stationary state with thermal Maxwell-Boltzmann distribution for those atoms which propagate deeply into the scattering region.
The present paper is devoted to a detailed presentation of the underlying bosonic many-particle scattering theory. Starting from the N -particle Hamiltonian, we derive a nonlinear transport equation for the average particle density. Since this transport equation amounts to a stationary version of the Boltzmann equation [32], our approach is, in this respect, comparable to previous works on quantum kinetic equations [33,34,35,36,37,38,39,40,41]. In contrast to these works, however, the additional presence of a disorder potential (apart from the atom-atom interactions) in our setup allows us to quantify the regime of validity of the transport equation in a more rigorous way. In the regime of weak disorder, the disorder average enables us to neglect correlations between atoms induced by collisions, which in turn is the basic assumption required for reducing a many-particle problem to an effective single-particle description. Moreover -and again in contrast to the above works -we go beyond the case of purely diffusive transport, and also incorporate quantum interference corrections leading to coherent backscattering into our theory.
Correspondingly, the paper is structured as follows: In Sec. 2, we set the stage by reviewing some important aspects of standard scattering theory for a single particle. The case of many interacting particles will be addressed in Sec. 3: Starting from the Hamiltonian including pairwise atom-atom interaction, we introduce a diagrammatic notation for the transition amplitudes of many particles, from which the scattered flux density can be calculated after taking the trace over the undetected particles. As we will see, this trace leads to a distinction of atom-atom collisions events into inelastic and elastic collisions, respectively, where the latter are shown to reproduce the mean-field description given by the Gross-Pitaevskii equation. Whereas the methods presented up to Sec. 3 are generally valid for an arbitrary scattering potential, we focus on the case of a weak random potential from Sec. 4 on. The assumption of weak disorder (k dis 1 with wavenumber k and disorder scattering mean free path dis ) is crucial, since it allows to reduce the -in principle infinitely complicated [42] -hierarchy of many-particle diagrams to a tractable subclass of diagrams, i.e. ladder and crossed diagrams [43], which are composed out of a small number of building blocks. As shown in Sec. 4, the sum of all ladder diagrams amounts to a Boltzmann-like equation for diffusive transport eventually leading to complete thermalization due to inelastic atomatom collisions in case of an infinitely large scattering region. Sec. 5 is devoted to the derivation of transport equations describing coherent backscattering based on crossed diagrams. Finally, in Sec. 6 we present the results of numerical solutions of the ladder and crossed transport equations exemplifying the behaviour of diffusive transport for a finitely large scattering region, and the effect of elastic and inelastic atom-atom collisions on coherent backscattering, respectively. Sec. 7 concludes the paper. Several technical aspects are relegated to Appendices A-E.

Scattering theory for a single particle
We write the Hamiltonian for a single particle in the following form: whereĤ 0 denotes free propagation andV the disorder potential. The eigenstates |k of H 0 are plane waves with wave vector k: and energy where we seth 2 /(2m) ≡ 1. The matrix elements ofV are given by the Fourier transform of the disorder potential V (r): In order to obtain a properly defined scattering scenario, we assume that V (r) is nonzero only inside a finite scattering region V. This allows us to define an asymptotically free initial state: with normalized wavepacket w(k), i.e. dk|w(k)| 2 = (2π) 3 , which we assume to be a quasi-monochromatic wavepacket, i.e., sharply peaked around the initial wavevector . Therefore, the spatial density resulting from the Fourier transform of w(k): is approximately constant inside the scattering region, i.e. for r ∈ V. If the state exp(−iĤ 0 T )|i 1 is prepared at time T → −∞, the wavepacket arrives at the scattering region at time t = 0, and a quasi-stationary scattering state is reached at that time. Here, the operatorΩ whereĜ with infinitesimally small > 0, denotes the (retarded) Green's operator associated to the HamiltonianĤ =Ĥ 0 +V . The operatorsΩ whereĜ 0 (E) denotes the vacuum Green's operator: The operatorΩ is closely related to the Møller operatorΩ Since, in the following, we will applyΩ (V ) + (E) only to such eigenstates -or quasi-eigenstates, as |i 1 in Eq. (7) -we will henceforth refer also toΩ (V ) + (E) as 'Møller operator'. Finally, the expectation value of an arbitrary observableÂ in the (quasi-)stationary scattering state results as Â = f +,1 |Â|f +,1 .
Let us note that, instead of using the Møller operator, a scattering process can also be characterized by the S-matrix,Ŝ = Ω (V ) − is defined in the same way asΩ (V ) + , but with T → +∞ instead of −∞). We could formulate the following N -particle scattering theory equally well in terms of the S-matrix. However, since the S-matrix maps incoming onto outgoing asymptotically free states, it does not allowin contrast to the Møller operator -to evaluate what is happening inside the scattering region, e.g. to calculate the (quasi-)stationary density or flux of particles inside V. For this reason, we prefer using the (quasi-)stationary scattering state |f +,1 , see Eq. (7) (or its N -particle counterpart |f + , see Eq. (22) below) in the following.
In contrast toĤ 0 andV , the interactionÛ acts on two particles: with atom-atom interaction potential U (r). In the following, a collision event between two particles will be described by the T -matrix [44]: According to Eq. (17), the matrix elements ofT U (E) with respect to two-particle states describe repeated application of the interactionÛ on the same pair of particles, interrupted by free propagationĜ 0 (E). Separating the center-of-mass from the relative coordinates, the two-body T matrix fulfills momentum conservation: is the T -matrix for a single particle (with reduced mass m/2) scattered by the potential U (r) at energy The single-particle T -matrix, in turn, fulfills the optical theorem [44]: expressing conservation of the particle and the energy flux (whereĜ 0,m/2 denotes the vacuum Green's operator for a particle with mass m/2 and corresponding dispersion relation E = 2k 2 ). Our many-particle scattering theory presented below, and in particular the transport equations in Secs. 4 and 5, are valid for an arbitrary interaction potential U (r) -as long as it is sufficiently weak in the sense specified below (mean distance between collision events larger than between disorder scattering events). Only for the numerical results presented in Sec. 6, we will assume a short-range potential with corresponding s-wave scattering approximation, see Eq. (E.1).
Finally, we note that, in principle, the vacuum T -matrix as defined in Eq. (17) is modified by the presence of the disorder potential. To take this into account, the vacuum Green's operatorĜ 0 (E) must be replaced by the disorder Green's operator G V (E), see Eq. (9), in Eq. (17). However, since the present paper assumes the case of a very weak disorder potential, we will neglect the disorder during each collision event in the following, and therefore use the vacuum T -matrix as introduced above. This approximation is valid if the range of the interaction potential U (r) is much smaller than the disorder mean free path dis introduced in Sec. 4.

Many-particle transition amplitudes
We now generalize the scattering scenario outlined in Sec. 2 to the case of many particles. For this purpose, we assume that, both, the disorder and the particle-particle interaction are non-zero only inside a finite region V (which, for simplicity, we assume to be the same forV andÛ ). Note that the introduction of a finite interaction region in principle breaks translational invariance, and therefore the δ-function expressing momentum conservation in Eq. (18) turns into an approximate δ-function. Since, however, we assume the size L of the scattering region V to be much larger than the disorder mean free path, i.e. L dis k −1 (see below), we can safely neglect the associated small width (∝ 1/L) of this δ-function, and still work with the T -matrix as given by Eq. (18).
Our initial state for N particles reads: where all N atoms are described by the same quasi-monochromatic single-atom wavepacket w(k) as given in Eq. (5). The factor 1/ √ N ! arises from the indistinguishability of bosonic particles. The corresponding density of particles reads: As mentioned above, this density is approximately uniform within the whole scattering region V for a wavepacket sharply peaked around the initial wavevector k i . Since, in this quasi-monochromatic limit, the density, Eq. (6), for N = 1 approaches zero (since the wave packet is spread over an increasingly large region of space), the number N of particles correspondingly must tend to infinity in order to obtain a finite density ρ 0 . The Møller operator, which yields the quasi-stationary N -particle scattering state is defined in the same way as above, see Eqs. (8,9) but withV +Û instead ofV . It therefore fulfills the Lippmann-Schwinger equation: which, using Eqs. (8,11), can be rewritten as: Iteration of Eq. (24) yields an expansion in powers ofÛ : (25) Remember that, according to Eq. (16), each operatorÛ annihilates and creates two particles. In contrast, the Green's operatorĜ V and the Møller operatorΩ (V ) + act on all N particles. However, since these operators describe non-interacting particles, they can be factorized into single-particle operators. As an example, we give here the factorization formulas for the case N = 2: and Figure 1. Example of a three-particle scattering process with initial state |k 1 , k 2 , k 3 and final state |k 4 , k 5 , k 6 . The three arrows associated with the initial state represent the Møller operatorΩ (V ) + (E i ) of the disorder potential, see Eq. (8), whereas the remaining arrows refer to the disorder Green's operatorĜ V , Eq. (9). Squares correspond to the two-body T -matrix of the particle-particle interaction, Eq. (18). The transition amplitude corresponding to this scattering process is given in Eq. (28).
is included in the T -matrix, see Eq. (17) (and the discussion at the end of Sec. 3.1). We hence replace two-particle matrix elements ofÛ by matrix elements ofT U (E) (with appropriately defined two-particle energy E, see below) in Eq. (25), and thereby obtain a sequence of collision events between different pairs of particles. An example of a three-particle scattering process is demonstrated in Fig. 1. As shown in Appendix A, this diagram gives rise to the following contribution to the transition amplitude: with E k 1 E k 2 E k 3 E i , according to our above assumption of a quasimonochromatic wavepacket.
In general, the rules for constructing an arbitrary N -particle scattering amplitude for a given diagram are as follows: (i) Apply the disorder Møller operatorΩ (8), to each initial single-particle state |k 1 , . . . , |k N . The energy associated to each initial particle is given by E i . (ii) Integrate over all intermediate particles (p 1 , . . . , p 8 in Fig. 1). (iii) Write down the corresponding two-body T -matrix element, see Eq. (18), for any collision between two particles. The energy argument ofT U is given by the sum of the two incoming single-particle energies. (iv) For eachT U (E), write down an integral ∞ −∞ dE /(−2πi) which determines the energy arguments of the Green's operatorsĜ V (E ) andĜ V (E − E ), see Eq. (27), for the two particles after the collision. (v) These two particles may then collide with other particles, and so on ... .
The total transition amplitude defining the stationary scattering state |f + , see Eq. (22), is then obtained by summing the contributions from all possible different diagrams. For example, in addition to the diagram shown in Fig. 1, eight more diagrams obtained by exchanging the initial and/or final wavevectors (k 1 , k 2 , k 3 ) and (k 4 , k 5 , k 6 ) also contribute to |f + .

Scattered flux
As the finally measured quantity, we determine the expectation value of the flux density operatorĴ with respect to the stationary scattering state |f + . SinceĴ(r) is a one-particle operator, this implies a partial trace of the density matrix |f + f + | over N −1 undetected particles: Placing the detector at position R in the far field of the scattering region (i.e. |R| |r| for r ∈ V), the scattered flux is finally expressed as a dimensionless quantity (the socalled 'bistatic coefficient' [45]): normalized with respect to the incident flux Aρ 0 √ E i , where A denotes the transverse area (with respect to the incident wave) of the scattering volume V, andk d = R/|R| is the direction of the detected particle's wavevector. The limit R → ∞ is to be taken after the quasi-stationary limit N → ∞, see the discussion after Eq. (21). Apart from the total flux density γ(k d ), we will also be interested in the spectral density γ E (k d ), i.e. the flux of particles scattered into directionk d with energy E, which is given by: . The factor 1/N ! in Eq. (30) arises from the indistinguishability of the bosonic particles. It turns out, however, that this factor -together with the factors 1/ √ N ! in Eq. (20) -is exactly counterbalanced once we sum the amplitudes of all processes where the initial and/or final particles are exchanged. In total, we get the same result as if the particles were distinguishable. This equivalence is generally valid if all particles are prepared in the same initial state, and if the Hamiltonian is symmetric under exchange of particles [46].
Remember that the number N of particles tends to infinity in the quasi-stationary limit, whereas, in case of a finite scattering region, only a finite number of particles will eventually interact with the finally detected particle. The evolution of the remaining . The half circle symbol denotes the detector, whereas the dots on the left-hand side of both equations represent the trace over the undetected particle. a) Inelastic scattering of two particles. On the right-hand side, the trace has been performed using Eq. (33). This results in the dashed-solid double arrow representing the spectral The energy of the detected particle is 2E i − E. b) Elastic scattering of two particles. The trace is performed using Eq. (34). The resulting diagram on the right-hand side is equivalent to a diagram obtained from the Gross-Pitaevski equation. Since, in contrast to a), the conjugate particles (dashed lines) do not undergo a collision, the energy of the detected particle is unchanged (E i ).
particles (which do not interact with the detected particle) does not influence the result of the partial trace, Eq. (30). This follows from the factorization property, Eq. (26), and the left-unitarity, Ω + †Ω + = 1 of the Møller operator. Consequently, in order to calculate the detection signal, we may disregard all scattering processes concerning those particles which do not interact (neither in |f + nor in f + |) with the detected particle. (The presence of these particles only leads to a prefactor giving rise to the correct dependence of a given scattering diagram on the density ρ 0 , see the discussion at the end of Appendix B.)

Trace over undetected particles
According to the recipe given above, the flux density for an arbitrary N -particle scattering process is obtained as follows: take a diagram contributing to |f + , a conjugate diagram contributing to f + |, apply the observableĴ(r) to one of the final particles of both diagrams, and trace over the undetected particles. An example for two particles is shown in Fig. 2a) (left-hand side). Since both conjugate diagrams (solid and dashed lines, respectively) exhibit a collision event, which redistributes the energy among the two particles according to the factorization formula, Eq. (27), the energy of the detected = + Figure 3. Elastic scattering diagram where the conjugate undetected amplitude originates from a previous elastic scattering event. The sum of the two processes shown on the left-hand side reproduces the Gross-Pitaevskii diagram (cf. [29]) on the right-hand side.
particle is different from the initial energy E i . For this reason, we call this scattering process 'inelastic'. This means that the energies of the single particles change -although their sum remains conserved. In contrast, Fig. 2b) shows an elastic scattering process. Here, the conjugate diagram (dashed lines) on the left-hand side does not exhibit a collision event. As shown below, this implies that the energies of both particles remain unchanged.
We will now demonstrate how to perform the trace over the undetected particle for inelastic and elastic collisions, respectively. The result is represented on the right-hand side of Fig. 2.
Inelastic collisions. The complete expression for the inelastic scattering diagram, Fig. 2a), is given in Eq. (B.1). Focusing on those terms which are relevant for the trace over the undetected particle, this trace can be written in the following general form: On the left-hand side of Eq. (33), k corresponds to the final state of the undetected particle, whereasĜ V (E) andĜ † V (E ) refer to the (single-particle) Green's operators expressing propagation from the collision event to the final state. According to the rules given in Sec. 3.2, the collision events are associated with integrals dE/(−2πi) and dE /(2πi) which determine the energies of the undetected (E and E ) and the detected particle (2E i − E and 2E i − E ). The brackets (. . .) Under these conditions -which do not only hold for the example shown in Fig. 2, but for all other inelastic scattering diagrams we will encounter in the following -the result of the trace is given on the right-hand side of Eq. (33). This general formula is proven in Appendix C. Graphically, the result is depicted on the right-hand side of Fig. 2a). As a consequence, the energies E and E are set equal to each other, and the two conjugate Green's functionsĜ † (which is also known as the 'spectral function', since the imaginary part of the Green's function determines the density of states [47]).
Elastic collisions. In a similar way, the trace in the elastic scattering diagram, see Eq. (B.2), is performed as follows: see Appendix C. According to Eq. (34) -which is graphically depicted in Fig. 2b)the outgoing solid arrow emitted from the two-body collision event is replaced by an incoming dashed arrow with energy E i . We note that precisely this diagram is the only interaction contribution generated by the Gross-Pitaevskii equation [29]. Thereby, we have shown that our N -particle scattering theory reproduces the Gross-Piatevskii equation if only elastic scattering is taken into account. In Eq. (34), the conjugate undetected particle originates directly from the initial state k i | propagated in the disorder potential through the Møller operator Ω (V ) The formula can be generalized, however, to the case where the undetected particle undergoes previous collisions with other particles before colliding with the detected particle. An example is depicted in Fig. 3. Also in this case, the corresponding Gross-Pitaevskii diagram is reproduced (i.e. the outgoing solid arrow is replaced by an incoming dashed arrow). In a similar way, also the inelastic trace formula, Eq. (33), is valid in the case where the undetected particle undergoes further collisions with other particles before the trace is taken -provided that none of these other particles, in turn, collides with the detected particle which, as discussed in Sec. 4, is the case for a weak disorder potential. This allows us to take the trace over the undetected particles directly after their last collision with the detected particle -without being obliged to follow their further evolution before finally leaving the scattering region.

Ladder diagrams
The N -particle scattering formalism outlined above is valid for an arbitrary potential V (r). Now, we consider V (r) as a random potential, and calculate the corresponding average density matrix |f + f + |. For this purpose, we assume a Gaussian white noise Figure 4. Example of a ladder diagram describing the propagation of three interacting particles in a slab with a random scattering potential. Pairs of conjugate amplitudes (solid and dashed arrows, respectively) undergo the same sequence of scattering events (encircled crosses) induced by the disorder potential, see Eq. (35), at r 1 , . . . , r 9 . Due to particle-particle collision events (squares), the particles redistribute their energies. Here, solid and dashed arrows correspond to disorder averaged single-particle Green's functions, Eq. (36), and their complex conjugates, respectively. Upon flux detection, one particle is annihilated, while the undetected particles are traced over (dots). potential, specified by the mean value V (r) = 0 and the two-point correlation function: Furthermore, the disorder potential is assumed to be weak, i.e. √ E dis 1 for all relevant single-particle energies E, see Eq. (3). Initially, this is the case if √ E i dis 1. Due to inelastic collisions, the energies will change, but, as we will see later, their distribution will still be centered close to E i , with only a negligible fraction of particles that reach single-particle energies E 0.
For the case of a single particle, the disorder average in the limit k dis 1 is well known [48,49]: first, the vacuum Green's functionĜ 0 , see Eq. (12), is replaced by the average single-particle Green's function: with wherek E = √ E + i/(2 dis ). In position representation, this leads to an exponential decay of the average density with 2 Imk E = 1/ dis as the decay constant, see Eq. (39) below. This establishes dis as the mean free path, i.e. the average distance between subsequent disorder scattering events. Second, when calculating the average density matrix |f + f + |, and representing both |f + and f + | as a sum of diagrams, only those combination of diagrams survive where both |f + and f + | undergo the same sequence of disorder scattering events. Here, a disorder scattering event is induced by the correlation function, Eq. Figure 5. Trace over the undetected particle for disorder-averaged diagrams (see correlators with V (r 1 ) and V (r 2 ) both acting in |f + or both in f + | are accounted for by the average Green's function (37) [48]). These combinations of diagrams give rise to so-called ladder diagrams for the average density [49]. We now apply the same procedure to the N -particle scattering processes presented in Sec. 3.2. First, we take a diagram contributing to |f + and another one (called 'conjugate diagram' in the following) contributing to f + |. Then, we replace all vacuum Green's functions by average Green's functions and correlate, using Eq. (35), each disorder scattering event with another one in the conjugate diagram such that both conjugate diagrams undergo the same sequence of disorder scattering events. Finally, we choose one of the final particles as detected particle, and trace over the remaining N − 1 particles, see Eq. (30). An example is shown in Fig. 4. In this figure, the trace over undetected particles is performed as soon as the corresponding particle (solid line) is re-united with its conjugate counterpart (dashed line) at a disorder scattering event. It turns out that the same result is obtained if the trace is performed before taking the disorder average according to Fig. 2. This leads to disorder-averaged trace formulas as depicted in Fig. 5, which we will use in the following to evaluate the trace over the undetected particles. Among the N -particle ladder diagrams thus constructed, we neglect all those where two particles which interacted once meet again. This approximation is equivalent to the neglect of recurrent scattering [50] for a single particle, which, alike the neglect of nonladder diagrams, is valid for k dis 1. It allows us to trace away the undetected particles directly after their interaction with the detected particle. Finally, we assume that at least one disorder scattering event occurs between two collision events. This is with σ denoting the scattering cross section of the atom-atom interaction potential U (r), defines the average distance between two inelastic collision events. For s-wave Figure 6. The three building blocks from which all ladder diagrams (see Fig. 4) are constructed. a) Single-particle propagation in the disorder potential, see Eq. (39). b) Elastic two-particle collision g E1,E2 (r 1 , r 2 , r 3 ), see Eq. (40). c) Inelastic two-particle collision f E1,E2,E3 (r 1 , r 2 , r 3 ), see Eq. (41).

Building blocks
The trace over the undetected particle allows us to decompose every ladder diagram (like the one shown in Fig. 4) into independent building blocks. These building blocks are shown in Fig. 6. The first one, Fig. 6a), represents a single average propagation step of a single particle with energy E and corresponding wave vector k in the disordered potential from r 1 to r 2 : for E ≥ 0. Note that, for a white noise potential as defined in Eq. (35), the mean free path dis is independent of E [53]. For E < 0, the propagation is exponentially suppressed; in this case, Eq. (39) is multiplied by an additional factor exp −2|r 1 − r 2 | |E| ). Since the typical distance between two scattering events is given by the mean free path dis , we can neglect the occurrence of negative energies if |E| dis 1. The second building block, Fig. 6b), represents an elastic collision event, where the energies of both particles are unchanged: The trace over the undetected particle was performed according to Fig. 5b), giving rise to an average Green's function G * E 1 (k 4 ). For reasons of clarity, the wave vectors k 1 , . . . , k 5 are not explicitly shown in Fig. 6. They can, however, be easily deduced from the phase factors exp(±ik · r) describing annihilation or creation of a particle k due to disorder scattering at r, see Eqs. (4,15), with the help of following rule: outgoing solid (dashed) arrows always contribute with negative (positive) sign, the opposite holds for incoming arrows. For example, k 5 -with phase factor exp[ik 5 · (r 2 − r 3 )] in Eq. (40)is associated with the dashed arrow pointing from r 2 to r 3 in Fig. 6b).
The first factor 2 in Eq. (40) originates from the fact that the solid and dashed incoming amplitudes can be grouped together in two different ways. It can be shown that this accounts for fluctuations of the atomic density inside the disordered slab [28]. The factor 1/2 in front of the integral originates from the indistinguishability of particles, see the discussion at the end of Appendix B.
Finally, the third building block, Fig. 6c), amounts to an inelastic collision event, where the energies of two particles E 1 and E 2 change to E 3 and where Fig. 5a) was used for the trace over the undetected particle.

Transport equation
The outgoing arrows of each building block may now be attached to the incoming arrows of the next building block, and so on. The sum of ladder diagrams resulting from all combinations of these building blocks is expressed by the following nonlinear integral equation: where represents the incoming wave propagating to r without being scattered. Correspondingly, z r,−k i denotes the distance from the surface of the scattering region V to r along a straight line parallel to the direction k i of the incident wavepacket. The quantity I E (r) can be interpreted as the average density of particles with energy E at position r (at least in the case dis √ E 1 of weak disorder where the spectral function [G * E (k) − G E (k)]/(2πi) → δ(E − k 2 ) approaches a δ-function, such that a particle with wavevector k possesses a well defined energy). In particular, I(r) = dEI E (r) = f + |ρ(r)|f + gives the disorder-averaged expectation value of the single-particle density operatorρ(r) =ψ † (r)ψ(r) with respect to the quasi-stationary scattering state |f + . From I E (r), the diffuse flux γ (L) (k d ) -i.e. the disorder average of γ(k d ), see Eq. (31), in ladder approximation -of particles scattered into directionk d is finally obtained as: where γ (L) Since we assume that disorder scattering events -represented by the term P E in Eq. (42) -are much more frequent than collision events, we may neglect the spatial dependence of the collision terms in Eq. (42) and approximate them by δ-functions: f E ,E ,E = dr dr f E ,E ,E (r , r , r) .
The transport equation (42) then reduces to: As shown in [51], this equation can also be derived from the nonlinear Boltzmann transport equation. Due to the above collision approximation, the spatial transport of particles in Eq. (48) is solely governed by the propagation P in the disorder potential, whereas the collision terms g and f lead to a redistribution of energies. As compared to Eqs. (40,41), these terms simplify as follows: with k 4 = k 1 + k 2 − k 3 , E 12 = E 1 + E 2 − E k 1 +k 2 /2, and |k 12 , |k 34 as defined after Eq. (18).

Thermalization
For a given form of the two-body T -matrix (e.g. s-wave scattering, see below), we can now calculate the collision terms g and f according to Eqs. (49,50), and then numerically solve the transport equation (48) by iteration. Before presenting the corresponding numerical results in Sec. 6, however, we will discuss, in the remainder of this section, some general properties of g and f , which, as shown below, lead to thermalization of the single-particle energies for an infinite system. As shown in Appendix D, the collision terms fulfill the following relations: Both relations follow from the fact that the T -matrix associated to the atom-atom interaction potential U (r) fulfills the optical theorem, Eq. (19), and express conservation of the particle and the energy flux, respectively. Moreover, from Eq. (41), one can show that: This equation expresses microscopic reversibility of the collision dynamics: given two particles with energy E 1 and E 2 , the collision process E 1 , E 2 → E 3 , E 4 occurs with the same probability as the reverse process E 3 , E 4 → E 1 , E 2 given two particles with energy E 3 and E 4 . The square roots in the denomimators of (53) result from the traces over the undetected particle with energy E 4 = E 1 + E 2 − E 3 (left-hand side) or E 2 (right-hand side), respectively. Using Eqs. (51,52) -and the fact that the linear propagator P E (r, r ) = P (r, r ) is independent of E -it follows that the quantities J(r) = ∞ 0 dE J E (r) and K(r) corresponding to the particle and energy flux, respectively, both fulfill the same linear transport equation: where the source terms J 0 (r) = √ E i I 0 (r) and K 0 (r) = E i √ E i I 0 (r) differ only by the constant factor E i . Due to Eqs. (51,52), the collision terms drop out from Eq. (48) when integrating over E. Since the linear transport equation fulfills flux conservation, this, in turn, implies that, both, particle and energy flux are conserved. Furthermore, since K 0 (r) = E i J 0 (r), the same relation holds for the solutions of the linear equations (55,56): Figure 7. Example of a crossed diagram contributing to coherent backscattering in the case of three interacting particles. It is obtained from the ladder diagram shown in Fig. 4 by reversing the direction of propagation of the dashed propagators along the path r 1 → r 3 → r 4 → r 5 → r 7 → r 8 → r 9 .
After these preparatory steps, we can now look for a solution of the transport equation (48) in case of a semi-infinite medium. Far away from its boundary, I E (r) = I E should become independent of r, and I 0 (r), Eq. (43), tends to zero. Hence, the constant solution I E must fulfill: Using Eqs. (51,53), one can show that I E = √ Ee −γE fulfills Eq. (58) for γ > 0. The constant γ, in turn, is determined by Eq. (57) as γ = 2/E i . Hence the normalized particle flux distribution is given by: This corresponds to a Maxwell-Boltzmann distribution the temperature of which is determined by the initial energy (E i = k B T /2). Thereby, we have demonstrated thermalization in case of a semi-infinite medium. In Sec. 6, we will study the transport behaviour predicted by Eq. (48) for a finite medium, and see how the thermal distribution is approached during propagation through a finite slab.

Crossed diagrams
Before turning to the numerical results, however, we will extend the general formalism of Sec. 4 in order to calculate the leading interference correction (in the weak disorder parameter 1/(k dis )) to the average scattered flux density. This correction is described by crossed diagrams [52], which are obtained from the ladder diagrams by reversing the direction of propagation of a single amplitude. Starting from the ladder diagram shown in Fig. 4, we can construct, for example, the crossed diagram shown in Fig. 7. It amounts to an interference between two amplitudes where the detected atom is emitted from r 9 and r 1 , respectively. For a given wavevector k d of the detected atom, the backscattering angle θ is defined by cos θ = −k i ·k d . Since annihilation and creation of atoms with wavevector k by the disorder potential at position r are associated with factors e ±ik·r , respectively, see Eqs. (4,15), this leads to a phase factor e iq·(r 1 −r 9 ) with respect to the ladder diagram, Fig. 4, where Since r 1 and r 9 refer to randomly chosen positions of scattering events, this phase factor vanishes on average unless q 0 ⇔ k i −k d , corresponding to exact backscattering (θ = 0). Therefore, this effect of interference between reversed amplitudes is called 'coherent backscattering' [8,9,10]. More precisely, one can show that the angular width of the coherent backscattering interference peak is approximately given by ∆θ 1/(k dis ) [53]. For a single particle, the height of this peak at θ = 0 equals the incoherent background as described by the ladder diagrams (except for single scattering which only contributes to the background), what amounts to an enhancement of the backscattered flux by a factor 2. We will show below how this enhancement factor changes as a consequence of elastic and inelastic atom-atom collisions.
For this purpose, we will derive a transport equation for the 'crossed density' which describes a pair of amplitudes propagating in opposite directions. In Fig. 7, the corresponding crossed scattering path is given by r 1 → r 3 → r 4 → r 5 → r 7 → r 8 → r 9 (where we define the direction of the path to be fixed by the solid arrows, whereas the dashed arrows propagate in the opposite sense). The remaining parts (r 2 and r 6 ) correspond to ladder diagrams already treated in Sec. 4. Due to energy conservation, the energies E and E associated with the two counterpropagating conjugate amplitudes always fulfill the following relation:

Crossed building blocks
This leads us to the first crossed building block: describing single-particle propagation with different energies E (wave vector k 1 ) and E (wave vector k 2 ) for the conjugate amplitudes, see Fig. 8a). For E = E, i.e. E = (E i + E d )/2 due to Eq. (61), it reduces to the ladder propagator P E , see Eq. (39).
The following building block, Fig. 8b) ,  65), obtained by reversing the solid arrow between r 1 and r 3 for g E1,E2 , Fig. 6b). e) Conjugate crossed collision h (C) * E2, E1 . Note that, due to different possibilities of reversing single-particle amplitudes, there are more crossed than ladder building blocks (see Fig. 6).
Similarly, Fig. 8c), where and |k 12 , |k 34 as defined after Eq. (18), represents the crossed counterpart of inelastic collision f E 1 ,E 2 ,E 3 , see Fig. 6c). It reduces to two times the corresponding ladder term, Eq. (50), for E 2 = E 3 = E 2 = E 3 . The factor 2 originates from the fact that we can reverse the single-particle amplitudes of the ladder building block, Fig. 6c), also in a different way (with the outgoing dashed arrow pointing to r 1 instead of r 2 ) giving rise to an identical term. The wavevectors k 1 , k 2 and k 3 in Eq. (64) are associated with Green's functions emitted from (or pointing towards) positions r 1 , r 2 and r 3 , respectively.
Similarly, there also exist two different possibilities for reversing the ladder building block g E 1 ,E 2 , Fig. 6b). Apart from g (C) E 1 ,E 2 , see Eq. (63) and Fig. 8b), this gives rise to a new building block, Fig. 8d): where |k 11 = |k 1 + |−k 1 / √ 2. The corresponding conjugate diagram, see Fig. 8e), is given by h . Note that the two colliding particles exhibit opposite wavevectors (k 1 and −k 1 ), and therefore the energy of the collision event is fixed toẼ 1 +E 1 = E i +E d due to Eq. (61). Each of the above crossed building blocks exhibits an incoming and an outgoing crossed density (defined by the direction of the solid arrow, as mentioned above). Additionally, the two-particle building blocks, Figs. 8b-e), exhibit an incoming ladder density. The latter is given by the solution I E (r) of the ladder transport equation (48).

Crossed transport equation
Propagation of the crossed density can now be described by an integral equation accounting for all possible combinations of the above crossed building blocks (see Fig. 8). An example is displayed in Fig. 9a). Here, the outgoing crossed density of the building block shown in Fig. 8d) serves as the incoming crossed density for the building block shown in Fig. 8e). The resulting combination, Fig. 9a), exhibits the following remarkable property: if we look at the outgoing arrows (solid arrow pointing to r 2 , dashed arrow pointing to r 1 ) corresponding to the detected particle, we see that the detected particle exhibits no collision with the other particles involved in Fig. 9a). The evolution of these undetected particles therefore has no impact on the detected particle and, consequently, as discussed at the end of Sec. 3.2, the process shown in Fig. 9a) may be disregarded when calculating the detection signal. The same remains true if -instead of attaching Fig. 8e) directly to Fig. 8d) -an arbitrary sequence of the remaining crossed building blocks, Figs. 8a), b) or c), is inserted in between. In contrast, for any other combination of building blocks, e.g. Fig. 9b), all involved particles turn out to be connected to each other (through collision events or partial traces), thus contributing to the propagation of the crossed density.
In order to exclude combinations of the former type from the transport equation, we split the crossed density into two parts, i.e. C E (r) = C    E (r). According to the rules mentioned above, the building block Fig. 8e) is excluded from the transport equation for C (2) E (r). In total, the transport equations therefore read as follows: Figure 9. Combinations of crossed building blocks. a) When attaching Fig. 8d) to Fig. 8e), the detected atom (solid arrow pointing to r 2 , dashed arrow pointing to r 1 ) exhibits no collision with the undetected atoms. Combinations of this type therefore do not contribute to the detection signal, and must be excluded from the crossed transport equation. b) In contrast, the inverse combination, i.e. attaching Fig. 8e) to Fig. 8d), contributes to the propagation of the crossed density, and must be taken into account in the transport equation. and In Eq. (66), the incoming crossed density is given by: where q = k i + k d , see Eq. (60), where the wave vector of the detected particle is determined by the energy E d and the position R of the detector (in the far field) as k d = √ E d R/R, and z r,−k i (or z r,k d ) corresponds to the distance an incoming (or outgoing) particle travels inside the scattering region, as defined in Eqs. (43,45). Finally, the coherently backscattered flux density results as with associated spectral density where the last term accounts for single scattering. The total average flux measured by a detector placed in directionk d , see Eq. (31), then corresponds to the sum of the ladder and the crossed component,

Numerical solutions of the transport equations
After having developed the general scattering formalism valid for arbitrary shapes of the interaction potential U (r) and the scattering region V, we will now focus on the case of a short-range potential U (r) and a slab geometry for V. As explained in Appendix E, the T -matrix is then described by a single parameter -the s-wave scattering length a s . The corresponding average distance between (inelastic) collision events is given by: In the following, we measure the interaction strength in terms of the ratio between dis and int : which, as explained in Sec. 4.1, should fulfill the condition α 1, and, due to Eq. (71), is proportional to a 2 s . Indeed, we see from Eqs. (E.2,E.3) that the ladder collision terms g and f both depend on α. The terms proportional to a s in g drop out as a consequence of flux conservation, see Eq. (51). The same is not true for the crossed collision terms g (C) and h (C) , see Eqs. (E.4,E.5), which depend on a second parameter proportional to a s : Since √ E i a s 1 for s-wave scattering, it follows that β α. The parameter β can also be expressed in terms of the healing length ξ = (8πρ 0 a s ) −1/2 [55], i.e. β = dis /( √ E i ξ 2 ), or in terms of the interaction parameter g = 8πa s appearing in the Gross-Pitaevskii equation, i.e. β = gρ 0 dis / √ E i . The Gross-Pitaevskii equation is valid in the limit a s → 0 and ρ 0 → ∞ [56] such that a s ρ 0 = const. Since α → 0 in this limit, our previously derived transport equations for nonlinear coherent backscattering based on the Gross-Pitaevskii equation [27,28,29] must be recovered from Eqs. (48,66,67) for α = 0, as it is indeed the case if we insert the s-wave expressions, Eqs. (E.2-E.6), evaluated at α = 0.
Concerning the geometry of the scattering medium, we choose a slab confined between two planes, z = 0 and z = L, respectively, with perpendicular incident wavevector, i.e. k i = (0, 0, k i ). The thickness of the slab in units of the disorder mean free path defines its optical thickness b = L/ dis . The slab geometry is very convenient from a numerical point of view, since the integration over x and y can be performed analytically in Eqs. (48,66,67), such that the resulting transport equations only depend on z [54]. Moreover, due to rotational symmetry around the z-axis, the backscattered flux g(k d ) = g(θ) depends only on the backscattering angle θ defined bŷ k i ·k d = − cos θ, and the distances appearing in Eqs. (43,45,68,70) simplify to z r,−k i = z and z r,k d = z/ cos θ, respectively. Finally, the integration over the scattering volume V in Eqs. (45,70) reduces to V dr/A → L 0 dz. The one-dimensional versions of the transport equations (48,66,67) can now be solved numerically, e.g. by iteration.  The kink of this distribution is recovered at the beginning of the slab (i.e., for z = 0). Deep inside the slab (i.e., for z = L/4 and z = L), the spectrum collapses onto a thermal Maxwell-Boltzmann distribution with average energy E av = E i (red solid).

Density inside the slab
. We see that, in spite of the weakness of the interaction (α = 1/100), the inelastic component dominates, especially deep inside the slab. This can be explained by the large number (≈ b 2 ) of scattering events required to traverse a slab with thickness b. The expected number of two-body collision events thus results as approximately αb 2 = 16. Let us note that the inelastic component of the flux is associated with a non-condensed fraction of atoms, since an N -fold product of a single-particle state (as required from the formal definition of a condensate via the stationary one-particle density matrix [56]) with fixed total energy implies fixed energies also for the individual particles.
The normalized energy distribution J (inel) E (z) of the inelastic component is shown in Fig. 10b), for different positions z inside the slab. We see that, far inside the slab, i.e. at z = 10 dis (green dashed) and z = 40 dis (black long-dashed), the energy distribution approaches a Maxwell-Boltzmann distribution J (59). Thereby, we confirm the analytical result derived in Sec. 4.4 for an infinite slab. In contrast, at the beginning of the slab (blue dotted), the distribution is not yet thermalized, and lies between the Maxwell-Boltzmann distribution and the distribution √ obtained after a single inelastic collision event (thin black line) and normalized according to Eq. (51).  Fig. 10). As already known from previous work on nonlinear coherent backscattering in the purely elastic Gross-Pitaevskii limit [27,28,29], the backscattered flux γ (C) (0) for α = 0 (red line) exhibits a transition from constructive to destructive interference for increasing β. This transition can be explained by the fact that the nonlinearity effectively introduces dephasing between two reversed scattering paths [27,28,29]. For larger β, the phase difference accumulated along the shortest scattering paths, i.e. those exhibiting only a few, but at least two scattering events (since, as mentioned above, single scattering does not contribute), may exceed π/2, and thus lead to destructive interference. This behavior changes if a finite amount of inelastic scattering, α = β/10 corresponding to a s √ E i = 1/10, see Eq. (73), is taken into account (blue line): At first, the interference drops faster than for α = 0, since inelastic scattering changes the frequency and thus leads to an additional dephasing mechanism. At larger values β, however, the decrease of the backscattered flux is slowed down as compared to the purely elastic case, such that a transition to destructive interference is not observed in Fig. 11a).
fact that the maxima of the peaks are slightly shifted with respect to E i , for which, at present, we are lacking a simple explanation. Please note, however, that the interference peaks exceed the background (most remarkably for larger values of β), corresponding to a coherent backscattering enhancement factor of the inelastic flux contribution η (0) > 2 within a narrow spectral window around the energy E d where the crossed contribution is maximal. This enhancement is a consequence of the many-wave interference character of nonlinear coherent backscattering [57,58] resulting from the fact that, as discussed in Sec. 5.2, there are several ways of reversing the scattering paths when constructing crossed from ladder diagrams. The number of these possibilities increases with increasing number of inelastic scattering events. Although, due to the small width of these peaks, the total inelastic component γ (C,inel) (0) (integrated over E d ) turns out to be smaller than the background γ (L,inel) (0) [cf. the dashed blue and dotted green lines in Fig. 11a)], this many-wave interference effect contributes to the above observed slowing down of the decrease of the backscattered flux.

Conclusion
Within this paper, we have derived a microscopic N -body scattering theory for interacting particles in a weak disorder potential in three dimensions. We have applied this diagrammatic theory to a stationary scattering scenario for an asymptotically noninteracting quasi-plane matter wave incident on a three-dimensional slab, with the disorder potential and inter-particle collisions confined to the slab region, and hereby verified the viability of our theory to address, on the one hand, very fundamental but, on the other hand, very timely questions of quantum transport for interacting particles in random environments. In a clear and precise manner, we demonstrated how one can bridge the gap between strictly unitary many-body evolution and its implications on the mesoscopic level governed by a transport equation similar to the nonlinear Boltzmann transport equation. Furthermore, we have determined the coherent corrections due to the wave nature of the particles leading to the effect of coherent backscattering. We have demonstrated that inelastic scattering slows down the decrease of the coherent backscattering peak as compared to the purely elastic case described by the Gross-Pitaevskii equation.
Let us briefly summarize the basic assumptions of our theory: first, we assume an optically thick scattering medium (b = L/ dis 1) allowing for multiple scattering in a weak disorder potential, where the mean free path is much larger than the wavelength of the incoming particles ( dis √ E i 1). Second, the interaction strength as, respectively, quantified by the parameters α and β for ladder and crossed collisions, see Eqs. (72,73) for the case of s-wave scattering, should fulfill the condition α, β 1, such that atomatom collisions occur less frequently than disorder scattering events. With view at future work, we expect that the condition β 1 can be dropped by summing the corresponding diagrams, see Fig. 8b,d,e), without intermediate disorder scattering [28]. If β 2 dis √ E i , another type of collision process -corresponding to scattering induced by the fluctuating background density -sets in which is described by diagrams similar to our ladder diagrams [59]. Furthermore, relaxing the contact approximation for the collision terms, Eqs. (46,47), will allow to determine the effect of attractive or repulsive interactions onto the spatial atomic density profile. It will hence be a worthwhile task to extend our theory to stronger interactions, although we surely expect to encounter certain limits (e.g., the regime of superfluidity) where other methods will be required.
Concerning an experimental verification of our results, the application to a stationary scattering setup with matter waves constitutes, on the one hand, a very timely scenario, as, e.g., the developments of atom lasers and matter wave interferometers on atom chips [22,23,60] rapidly progress. On the other hand, many years of expertise have been gathered within the field of wave-packet spreading upon releasing the condensate Figure A1. The dashed lines split the diagram of Fig. 1 into 4 subdiagrams (i), (ii), (iii), and (iv). Note that the subdiagrams (i) and (ii) -and likewise (iii) and (iv) -are not connected to each other. This allows us to factorize the 3-particle diagram into 1and 2-particle diagrams.
from a trap into a new environment, where, e.g., the first experiments on coherent backscattering of (non-interacting) matter waves have been reported recently [11,12]. Consequently, an extension of our theory to time-dependent scenarios based on recent progress in this field [51,61,62] presents a significant and feasible task. Similarly, also a finite correlation length of the disorder potential can be taken into account, see [24,25] for the non-interacting case.
In conclusion, we are confident, that our present theory and the rather straightforward extensions discussed above will substantially foster a more complete understanding of quantum transport under the interplay of disorder and inter-particle interaction and can contribute to a unifying picture from microscopic to macroscopic scales.

Appendix A. Factorization of the transition amplitude
In this appendix, we show how an arbitrary N -particle scattering diagram can be factorized into single-particle propagators and two-body collisions. We first look at the example diagram shown in Fig. 1. As shown in Fig. A1, this diagram can be split into four independent subdiagrams. The two subdiagrams connected to the initial state -(i) and (ii) in Fig. A1 -correspond to Møller operators, and the remaining ones -(iii) and (iv) -to Green's operators. This gives rise to the following matrix elements: + (E) = 1 2 dp 2 dp 3 (2π) 6 p 4 , p 7 |T U (E)|p 2 , p 3 p 2 , p 3 |Ω (V ) (E)|k 2 , k 3 ,(A.2) G (iii) (E) = 1 4 dp 1 dp 5 dp 6 dp 8 (2π) 12 k 5 , k 6 |Ĝ V (E)|p 5 , p 6 (A.3) × p 5 , p 6 |T U (E)|p 1 , p 8 p 1 , p 8 |Ĝ V (E)|p 7 , p 9 , Note that the diagrams (i) and (ii) are not connected to each other in Fig. A1. The corresponding Møller operators can therefore be factorized as in Eq. (26). Likewise, the Green's functions corresponding to (iii) and (iv) are factorized according to Eq. (27). The prefactors 1/2 and 1/4 in Eqs. (A.2,A.3) originate from symmetrization in the twoparticle subspace (e.g. the states |p 2 , p 3 and |p 3 , p 2 are identical and therefore must not be summed over twice). It turns out that these factors are compensated, however, by the two possibilities to associate the initial and final single-particle states with each other in Eqs. (26,27).
The total transition amplitude results as: k 4 , k 5 , k 6 |Ω ( fig.1) + (3E i )|k 1 , k 2 , k 3 = 1 2 dp 4 dp 7 dp 9 (2π) 9 Ω (i) Now, we again apply Eqs. (26,27) to factorize the two-particle Møller and Green's operators on the right-hand side of Eqs. (A.2,A.3) into single-particle operators. In this way, we recover most of the terms in Eq. (28). The only ones which appear to differ from Eq. (28) are those associated to k 1 , p 1 , p 7 and p 8 , which we reformulate as follows: 1 2 dp 9 (2π) 3 p 1 , p 8 |Ĝ V (3E i − E 4 )|p 7 , p 9 p 9 |Ω (V ) where we again applied Eq. (27), used the completeness relation dp 9 |p 9 p 9 | = (2π) 3 , and defined: Note that G(−E 1 ) is a complex analytic function of E 1 with poles only in the upper half of the complex plane. This, again, is due to the fact thatĜ V (E) as a function of E exhibits poles only in the lower half, whereas E 1 enters with negative sign in the right hand side of Eq. (A.7). We now reformulate some terms in Eq. (A.6) as follows: where we used Eq. (8), the alternative but equivalent expressionĜ V (E) =Ĝ 0 (E) + G V (E)VĜ 0 (E) with respect to Eq. (11), and the identitŷ resulting from 1 ab = 1 b−a 1 a − 1 b (where we set the imaginary part in the denominator ofĜ V (E 1 ), see Eq. (9), equal to 2 instead of , and used the fact thatĜ V (E 1 ) and G V (E i ) have the same set of eigenvectors). After inserting Eq. (A.8) into Eq. (A.6), we perform the integral over E 1 by closing the integration contour in the lower half of the complex plane. (Remember that G(−E 1 ) has no poles in the lower half!) Thereby, the energy E 1 is set to E i , and we finally recover the missing terms in Eq. (28): The above procedure can be generalized to an arbitrary many-particle scattering diagram: We first divide the whole diagram into independent subdiagrams. Then, some of these subdiagrams turn out to be connected to each other by single-particle propagators. In the above example, Fig. A1, this is the case for the subdiagram (i) and (iii), which are connected by the single-atom propagators from k 1 to p 9 with energy E i and from p 9 to p 1 with energy E 1 . We have to show that these propagators merge into a single propagator (from k 1 to p 1 with energy E i ). For the case that one of the propagators is connected to the initial state (and thus corresponds to a Møller operator), the corresponding general identity is given by Eq. (A.10). If both propagators correspond to Green's operators, e.g.Ĝ V (E 1 ) andĜ V (E 2 ) below, the required identity is proven as follows: where we again used Eq. (A.9) and the fact that G (1) (−E 1 ) and G (2) (−E 2 ) (which correspond to arbitrary other subdiagrams where the energy E 1 or E 2 enters with negative sign) exhibit no pole in the lower half of the complex plane. Note that the term withĜ V (E 1 ) in the second line of Eq. (A.11) vanishes after integrating over E 2 , since no pole remains in the lower half. In total, the concatenation of two Green's operators, i.e.Ĝ V (E 1 )Ĝ V (E 2 ) on the left-hand-side of Eq. (A.11), reduces to a single Green's operator, i.e.Ĝ V (E 2 ) on the left-hand-side, whereas the energy E 1 is set equal to E 2 .
Hence, their evolution factorizes from the one of the detected particle and need not be taken into account. The prefactors are then generalized as follows: 1/2 → N (N − 1)/4 in Eq. (B.1) and 1 → N (N − 1)/2 in Eq. (B.2). Let us now compare these prefactors with the ones obtained from the iterative procedure based on the connection of building blocks in Secs. 4 and 5: The factors N (N − 1) N 2 (for N 1, since N → ∞ in the quasi-stationary limit) are accounted for by the source term ρ 0 in Eq. (43), which is proportional to N , see Eq. (21), and occurs two times for a two-particle process proportional to the density squared. What remains is a factor 1/2 for each collision event [twice in Fig. 2a) and once in Fig. 2b)] which is included in the definition of the building blocks, Eqs. (40,41). The origin of this factor can be traced back to the indistinguishability of bosonic particles. Indeed, as argued at the end of Sec. 3.2, all factors related to indistinguishability finally drop out in the case where all particles are initially in the same state. Since the T -matrix for indistinguishable particles, see Eq. (18), differs by a factor 2 from the one for distinguishable particles, this must be counterbalanced by the above factor 1/2.