Quench dynamics of a weakly interacting disordered Bose gas in momentum space

We theoretically study the out-of-equilibrium dynamics in momentum space of a weakly interacting disordered Bose gas launched with a finite velocity. In the absence of interactions, coherent multiple scattering gives rise to a background of diffusive particles, on top of which a coherent backscattering interference emerges. We revisit this scenario in the presence of interactions, using a diagrammatic quantum transport theory. We find that the dynamics is governed by coupled kinetic equations describing the thermalization of the diffusive and coherent components of the gas. This phenomenon leads to a destruction of coherent backscattering which can, at short time, be approximated by an effective dephasing mechanism, whose rate is controlled by the particle collision time. These predictions are well reproduced by numerical simulations.


I. INTRODUCTION
When perturbed from an equilibrium situation, isolated many-body systems generally experience a thermalization process and eventually return to equilibrium at sufficiently long time [1].This process arises because the interacting system serves as a "bath" for all its subparts, the final state being characterized by a Gibbs ensemble.A number of works have explored the formation of this thermalized state, with special attention dedicated to the dynamical emergence of a Bose condensate [2][3][4].The out-of-equilibrium dynamics leading to thermalization can also follow a rich variety of scenarios.Those have recently raised considerable interest in the cold-atom community, where the conditions of truly isolated quantum gases can be achieved at an unprecedented level.Non-integrable systems, for instance, usually display an intermediate "prethermal" stage where the system evolves rather slowly and looks approximately thermalized [5,6].Prethermalization can be modeled by a generalized Gibbs ensemble, characterized by a small set of parameters [7].The many-body dynamics may also exhibit universal scaling properties when the system is quenched through or in the vicinity of a quantum phase transition [8,9], or when it is initially prepared in a far off-equilibrium state [10][11][12][13].
Much less is known about the out-of-equilibrium of interacting disordered systems.At the many-body level, the competition between disorder and interactions may lead to many-body localization (for recent reviews see [14,15]), initially addressed in the context of electron conduction [16].A consequence of many-body localization is the absence of thermalization.When quenched out-of-equilibrium, many-body disordered systems may or may not reach a universal thermalized state, depending on the magnitude and properties of the disorder and interactions.In weakly interacting Bose gases, which can be described at a mean-field level with a nonlinear Schrödinger equation, many-body localization does not occur and thermalization is the rule.Notwithstanding the relative simplicity of this limit, the combination of weak interactions and disorder still gives rise to a number of puzzling phenomena, such as thermalization via weakly coupled localized states [17,18] or subdiffusive spreading of wave packets [19][20][21][22][23][24].
In this article, we study the interplay between disorder and interactions in a dilute Bose gas, in the limit of weak disorder.This regime has been under the focus of a number of cold-atom experiments probing, e.g., one-dimensional Anderson localization [25,26], coherent backscattering (CBS) [27][28][29] as well as its control over external dephasing [29,30].CBS of cold atoms, in particular, was probed in an optics-like configuration where an ultracold Bose gas was initially given a finite mean velocity, and its subsequent dynamics in the presence of disorder probed in momentum space.This configuration, originally introduced in [31], turned also useful to explore other interference phenomena like coherent forward scattering [32], to achieve an echo spectroscopy of coherent transport in disorder [30,33] or to monitor the thermalization and dynamical formation of condensates in momentum space [34].In most of these works, disorderalbeit weak -was the main ingredient driving the atomic dynamics, so that interactions could be neglected in first approximation.As is well known, however, particle interactions generally affect significantly coherent transport and, in particular, coherent backscattering.This question was previously addressed in the context of nonlinear optics of continuous beams [35][36][37][38] or of atom lasers [39].A theoretical description of particle interactions in an out-of-equilibrium regime where mesoscopic effects like CBS occur is, on the other hand, still absent.It is the goal of the present work to fill this gap.
Following [27,28,30,31], we consider the out-ofequilibrium dynamics of a two-dimensional, weakly interacting Bose gas initially prepared in a plane-wave state with finite velocity in a random potential.In this config-uration, the momentum distribution quickly acquires a ring-shape profile associated with classical particle diffusion, with an interference, CBS peak emerging on the top.In the presence of interactions, this picture slowly evolves in time as the whole distribution thermalizes.Well before thermalization fully develops however, we observe a rather fast contrast loss of the CBS peak.To explain this phenomenon, we develop a microscopic, diagrammatic perturbation theory of coherent particle transport including both disorder and interactions.This allows us to derive coupled kinetic equations for the diffusive and coherent components of the momentum distribution.By solving the latter numerically at short time, we achieve a precise description of the time-evolution of CBS in the presence of interactions.Although the fundamental mechanism responsible for the decay of the CBS peak is thermalization, we show that it can be approximately described in terms of an effective dephasing effect, of rate controlled by the particle collision time.
The article is organized as follows.In Sec.II we formulate the problem and illustrate it through a numerical simulation.In Sec.III, we recall the main elements of quantum transport theory in disorder when interactions are neglected.This approach is extended to the interacting regime in Secs.IV and V, and confronted with numerical simulations in Sec.VI.Main results are summarized in Sec.VII, and technical details are collected in two appendices.

II. MOMENTUM-SPACE DYNAMICS
We consider the out-of-equilibrium evolution of a disordered interacting Bose gas.Interactions, assumed weak, are treated at the mean-field level, on the basis of the Gross-Pitaevskii equation (GPE) for the Bose field Ψ(r, t).In Eq. ( 1) and in the following, we set = 1.From now on we focus on the two-dimensional geometry, although most results of the article are valid in dimension 3 as well.V (r) is a random potential, assumed to follow a Gaussian statistics with zero mean and no correlation : where γ > 0 sets the disorder strength and the overbar refers to configuration averaging.The assumption of uncorrelated disorder does not imply any loss of generality.The mean free time (defined below) being the only relevant disorder parameter for the dynamics, the results of the paper hold as well for a short-range correlated disorder.To illustrate the problem we are interested in, we first study the numerical propagation of an initial plane wave φ(r) ≡ r |Ψ(t = 0) = √ ρ 0 exp(ik 0 • r) in the random potential (ρ 0 is the particle density), by computing the wave function Ψ(r, t) at different times on a discretized regular grid of 200 × 200 sites with step a, for periodic boundary conditions.The temporal propagation is performed using a split-step algorithm of time step ∆t, alternating propagations of the linear part of the GPE, exp −i(−∇ 2 /2m + V (r))∆t , and of the nonlinear part, The linear part of the evolution operator is expanded in a series of Chebyshev polynomials, as described in [40][41][42][43].From the wave function, we compute its Fourier transform from which we infer the disorder-averaged momentum distribution, |Ψ(k, t)| 2 , normalized according to The computed momentum distribution is shown in Fig. 1 at three different times.We choose as unit of length the discretization step a in FIG. 1. Momentum distribution |Ψ(kx, ky, t)| 2 computed numerically at three successive times using the GPE (1), starting from a plane wave of initial momentum k0 = (π/5, 0) in a Gaussian, uncorrelated random potential.The nonlinear term in the GPE leads to an early-time decay of the CBS peak.Here gρ0 = 0.004 and γ = 0.038.Data are averaged over about 14000 disorder realizations.our numerics.In order for the discretization to be a good approximation to the continuous equation ( 1), one must simply ensure that the de Broglie wavelength is significantly larger than the grid spacing, i.e. 2π/k 0 a.In the following we typically use k 0 a = π/5, so that discretization effects are small.In Fig. 1 and in all other simulations based on Eq. (1), we express momenta in units of 1/a, lengths in a, times in ma 2 and energies in 1/(ma 2 ).The disorder amplitude, γ, is then in units of 1/(m 2 a 2 ).In Fig. 1, the three times shown are given in units of the mean free time τ , i.e. the typical collision time on the scattering potential.In the Born approximation, τ = 1/(mγ) [see Eq. ( 11) below].We choose τ 26.3, so that the product of k 0 with the mean free path ≡ k 0 τ is k 0 10.4,i.e. much larger than one.18) and (19), respectively.They both saturate at a value slightly above 100, which is close to the analytical estimate τ /(πν 0 ) 98.9 for a numerically computed mean free time τ 51.84 and on-shell density of state ν 0 0.167.This is the so-called limit of weak disorder, where, in the absence of interactions, the momentum distribution is essentially the sum of two contributions.The first is a background of diffusive particles scattered elastically on the random potential.In Fig. 1 this contribution manifests itself as a ring of radius k 0 .The second is a narrow interference peak centered around k = −k 0 , the coherent backscattering (CBS) effect.CBS in this configuration was first described theoretically in [31], and experimentally measured with cold atoms in [27,28,30].
In the presence of weak interactions, a main change compared to this picture is a decay of the CBS peak amplitude at short time.This decay appears as soon as interactions are nonzero, as shown in the lower panel of Fig. 2 for even weaker values of gρ 0 .Also shown in the upper panel is the height of the diffusive ring, which decays as well, albeit more slowly than the CBS peak.This results in a decay of the CBS contrast, well visible in Fig. 1.In the simulations of Figs. 1 and 2, the magnitude of interactions is chosen such that the mean free path associated with the collisions on the nonlinear potential g|Ψ| 2 is larger than (this condition will be clarified in Sec.IV).In the rather small time-range of these figures, this nonlinear potential thus plays the role of a perturbation for the dynamics, which remains mainly governed by the disorder.
The observed decays of the diffusive background and CBS peak constitute the early-time manifestations of a thermalization process of the whole momentum distribution.The long-time thermalization properties of the diffusive background have been previously addressed in [34].In the sequel of the paper, we provide a theoretical basis for the formalism used in [34], and go one step forward by constructing a kinetic theory which encompasses both the incoherent diffusive component and the CBS contribution (neglected in [34]).Equipped with this theory, we then reproduce and explain the temporal evolutions observed in Figs. 1 and 2, and in particular the decay of the CBS peak, which turns out to resemble the one associated with a decoherence process.

III. LINEAR REGIME: THEORY
The theory of diffusion and CBS in momentum space has been presented in [31].Here we only recall the main elements required to introduce the nonlinear diagrammatic theory in the next sections.We also adopt a slightly different point of view than in [31], focusing more on the energy distribution of the particles, which plays a major role in the thermalization process at work when interactions are nonzero.
When g = 0, the momentum distribution can be expressed as [31,44,45] where the density kernel is defined in terms of the energydependent, retarded and advanced Green's operators G R/A and of the initial state |φ as In the following, we will also work with the disorderaveraged occupation number f (t), defined as where ν is the density of states per unit volume at energy .As a consequence of particle conservation, this quantity is normalized according to: This condition identifies the product ν f (t) as the energy distribution of the Bose gas.8).Upper solid lines with arrows refer to k|G R +ω/2 |k0 , and lower dashed lines to k0|G A −ω/2 |k .Note that for solid (dashed) lines, arrows coincide with the (opposite of the) direction of propagation.Dotted vertical lines symbolize the correlation function in Eq. ( 2).(b) Series of diagrams describing the interference between time-reversed paths, responsible for CBS ("crossed diagrams").The second equality follows from time-reversal invariance.(c) In the long-time limit t τ where low scattering orders have a negligible weight, the series of ladder and crossed diagrams coincide.

A. Diffusive background at long time
From now on, we focus on the case where the initial state is a plane wave, |φ = |k 0 .We also assume disorder to be weak, k 0 1, so that perturbation theory can be used.The main contribution to I ,ω (k) is then given by the series of ladder diagrams ("Diffuson"), which we denote by I D ,ω (k).The latter obeys the Bethe-Salpeter equation [31,44,45] which is shown diagrammatically in Fig. 3a.The average Green's function is given by In the weak-disorder limit, the self energy Σ( , k) can be evaluated by perturbation theory.We restrict ourselves to the leading-order contribution provided by the Born approximation, which coincides with the Fermi golden rule: where B(k) = γ is the Fourier transform of the disorder correlation function (2).Equation (10) defines the mean free time, with, in the Born approximation, ν = m/(2π).In the long-time limit t τ , the contribution to the momentum distribution (4) due to I D ,ω can be obtained in the following way.First we integrate Eq. ( 8) over k, and expand the second term in the right-hand side for ωτ 1 (hydrodynamic limit).This yields where we used that G We thus find, according to Eq. ( 6), for the occupation number, where the superscript (0) refers to the non-interacting limit and we introduced the spectral function: Second, we integrate Eq. ( 8) over and take the Fourier transform with respect to ω.Using the result (13), we infer Eq. ( 15) is an isotropic function of k, centered at |k| = k 0 , which corresponds to the diffusive ring in Fig.
The radial profile of this ring is essentially the one of the energy distribution ν f (0) = A (k 0 ), which for g = 0 coincides with the spectral function.In the absence of interactions, this energy distribution does not change in time.This is of course expected, as the only process at work is elastic multiple scattering, which does not involve any energy transfer.

B. Coherent backscattering at long time
The CBS contribution is deduced from the diffusive one by making use of time-reversal invariance.The CBS peak stems from the interference between time-reversed multiple scattering paths, described by the diagrammatic series in Fig. 3b ("crossed diagrams").In the long-time limit t τ where low scattering orders have a negligible weight, the ladder and crossed series exactly coincide at k = −k 0 due to time-reversal symmetry (see Figs. 3b  and c): such that where 0 ≡ k 2 0 /(2m).In the absence of interactions, the diffusive and CBS amplitudes at −k 0 thus coincide, see Eqs. ( 15) and ( 17), and are independent of time.The full k dependence of the CBS profile can be calculated as well, as was done in [31].In the rest of the article however, we will essentially focus on its amplitude, |Ψ C (−k 0 , t)| 2 .

C. Full time evolution
Eqs. ( 15) and ( 17) have been obtained in the regime of long times t τ , where low scattering orders can be neglected.While an exact calculation of the diffusive and CBS contributions at any time is a difficult task in general (see, for instance, [46] where this problem was tackled for a speckle potential), for the particular model of uncorrelated disorder we have found that the Bethe-Salpether equation can be exactly solved, giving: for the diffusive background, and for the amplitude of the CBS peak, with ρ max ≡ d A 2 (k 0 )/ν τ /(πν 0 ).Note that, as expected, the CBS and diffusive amplitudes coincide at long time, but not at short time t ∼ τ where low scattering orders -described by the terms within square brackets -come into play.Eqs. ( 18) and ( 19) are shown in Fig. 2 (dashed curves) on top of the results of numerical simulations for g = 0.The agreement is very good at all times.

IV. INTERACTING DIFFUSIVE PARTICLES: THEORY
We now turn to the case g = 0.While an exact theory accounting for both interactions and disorder is a formidable task, even at the level of the GPE (1), relatively simple results can be obtained when interactions are "weak" compared with the disorder.Indeed, in this regime the effect of interactions can be treated as a perturbation of the series of crossed and ladder diagrams.
This approach was previously used in [35][36][37]39] to describe the stationary coherent backscattering effect of continuous waves in finite media, and in [47][48][49] to model the dynamics of interacting wave packets in the diffusive limit.The latter configuration was later extended to the localization regime in [22,23], but by taking into account first-order corrections in g only, see Sec.IV B, while neglecting second-order corrections responsible for thermalization.
In this section and the next one, we develop a quantum transport theory describing the effect of interactions on both the diffusive and CBS signals in the dynamical scenario of Fig. 1.We show that, because the average density of the Bose gas is uniform, linear corrections in g reduce to an irrelevant shift of the mean energy, so that the physics in momentum space is governed by secondorder corrections.The latter are responsible for two coupled thermalization processes of the diffusive and CBS components.

A. Weak interactions
A treatment of the nonlinear potential in Eq. ( 1) as a perturbation of the ladder and crossed series requires that scattering events on the nonlinear potential g|Ψ(r, t)| 2 are rare compared to scattering events on the random potential V (r).In terms of time scales, this condition reads τ NL τ , where τ NL is the particle collision time.To estimate this quantity, we use the Fermi golden rule where we introduced the power spectrum of the nonlinear potential, To leading approximation, the density-density correlator is not modified by interactions, and reads [50] |Ψ(r, t) By defining the mean free path for particle collisions as NL ≡ k 0 τ NL /m we then rewrite the criteria of rare particle collisions and weak disorder as In the following we will assume these conditions to be fulfilled.

B. Leading-order nonlinear corrections
When g = 0, the notion of Green's function can no longer be utilized to express the momentum distribution, as we did in Eq. ( 5).Nevertheless, Eq. ( 4) can still be written, with the density kernel defined as where Our diagrammatic quantum transport theory in the presence of interactions is constructed from the Lippmann-Schwinger equation associated with the GPE (1): where G 0 (k) = ( − k 2 /2m + i0 + ) −1 is the free-space (retarded) Green's function.Iteration of Eq. ( 26) leads to an expansion of Ψ in powers of V and g known as the Born series.In addition to the usual scattering processes on the random potential, this Born series also involves particle collisions.These two elementary processes are diagrammatically shown in Fig. 4, together with the conservation rules for energies and momenta.
FIG. 4. The Born series obtained by iterating the Lippmann-Schwinger equation (26) generates terms built on scattering processes on the random potential V (left) and on the nonlinear potential g|Ψ| 2 (right).The solid and dashed lines with arrows refer to the free-space, retarded and advanced Green's functions, respectively.The cross refers to the random potential V and the wavy line to the interaction parameter g.For each vertex, the lower-right line ( , k) is the outgoing field, and integrations over 1, 2, k , k1 and k2 are understood.
From the Born series for Ψ , one can construct a Bethe-Salpeter equation for the density kernel I ,ω (k).Insofar as particle collisions are less frequent than collisions on the random potential -remember condition (23) -two distinct iterative equations for the diffusive and CBS contributions, I D ,ω and I C ,ω , can still be identified.In this section we first focus on the Bethe-Salpeter equation for I D ,ω .The latter is obtained by adding extra terms to the right-hand side of the iterative equation in Fig. 3a for g = 0, in which the ladder sequence can be interrupted by one or several particle collisions.The leading-order, iterative correction to the Bethe-Salpeter equation for I D ,ω is given by the diagram in Fig. 5a (plus its complex conjugate).It reads where the factor 2 stems from the Wick decomposition of the average Ψ 1 Ψ * 2 Ψ − 1 + 2 Ψ * arising when Eq. ( 26) is multiplied by Ψ * , and the second equality follows from particle conservation, which imposes that From Eq. ( 27), it is easy to show that the contribution of the diagram in Fig. 5a boils down to a constant energy shift −2gρ 0 of G R in the linear Bethe-Salpeter equation (8).Indeed, if we perform the substitution → − 2gρ 0 in the right-hand side of Eq. ( 8) and expand for small g, using we get that the right-hand side of Eq. ( 8) is modified by an extra term which exactly coincides with Eq. ( 27).
In other words, first-order nonlinear corrections to the Bethe-Salpeter equation do not quantitatively affect the diffusive dynamics.Note that this result is in stark contrast with the scenario where one follows the spreading of a wave packet in position space.In this case, diagrams of the type of Fig. 5a were shown to significantly modify the wave-packet density distribution [47][48][49].The difference lies in the behavior of the mean density |Ψ(r, t)| 2 , which evolves in time for an initial wave packet, whereas it always remains uniform for an initial plane wave.
As explained in appendix A, the energy shift obtained here can in turn be described in terms of a modification of the real part of the self energy Σ( , k) appearing in average Green's functions, Eq. ( 9).In addition to this shift, there also exist first-order nonlinear corrections shifting the imaginary part of Σ( , k).These corrections stem from correlations between the disorder and nonlinear potentials in the GPE equation ( 1), but turn out to be very small in the weak-disorder limit.From now, we will thus neglect these self-energy corrections, and always evaluate average Green's functions using Eq. ( 9), with Σ( , k) = −i/2τ .

C. Second-order corrections: thermalization
We now examine interaction corrections to the ladder Bethe-Salpeter equation ( 8) that are of second order in g.Since each vertex g is connected to four field amplitudes (see Fig. 4), these corrections involve six incoming field amplitudes, i.e they are proportional to the third power of the density.Due to the condition (23), we also know that at least one disorder scattering event occurs before every particle collision event.Since the disorder scattering events are described by ladder diagrams (for weak disorder), we group the six incoming arrows into three incoming ladder sequences I D i,ωi , each of them originating from different disorder scattering events.Analyzing all possible non-trivial ways (i.e.those which do not reduce to a mere energy shift) of connecting the incoming arrows to the g vertices, we arrive at the diagrams shown in Fig. 6.The corresponding Bethe-Salpeter equation reads: where we defined for energies, and A closed equation for the distribution f (t) -remember the definition (6) -can be obtained by integrating Eq. ( 30) over k and taking the Fourier transform with respect to ω, in the hydrodynamic limit ωτ, ω i τ 1 (i = 1, 2).The details of the calculation are presented in Appendix B for clarity.They lead to the kinetic equation where we recall that 3 ≡ 1 + 2 − .The interaction kernel is given, in two dimensions, by where K is the complete elliptic integral of the first kind.
The kinetic equation should be complemented by an initial condition, which is here provided by the "coherent mode", i.e. the first term in the right-hand side of Eq. (30).The latter gives (see Appendix B): FIG. 6. Diagrammatic representation of the Bethe-Salpeter equation for the diffusive contribution to the momentum distribution, I D ,ω (k), taking into account second-order interaction corrections (first-order interaction diagrams are discarded, as explained in Sec.IV B).Symbols have the same meaning as in Fig. 5.The numerical prefactors account for the possible combinations of propagation lines connecting to the vertex g.
which is nothing but the occupation number for g = 0, Eq. ( 13).The fact that the non-interacting value of f plays the role of the initial condition for the interacting problem is due to our assumption that particle collisions are less frequent than scattering events on the disorder, Eq. (23).Indeed, in this regime the diffusive ring first establishes, and only then do interactions come into play.
We obtain the diffusive contribution to the momentum distribution for g = 0 by integrating Eq. ( 30) over and taking the Fourier transform with respect to ω.At leading order in the interaction strength, the term proportional to g 2 in the right-hand side can be neglected, such that, in the hydrodynamic regime ωτ 1, we have: This formula quite naturally generalizes the noninteracting Eq. ( 15).In the presence of interactions, the distribution f (t) is no longer constant in time, which leads to an evolution of the diffusive background.Equations (31) and (34) have been used, in particular, in [34], to qualitatively discuss the emergence of a Bose condensate at very long time, but without microscopic justification.We note that the kinetic equation ( 31) has also been derived in [51] by means of a quantum field theory in the presence of disorder.By multiplying the kinetic equation ( 31) by ν and integrating over , we readily obtain ∂ z d ν f (t) = 0.This implies that d ν f (t) = d ν f (0) = 1, which is nothing but the normalization condition (7).It follows that the diffusive contribution (34) is normalized, very much like in the noninteracting limit [50].
The attentive reader will notice that Eq. ( 31) in fact coincides with the free-space Boltzmann kinetic equation for Bose gases in the limit of large occupation numbers [52].Indeed, the kernel (32) is independent of any disorder parameter (in particular, the diffusion coefficient or even the mean free path do not appear).This result is different from the case of electrons in low-temperature disordered conductors, where the diffusive motion strengthens the effect of interactions [53].This difference stems from the mechanism of dynamical screening of electronelectron interactions, which is absent for low-temperature bosons.[54] For weakly interacting diffusive bosons, disorder thus only manifests itself through the spectral function, involved both in the initial condition f (0) and in Eq. (34).The situation would of course change at stronger disorder or in the localization regime [18].

V. CBS OF INTERACTING PARTICLES: THEORY A. Leading-order nonlinear corrections
We now come to the central part of our work and examine the effect of interactions on the series of time-reversed paths, responsible for coherent backscattering.As for the diffusive background, we first address the first-order nonlinear corrections to the Bethe-Salpeter equation of As for diagram (a), the collision processes involve both the direct and the reversed paths.This imposes them to occur almost simultaneously, at t coll t/2, which is very unlikely.This diagram is therefore negligible.Fig. 3b.These corrections are displayed in Fig. 5b and  5c.The diagram 5b has the very same property as its incoherent counterpart 5a: it can be recast as an energy shift −2gρ 0 of the linear Green's function, and thus does not play any role in the dynamics.The building block 5c, on the other hand, turns out to cancel with its conjugate counterpart.At this stage, an important comment is in order.In the stationary scenario considered in [35,36,39], it was shown that specific concatenations of the diagram 5c was leading to a dephasing between the reversed amplitudes, which could even change the sign of the coherent backscattering cone.It turns out, however, that in the present dynamical setup these combinations have a negligible weight.To see this, we show in a more visual fashion in Fig. 7a one example of interference sequence between time-reversed paths built from diagram 5c.The peculiarity of this sequence is that both the direct (solid) path and its time-reversed (dashed) partner are involved in the particle collision process.If the direct path undergoes the collision at a certain time t coll , and thus the time-reversed path at time t−t coll , the temporal locality of the collision imposes that t coll = t − t coll , i.e. that t coll = t/2.In other words, the collision must occur at a very specific time (more precisely, within a time window of width τ , centered around t/2).Within such a short time window, and given the condition (23), it is highly unlikely that two (or more) collisions occur.Any concatenation of diagrams 5c can thus be safely neglected here, and an examination of second-order corrections is again required.

B. Second-order corrections
All non-negligible second-order corrections to the Bethe-Salpeter equation for time-reversed sequences are depicted in Fig. 8.Note that, as compared to the dif-FIG.8. Second-order corrections to the Bethe-Salpeter equation for the CBS contribution to the momentum distribution, I C ,ω (k).Symbols have the same meaning as in Fig. 5. Recall that the two wave paths involved in the blue boxes propagate in opposite directions (the incoming and outgoing momenta k0 and k are explicitly indicated for clarity).The numerical prefactors account for the possible combinations of propagation lines connecting to the vertex g.
. fusive corrections in Fig. 6, there are only four different topologies, and not five.These topologies are the same as those of the four last diffusive diagrams in Fig. 6.It turns out, indeed, that the interference counterpart of the upper-right diagram in Fig. 6 is negligible in the long-time limit.This can be understood from Fig. 7b, which shows the interference sequence that would correspond to the upper-right diagram in Fig. 6 in which an amplitude is time-reversed.Because the two collision processes involve both the direct and the reversed paths, they must occur almost simultaneously, at t coll t/2, which is extremely unlikely given the rarity of particle collisions assumed here.This type of diagram is thus negligible.
The Bethe-Salpeter equation for I C ,ω is similar to Eq. ( 30) except for the missing diagram.The latter is responsible for the term ∝ f 1+ 2− f 1 f 2 in Eq. (31).Since this term is now absent, the f (t) present in each of the other terms can be factored out.We thus obtain, for the "coherent" occupation number the kinetic equation Distributions f (t) and f C (t), obtained by solving the kinetic equations ( 31) and ( 36) for k0 = (π/5, 0), gρ0 = 0.002, and γ = 0.0182, with the initial condition (33) evaluated numerically for g = 0.The different evolutions reflect the asymmetry of the kinetic equations: as time grows the distributions f broaden, whereas the f C flatten out.
where the kernel is still given by Eq. (32), and the initial condition is again set by the non-interacting limit, f C (t = 0) = f (0) .The amplitude of the coherent backscattering peak then follows from: The asymmetry between the kinetic equations for f C and f explains the different dynamic evolution of the CBS peak and diffusive background observed in numerical simulations, Fig.
2 [55].This asymmetry is a major difference with the non-interacting regime.We show in Fig. 9 the evolution of f (t) and f C (t) at short time, obtained by solving the kinetic equations ( 31) and ( 36) with the initial condition (33) evaluated numerically for g = 0.
To find the latter, we have computed the spectral function and the density of states numerically as explained in [56,57].The thermalization mechanism is well visible in the upper graph: the distribution f (t) broadens as time grows.A different behavior is observed for f C (t), which does not broaden but rather flattens out.This phenomenon resembles a dephasing process, a point that will be made more quantitative in the next section.Note, in passing, that the energy corresponding to the maximum of the distributions in Fig. 9 lies always slightly below = 0 .This effect, also seen in experiments [58], is mainly due to the real part of the self-energy (see appendix A).

VI. COMPARISON WITH NUMERICAL SIMULATIONS
In order to test our theoretical approach, we now confront the predictions of the last section to numerical simulations of plane-wave propagation based on the GPE (1).For this purpose, we integrate numerically the collision integrals in the kinetic equations (31) and (36).

A. Diffusive ring and CBS peak amplitudes
In Fig. 10, we reproduce the simulation results of Fig. 2 for the heights of the diffusive ring and CBS peak.For g = 0, we fit them with Eqs. ( 34) and ( 37), with f and f C computed from the kinetic equations ( 31) and (36), using gρ 0 as a fit parameter.To describe the times t ≤ τ , we multiply the right-hand side of Eqs. ( 34) and ( 37) by the same short-time corrections as in the linear case [terms within the square brackets in Eqs. ( 18) and (19)].This is a very good approximation in the regime τ NL τ considered here.As seen in Fig. 10, the agreement between theory and simulations is excellent at all times.We note, however, that for the ring height the fitted values of gρ 0 differ by ∼ 10% compared to those chosen in numerical simulations, and for the CBS height by ∼ 30%.One possible reason for this difference might come from the numerical uncertainties in the resolution of the collision integrals.Evaluating the latter indeed turns particularly challenging in two dimensions, where the kernel W exhibits a number of logarithmic singularities over the integration domain.This discrepancy could also stem from higher-order interaction contributions that renormalize the interaction strength g and not taken into account in the present work.Further investigation would however be required to clarify this point.

B. Decay rates
An important information one may extract from the time-dependent evolutions in Fig. 10 are the characteristic time scales τ D,C NL governing the decay of the diffusive ring and of the CBS amplitude.Theoretically, these characteristic times can be obtained from a short-time expansion of the solution of the kinetic equations ( 31) and (36), f D,C (t) f (t = 0) − α D,C ( )t + O(t 2 ), from which we obtain, using Eq. ( 34) and (37):  34) and (37), and dotted curves are the approximate formula (41).The same value g fit 1.1g (resp.g fit 1.3g) was used for all the diffusive (resp.CBS) curves.Solid curves at g = 0 are obtained from Eqs. (18) and (19). where Here we used -see Eq. ( 17) -that d A (k 0 )f (t = 0) = τ /(πν 0 ).By inserting the Taylor expansions for f and f C in the kinetic equations, we find the functions α D,C ( ) numerically and, from Eq. ( 39), infer the sought out time scales by numerical integration over .This leads to where β D 2.27 and β C 7.17 are numerical prefactors which include the adjustment of the interaction strength used for the fits in Fig. 10.The characteristic times governing the decay of the diffusive background and CBS peak are therefore both proportional to the particle collision time, which we previously introduced in Eq. ( 20).This is quite a natural result, but it should be noted that the decay time for CBS is approximately three times smaller than the decay time for the diffusive background.In other words, the CBS peak is way more sensitive to interactions, as clearly seen in Fig. 10.
To confirm these results, we also compare the theoretical prediction (40) to the decay rates extracted from numerical simulations based on the GPE.For this purpose, we compute numerically the CBS and diffusive signals versus time for several values of the disorder amplitude γ.We then extract the slope of these curves within a narrow time window following the curve maxima (located near t = 15τ in Fig. 10).The results are shown in Fig. 11 as a function of the disorder parameter k 0 (lower points are obtained from the decay of the diffusive ring and upper points from the decay of the CBS peak).The theoretical predictions (40) are shown on the same graph, and they match very well the simulations.These results confirm, in particular, the independence of the collision time on the disorder.By computing numerically the slopes for several values of g, we have also verified the g 2 dependence of (τ D,C NL ) −1 .
While we have not been able to find an exact analytical prediction for the whole time decay of the CBS peak, a simple, approximate expression can be inferred from the kinetic equation for f C , Eq. (36).Indeed, since the timedependence of the diffusive background -encoded in the f i functions in the right-hand side of Eq. ( 36) -is rather slow, in first approximation the occupation number f C (t) decays exponentially.This suggests the simple form × 1 − exp(−t/τ ) 1 + t/τ + t 2 /2τ 2 for the CBS peak amplitude (the second line is the shorttime evolution, which, we recall, is not modified by interactions).Although the fundamental mechanism responsible for the decay of the CBS peak is thermalization, this formula shows that it can also be described in terms of an effective dephasing effect, of rate controlled by the particle collision time.Eq. ( 41) is shown in Fig. 10 (dashed curves) on top of the exact solutions of the kinetic equation, and turns out to be a rather good approximation.

VII. CONCLUSION
In this article, we have constructed a microscopic diagrammatic theory describing the out-of-equilibrium evolution of a weakly interacting disordered Bose gas in momentum space.Assuming weak disorder and rare particle collisions, we have derived coupled kinetic equations for the two main physical processes at work in this regime, particle diffusion and coherent backscattering.Our approach has revealed a noticeable asymmetry in the kinetic equations for these two contributions, implying a faster decay of the CBS peak at short time.Although the relevant physics is the one of thermalization of the particle energy distribution, we have shown that the decay of the CBS peak can be approximated by an effective dephasing mechanism, governed by the particle-particle collision time.Natural extensions of this work concern the role of interactions in the localization regime, where the phenomenon of coherent forward scattering shows up in momentum space [32], and the properties of the Kosterlitz-Thouless transition expected in the equilibrium state reached at long time [59].Related open questions include the exploration of the opposite regime of a disorder weaker than interactions -which then mainly drive the dynamics -, the possible existence of nonthermal fixed points in the presence of disorder, or the outof-equilibrium dynamics in the many-body regime.

FIG. 2 .
FIG. 2. Heights of the diffusive ring and CBS peak versus time for increasing values of g.Data for the ring height are obtained by computing the maximum value of |Ψ(k, t)| 2 at (kx, ky) = (0, ±k0), where no CBS peak is present.Data for the CBS peak height are obtained by subtracting the momentum distribution rotated by 90 • from |Ψ(k, t)| 2 .Here k0 = (π/5, 0), γ 0.0182, and data are averaged over about 16000 disorder realizations.The solid smooth black curves are Eqs.(18) and(19), respectively.They both saturate at a value slightly above 100, which is close to the analytical estimate τ /(πν 0 ) 98.9 for a numerically computed mean free time τ 51.84 and on-shell density of state ν 0 0.167.

FIG. 3 .
FIG. 3. (a) Diagrammatic representation of the Bethe-Salpeter equation for ladder diagrams, Eq. (8).Upper solid lines with arrows refer to k|G R +ω/2 |k0 , and lower dashed lines to k0|G A −ω/2 |k .Note that for solid (dashed) lines, arrows coincide with the (opposite of the) direction of propagation.Dotted vertical lines symbolize the correlation function in Eq. (2).(b) Series of diagrams describing the interference between time-reversed paths, responsible for CBS ("crossed diagrams").The second equality follows from time-reversal invariance.(c) In the long-time limit t τ where low scattering orders have a negligible weight, the series of ladder and crossed diagrams coincide.

FIG. 5 .
FIG. 5. First-order diagrammatic corrections to the Bethe-Salpeter equation, involving one particle collision (for each diagram, one has to add the conjugate version).(a): Correction to the Bethe-Salpeter equation for I D ,ω .Red boxes refer to an incoming, ladder-type sequence I D ,ω , where, apart from particle collision processes, the two paths propagate along the same sequence of scatterers in the same direction.(b) and (c): Corrections to the Bethe-Salpeter equation for I C ,ω .Blue boxes refer to time-reversed scattering sequences.Diagrams (a) and (b) boil down to an irrelevant energy shift, while diagram (c) is compensated by its complex conjugate.The solid and dashed lines symbolize G R and G A , respectively, the vertical lines the correlation function in Eq. (2), and the wavy lines a particle collision, see Fig. 4.

FIG. 7 .
FIG. 7. (a) Example of interference sequence generated by the diagram in Fig.5c(for a better visualization we momentarily change the definition of arrows, which here always indicate the direction of propagation).The direct path, starting at t = 0, undergoes a particle collision at t col , and the time reversed-path at t − t coll .Since the collision is local in time, we must have t coll = t/2.(b) Example of interference sequence between time-reversed paths involving two collisions.As for diagram (a), the collision processes involve both the direct and the reversed paths.This imposes them to occur almost simultaneously, at t coll t/2, which is very unlikely.This diagram is therefore negligible.
FIG. 9.Distributions f (t) and f C (t), obtained by solving the kinetic equations (31) and (36) for k0 = (π/5, 0), gρ0 = 0.002, and γ = 0.0182, with the initial condition (33) evaluated numerically for g = 0.The different evolutions reflect the asymmetry of the kinetic equations: as time grows the distributions f broaden, whereas the f C flatten out.

FIG. 10 .
FIG.10.Heights of the diffusive ring and of the CBS peak versus time for increasing values of gρ0.Colored curves are the same simulation data as in Fig.2.Solid black curves are fits to the theory, Eqs.(34) and(37), and dotted curves are the approximate formula(41).The same value g fit 1.1g (resp.g fit 1.3g) was used for all the diffusive (resp.CBS) curves.Solid curves at g = 0 are obtained from Eqs.(18) and(19).

5 FIG. 11 .
FIG. 11.Decay rates (τ D,CNL ) −1 for the diffusive background (lower orange dots) and CBS (upper blue dots) amplitudes, at fixed gρ0 = 0.001.Dots are obtained from numerical simulations of the Gross-Pitaevskii equation, by computing the slope of the simulation curves in Fig.10for several values of γ (error bars originate from the fitting of the slopes by a straight line).From right to left: γ = 0.0036, 0.0056, 0.0081, 0.0121, 0.0182, 0.0256, 0.0324, 0.0506.Data are displayed as a function of the dimensionless disorder parameter k0 .They confirm the theoretical predictions(40), shown as solid lines, and in particular the independence of the decay rates on the disorder strength.