Dynamics of Many-Body Delocalization in the Time-dependent Hartree-Fock Approximation
Paul Pöpperl, Elmer V. H. Doggen, Jonas F. Karcher, Alexander D. Mirlin, Konstantin S. Tikhonov
DDynamics of Many-Body Delocalization in the Time-dependentHartree-Fock Approximation
Paul P¨opperl , Elmer V. H. Doggen, Jonas F. Karcher Institute for Quantum Materials and Technologies, Karlsruhe Institute of Technology, 76021 Karlsruhe, GermanyInstitut f¨ur Theorie der Kondensierten Materie, Karlsruhe Institute of Technology, 76128 Karlsruhe, Germany
Alexander D. Mirlin
Institute for Quantum Materials and Technologies, Karlsruhe Institute of Technology, 76021 Karlsruhe, GermanyInstitut f¨ur Theorie der Kondensierten Materie, Karlsruhe Institute of Technology, 76128 Karlsruhe, GermanyL. D. Landau Institute for Theoretical Physics RAS, 119334 Moscow, RussiaPetersburg Nuclear Physics Institute, 188300 St. Petersburg, Russia
Konstantin S. Tikhonov
Skolkovo Institute of Science and Technology, Moscow, 121205, RussiaCondensed-matter Physics Laboratory, National Research University Higher School of Economics, 101000 Moscow,Russia
Abstract
We explore dynamics of disordered and quasi-periodic interacting lattice models using a self-consistent time-dependent Hartree-Fock (TDHF) approximation, accessing both large systems(up to L =
400 sites) and very long times (up to t = ). We find that, in the t → ∞ limit,the many-body localization (MBL) is always destroyed within the TDHF approximation. Atthe same time, this approximation provides important information on the long-time characterof dynamics in the ergodic side of the MBL transition. Specifically, for one-dimensional (1D)disordered chains, we find slow power-law transport up to the longest times, supporting therare-region (Gri ffi ths) picture. The information on this subdi ff usive dynamics is obtained bythe analysis of three di ff erent observables— temporal decay ∼ t − β of real-space and energy-space imbalances as well as domain wall melting—which all yield consistent results. For two-dimensional (2D) systems, the decay is faster than a power law, in consistency with theoreticalpredictions that β grows as log t for the decay governed by rare regions. At longest times andmoderately strong disorder, β approaches the limiting value β = ff usion.In quasi-periodic (Aubry-Andr´e) 1D systems, where rare regions are absent, we find considerablyfaster decay that reaches the ballistic value β =
1, which provides further support to the Gri ffi thspicture of the slow transport in random systems. Keywords:
Low-dimensional systems, Many-body localization, Hartree-Fock approximation,spinless Fermi-Hubbard model, Gri ffi ths e ff ects [email protected] a r X i v : . [ c ond - m a t . d i s - nn ] J a n . Introduction Many-body localization (MBL) concerns the interplay of disorder and interactions in highlyexcited states of a many-body system. It is understood that a disordered system that would beAnderson localized in the absence of interaction [1] remains localized upon the introduction ofinteractions in a certain range of parameters [2, 3, 4, 5, 6, 7, 8]. The transition between theergodic and MBL phases is called the MBL transition. For one-dimensional (1D) systems withshort-range interactions, the MBL transition is believed to exist in the thermodynamic limit atfixed (system-size-independent) values of disorder and interaction. In other situations (includinghigher-dimensional systems and models with long-range interactions), the MBL transition stillexists but in a more sophisticated sense, requiring appropriate scaling of the disorder with thesystem size.The physics of the MBL transition and of both phases surrounding it is of great interest. Re-markably, not only the MBL phase is highly unusual but also the delocalized phase exhibits dis-tinct and unconventional properties. Specifically, it has been found numerically that the ergodicside of the MBL transition features slow, subdi ff usive transport in a broad range of parameters[9, 10]. The nature of this slow dynamics is still under debate. In particular, it was proposed thatthe slow transport is due to rare regions of anomalously high disorder leading to Gri ffi ths-typephysics.The problem that one deals with in this context—a many-body problem involving both disor-der and interactions—turns out to be notoriously di ffi cult, both for analytical studies and numer-ical simulations. While a number of analytical approaches to the MBL problem have been de-veloped, all of them either involve approximations or assumptions that are di ffi cult to control, ortreat some simplified toy models. In this situation, numerical simulations are of particular impor-tance. However, the computational treatment of the MBL problem is also an extremely challeng-ing problem. An exact numerical study of even the simplest quantum many-body systems—1Dlattice systems with a single, binary degree of freedom per lattice site—requires a computationalcomplexity that scales as 2 N , where N is the number of sites in the lattice. While clever numer-ical techniques have been employed to push this as far as possible [11], this exponential scalingof exact diagonalization (ED) methods is unavoidable, and systems beyond N ≈
24 are inac-cessible by such methods using reasonable computational resources. While for simpler physicalphenomena systems of this size might be already reasonably su ffi cient, they are much too smallfor understanding the large- N physics of high-energy many-body states of strongly disorderedinteracting systems. Thus, alternative computational approaches that give access to long-timedynamics in much larger systems are of vital importance.One approach that has been successfully applied to study the MBL problem for larger sys-tem sizes is based on the use of tensor-network-based algorithms, such as the density matrixrenormalization group [12, 13] and the time-dependent variational principle (TDVP) with ma-trix product states [14, 15, 16, 17, 18]. The TDVP approach has allowed to explore the drift ofthe MBL transition point W c ( L ) with the length L for 1D chains [15]. It was found that, whilethis drift is rather substantial between L ≈
20 (available for ED) and L =
50, it saturates at L ≈
50 — 100, supporting the existence of the MBL transition at finite W c in the limit L → ∞ in1D. Further, the TDVP approach has demonstrated a qualitative di ff erence between 1D and 2Dproblems in this respect: in the 2D case, the results support an increase of W c with the systemsize without bounds, as expected from the avalanche theory [17]. At the same time, the TDVP Preprint submitted to Annals of Physics January 19, 2021 pproach has its limitations. Specifically, for the large systems, it is limited by relatively mod-est times (typically hundreds of inverse hopping matrix elements) in the vicinity of the MBLtransition. For weaker disorder, the limitation is even stronger in view of a fast growth of theentanglement. While these times are of the same order as those probed in many experiments,there are important questions with respect to the behavior at considerably longer times. As anexample, Ref. [16] has provided some evidence of the di ff erence in the long-time dynamics ofrandom and quasi-periodic systems—which is of fundamental importance for understanding therole of rare events. However, it was clear that data for substantially longer times t are needed toverify the significance of the observed trend.It is therefore important to explore further computational approaches to the physics aroundthe MBL transition that may help to access simultaneously large systems (with N ≈
100 sitesand larger) and long times, t > . Such approaches necessarily involve additional approxi-mations, and one has to investigate whether they provide at least qualitative (or, perhaps, semi-quantitative) understanding of the relevant physics. Various approaches to quantum dynamics inthis type of systems have been suggested in recent years [19, 20, 21]. In the present paper, wefocus on the time-dependent Hartree-Fock approximation (TDHF), which has been put forwardas a potentially fruitful method for describing MBL systems in the recent paper [22]. The TDHFis exact in the non-interacting limit and treats the interaction self-consistently.The goal of this work is to explore systematically the dynamics in MBL systems withinthe TDHF approach. We first compare the approach to state-of-the-art “exact” methods suchas ED and the TDVP (for system sizes and times accessible to these methods). While we doobserve clear deviations, we see that the TDHF approach does capture the key ingredient of theproblem—which is in the focus of our work—the slow dynamics in the ergodic phase. We thusproceed and perform TDHF numerical simulations up to very long times, t ∼ , with systemsizes up to 400 sites. To study the slow, subdi ff usive dynamics of many-body delocalizationat long times, we use three di ff erent observables: (i) temporal decay of real-space imbalance,(ii) decay of energy-space imbalance, and (iii) melting of a domain wall. Importantly, all thesemethods yield consistent results for the flowing (time-dependent) power-law exponent β that weuse to characterize the numerical data.One of the central questions that we address is a comparison between the dynamics for ran-dom and quasi-periodic 1D systems. We find qualitatively di ff erent behavior of β ( t ) in thesetwo cases: while β ( t ) remains relatively small and experiences long-time saturation in randomsystems, it crosses over to the ballistic value β = ffi ths (rare-region) mechanism of slow dynamics in random systems; themechanism that is not operative in quasi-periodic systems [23, 24].To elucidate the role of spatial dimensionality for MBL, we consider also two-dimensional(2D) systems. Our results show that the exponent β in this case does not saturate at a subdi ff u-sive value but rather grows as β ( t ) ∼ log t , again in consistency with expectations based on theGri ffi ths mechanism of slow transport.Our results thus show that the TDHF approach characterizes remarkably well the ergodic sideof the MBL transition. A natural and important question to ask is whether the MBL phase (andthe MBL transition) are also captured by this approximation. This question is also addressed inthe present work. We find (at variance with a proposal in Ref. [22]) that the MBL phase is alwaysdestroyed, in the t → ∞ limit within the TDHF approximation. The corresponding numericalresults are supported by analytical arguments.While this manuscript was in preparation, we learned of a related work by Nandy et al. [25],where the TDHF approximation is used, in combination with numerically exact solution, for the3nalysis of MBL in relatively small 1D systems. Our findings about the absence of a true MBLtransition in the TDHF approximation in 1D is consistent with the results of Ref. [25].
2. Model and Method
We consider a lattice model of interacting spinless fermions described by the Hamiltonian H FH = N (cid:88) i = h i n i − J N (cid:88) i , j = δ (cid:104) i , j (cid:105) c † i c j + U N (cid:88) i , j = δ (cid:104) i , j (cid:105) n i n j , (1) n i = c † i c i , (2)where the c † i and c i are fermionic creation and annihilation operators in site space and N isthe number of sites. Further, δ (cid:104) i , j (cid:105) is equal to unity if the sites i , j are nearest neighbors on theconsidered lattice and zero otherwise. We consider two di ff erent types of on-site fields h i , i ∈ [1 , N ]. The first one is random disorder : h i are taken as random energies uniformly distributedin the interval [ − W , W ], yielding an interacting Anderson model. The second case is a quasi-periodic field, leading to an interacting Aubry-Andr´e model: h i = W π Φ i + φ ) , (3) φ ∈ [0 , π ) . (4)Here Φ is an irrational number rendering the period of the potential incommensurate with thelattice; in this paper we choose Φ = ( √ − /
2. Further, φ ∈ [0 , π ) is a constant phase shiftthat is taken as a random number over which the averaging is performed. In one dimension, themodel (1) maps by the Jordan-Wigner transformation to a spin-1 / W of the quasi-periodic potential as defined inEq. (4). We employ the time-dependent Hartree-Fock (TDHF) approximation to obtain an equationof motion for the lesser Green’s function G < i , j ( t , t (cid:48) ) = i Tr (cid:104) ρ c † j ( t (cid:48) ) c i ( t ) (cid:105) , (5)where ρ is the density operator corresponding to the initial state and the time-dependent opera-tors are in the Heisenberg picture. Given a Hamiltonian of the form H = H + N (cid:88) i , j = V i , j n i n j , (6) H = N (cid:88) i = ε i n i + N (cid:88) i , j = J i , j c † i c j , (7)4he TDHF equation of motion for G < i , j ( t , t (cid:48) ) readsi ∂ t ˆ G < ( t , t (cid:48) ) = (cid:104) ˆ H − ˆ Σ HF ( t ) (cid:105) ∗ ˆ G < ( t , t (cid:48) ) , (8)where Σ HF i , j is the Hartree-Fock self-energy, Σ HF i , j ( t ) = − i δ i , j (cid:88) k V i , k G < k , k ( t , t ) + i V i , j G < i , j ( t , t ) . (9)In Eq. (8) and below we denote matrices in site space with a hat and the corresponding matrixproduct with a star ( ∗ ). Equation (8) is a self-consistent approximation for description of the dy-namics of ˆ G < ( t , t (cid:48) ) in an interacting system [26]. In Ref. [22], this approximation was introducedin the context of MBL.Observables considered in this paper are expressed in terms of the density expectation valuesat individual sites j at time t , which are related to the lesser Green’s function as (cid:68) n j ( t ) (cid:69) = − i G < j , j ( t , t ) : = − i g j , j ( t ) . (10)Therefore it su ffi ces for our purposes to consider the lesser Green’s function with same timearguments, abbreviated by g i , j ( t ) in the last line. According to Eq. (8), the time evolution of g i , j ( t ) is determined byi ∂ t ˆ g ( t ) = (cid:16) ˆ H − ˆ Σ HF ( t ) (cid:17) ∗ ˆ g ( t ) − ˆ g ( t ) ∗ (cid:16) ˆ H − ˆ Σ HF ( t ) (cid:17) : = (cid:104) ˆ H − ˆ Σ HF ( t ) , ˆ g ( t ) (cid:105) . (11)We integrate this equation for successive time points by using an interface to ODEpack [27]solvers contained in SciPy [28].In the process of solving Eq. (11), we have to compute the commutator on its right hand-side at every time step. For the considered models with short-range hopping and interaction,the number of entries in the hopping matrix ˆ J and the interaction matrix ˆ V scale linearly withthe number of sites N . Thus, the necessary matrix multiplications can be performed with O ( N )operations. For the considered system sizes, the complexity of the solving process dependson the size mainly through the evaluation of Eq. (11). Therefore, the computation time scalesapproximately as O ( N · N time ), where N time is the number of time steps. For the solving process,we generally split a unit physical time into 10 steps. We have checked that the chosen time stepis su ffi ciently small to ensure that observables under consideration depend only weakly on it; seeAppendix A for details.For the purpose of benchmarking, we will compare the TDHF to two di ff erent establishedmethods. Exact results for the time evolution of an arbitrary initial state can be obtained bydirectly applying the time evolution operator and calculating matrix elements. The run time andmemory consumption of this process consequently scale exponentially with the system size N .Using the QuSpin package [29, 30] we can simulate up to N =
22 using modest numericalresources. Furthermore, we compare to results from Refs. [15, 31] obtained using the TDVPwith matrix product states [32]. 5 .3. Observables
For 1D systems, a central observable to be studied is the real-space imbalance as a functionof time, which is defined as I ( t ) = (cid:104) n even ( t ) (cid:105) − (cid:104) n odd ( t ) (cid:105) N . (12) n even = (cid:88) i even n i , n odd = (cid:88) i odd n i . (13)Angular brackets here and in the following are understood as denoting averaging over states aswell as over disorder configurations. In the quasi-periodic case, the disorder averaging is replacedby averaging over the phase shift φ , Eq. (4). Initially, we prepare the system in a staggered statewhere all even sites are occupied and all odd states are empty, corresponding to the maximalpossible imbalance, I ( t ) = g ( t ) and will be introduced in respectivesections of the paper below. Errors are estimated using a bootstrapping procedure. Decay of the imbalance.
Previous work suggested that the decay of the imbalance in disordered1D systems is of power-law character [9], I ( t ) ∼ t − β . (14)To characterize the decay found in the simulations, it is convenient to define a flowing (time-dependent) exponent β ( t ) via β ( t ) = − ∂ log( t ) log[ I ( t )] . (15)We evaluate Eq. (15) numerically by determining the slope from a window of intermediate sizethat permits to average out fast fluctuations, cf. Ref. [31].
3. One-dimensional random system
We consider the dynamics of a one dimensional randomly disordered system as described byHamiltonian (1), where we set J = . U = . We compare results obtained with the TDHF to results of a brute-force exact calculation upto 10 hopping times in a chain of L =
22 sites as well as to TDVP results from Ref. [15] ina system of L =
100 sites up to 10 hopping times. Our choice of units corresponds to theconvention from Ref. [15]. In accordance with Ref. [15] we use open boundary conditions forthe comparison. Disorder strengths are between W = W = .1.1. Small systems First, we consider a small system with L =
22 sites. Figure 1 shows the imbalance I ( t ) asobtained from exact calculation and TDHF for W = W =
4. Both the exact and the TDHFcurves show slow dynamics but the long-time decay within the TDHF approach is clearly faster.Clearly, the TDHF is only an approximate method. To inspect the character of the dynamics, weshow in Fig. 2 the time dependence of the slope β ( t ) defined by Eq. (15) for the disorder inter-val W ∈ [2 , β below 0.5 mean slow, subdi ff usive transport. While a quantitativedi ff erence between the exact results and those of the TDHF is again manifest, we see that withinthe TDHF the transport remains clearly subdi ff usive, with β ( t ) < .
3. This is a first indicationof the fact that the TDHF captures the subdi ff usive character of the dynamics, even though it isobtained for rather small systems. We will see below that this remains true for large L and, more-over, that the performance of the TDHF is even better in larger systems. (The latter observationis consistent with Ref. [22] where it was pointed out that field-theoretical approaches, like self-consistent TDHF approximation, tend to e ff ectively reduce finite-size e ff ects, thus mimickinglarger systems.) time t i m b a l a n ce I ( t ) W = 2 exactTDHF 10 time tW = 4 Figure 1: Imbalance as a function of time obtained from exact calculation (blue solid line) and TDHF (orange dashedline) in a 1D random system of L =
22 sites with open boundary conditions (OBC). Exact and TDHF results wereaveraged over ≈
60 and ≈
200 samples respectively. The left and right panel show the imbalance at W = W = time t s l op e ( t ) time t W = 2 W = 3 W = 4 W = 5 W = 6 W = 7 Figure 2: Time dependence of the imbalance exponent (slope) for the disorder W =
2, 3, 4, 5, 6, and 7 in a 1D randomsystem of L =
22 sites with OBC. Left and right panels present exact and TDHF results, respectively. The shaded regionsindicate the bootstrap errors on the individual fits. .1.2. Large systems We proceed now to much larger systems with L =
100 sites. Clearly, brute-force numericsis not possible any more in this case. However, essentially exact results for not too long times, t (cid:46) , can be obtained by using the TDVP approach with matrix product states, Ref. [15].Figure 3 shows a comparison between imbalances at for W = W = ff er noticeably in both plots, the slopesare remarkably close such that a di ff erence is barely visible. To compare the slopes (i.e., theexponents β ) obtained by the two methods more accurately, we plot them (calculated from datain the time window t ∈ [50 , W = W =
8. We see that the dependence β ( W ) given by TDHF is qualitatively verysimilar and numerically quite close to the “quasi-exact” one (obtained by TDVP). At the sametime, it is seen that TDHF somewhat overestimates β (i.e., the dynamics is “more delocalized”within TDHF than it actually is). In particular, β given by TDHF remains positive within errorbars well above the critical disorder strength W c ≈ . time t i m b a l a n ce I ( t ) W = 2 time t W = 4 TDVPTDHF2 3 4 5 6 7 8disorder strength W s l op e ( W ) TDVPTDHF
Figure 3: Comparison of imbalances as obtained from TDVP (blue) and TDHF (orange) in a 1D random system of L =
100 sites with OBC for disorder W = W = ∼ and ∼ disorder realizations, respectively. In the bottom figure, exponents β obtained from fits in the time interval t ∈ [5 · , ] are shown. The error bars (one sigma) are obtained by abootstrapping procedure. Due to the computational e ffi ciency of the TDHF, we can extend our results for large systems( L =
100 sites) up to much longer times. The upper left panel of Fig. 4 shows the time dependent8DHF imbalance I ( t ) up to time t = for disorder from W = W =
8. The figures indicatesthat at these very long times the curves have straight-line asymptotics (on the log-log scale),which corresponds to a power-law decay of the imbalance. In order to reveal the long-timebehavior of I ( t ) in the clearest form, we plot in the upper right panel the exponent β ( t ) (i.e., theslope of the curves from the left panel) in the time interval from t = till 10 . It is seen thatfor the relatively weak disorder, W =
2, the exponent β ( t ) is essentially constant in the wholerange of times, implying a clear power-law behavior. For stronger disorder, W =
3, 4, 5, and6, the running exponent first increases with time but eventually saturates. The saturation timeincreases with disorder strength, so that for strongest disorder in this plot ( W = β ( t ) as obtained by TDHF for relatively short times ( t ∈ [50 , t ∈ [5 · , ]) in a still broader range of disorder (from W = W = time t i m b a l a n ce I ( t ) time t s l op e ( t ) W = 2 W = 3 W = 4 W = 5 W = 6 W = 7 W = 8 W s l op e ( W ) fit in t [50, 10 ] fit in t [5 10 , 10 ] Figure 4: Imbalance I and exponent (slope) β up to 10 hopping times in a 1D random system of L =
100 sites atdi ff erent values of disorder, computed using the TDHF with OBC. The data is averaged over ∼ disorder realizations. Upper left: time dependence of the imbalance I ( t ), for disorder strength from W = W =
8, with the shaded regionsindicating the standard deviation within the disorder sample.
Upper right: time dependence of the running imbalanceexponent β ( t ) (slope of the curves in left panel), with the shaded regions indicating bootstrap errors on the individual fits. Bottom:
Exponents β obtained by fits in two time windows, [50 , · , ], as functions of disorder strengthsin the range from W = W =
14, with bootstrap errors (one sigma error bars).
Two main qualitative conclusions from Fig. 4 are as follows. First, as we have just pointedout, the observed long-time saturation of β ( t ) implies a power-law asymptotic behavior of theimbalance decay, I ( t ) ∼ t − β . Importantly, even for the weakest considered disorder, W = β ≈ .
25, i.e., well below the9i ff usive value 0.5 [34]. Therefore, the TDHF approach reveals anomalous di ff usion in a broadinterval of disorder on the ergodic side of the MBL transition up to the very long times t = .This is in agreement with the result of Ref. [22] (where times up to t ∼ were considered).Secondly, we do not observe a sharp MBL transition within the TDHF approximation, atvariance with the conclusion of Ref. [22]. Specifically, even for the strongest disorder that wehave considered, while the imbalance seems to stay constant at relatively short times ∼ , itstarts to drop then, with β ( t ) becoming distinctly di ff erent from zero. Therefore, our numericalobservations indicate that the TDHF approximation destabilizes the MBL phase, producing veryslow, strongly subdi ff usive dynamics for those values of disorder where the dynamics should becompletely frozen due to MBL. We will return to this question below and explain why this isalso expected from the analytical point of view. As a complementary probe of the dynamics of delocalization, we consider the broadeningprocess of a domain wall initially situated in the middle of a system of L =
50 sites. To this endwe calculate the first moment of the particle density as a function of time, x ( t ) = L (cid:88) i = i (cid:104) n i ( t ) − n i ( t = (cid:105) . (16)In the initial state, the sites 1 ≤ i ≤ L / x (0) =
0. As x ( t ) scales as the square of the domain-wall width, we introducethe running exponent β dw ( t ) characterizing the domain-wall broadening according to β dw ( t ) = ∂ log( t ) x ( t ) . (17)With this definition, one can directly compare the exponent β dw ( t ) to the imbalance exponent β ( t )as both of them describe a scaling of a length scale with time [9].The results for x ( t ) and β dw ( t ) are shown in Fig. 5 for disorder strengths in the interval from W = W =
8. The system size in this plot is L =
50, so that homogeneous occupation of thewhole system would correspond to x ( t ) → L / ≈ W = t = , the moment x ( t ) reaches the value ≈
50. For other disorder strengths, W ≥ x ( t ) (cid:46)
10 in the whole time range. We thus may expect some (relatively weak)finite-size e ff ects at longest times for W = ff ects for W ≥
3. Thisis indeed what is observed: a small decrease of β dw ( t ) for W = t (cid:38) can be presumablyattributed to finite-size e ff ects.In general, the behavior of the domain wall exponent β dw ( t ) is in a good agreement with theimbalance exponent β ( t ) from Fig. 4 although statistical fluctuations in β dw ( t ) are stronger. Thesaturation values of the exponent β dw ( t ) reached for the disorder W =
2, 3, 4 are in the range ≈ . β ( t ). The data on the domain-wallbroadening thus support the conclusion on subdi ff usive asymptotic behavior obtained from theanalysis of the imbalance. 10 time t f i r s t m o m e n t x ( t ) time t s l op e d w ( t ) W = 2 W = 3 W = 4 W = 5 W = 6 W = 7 W = 8 Figure 5: Melting of a domain wall in a 1D random system of L =
50 sites with OBC. The data is averaged over ∼ x ( t ) of the particle density as a function of time, Eq. (16),and the right panel the corresponding exponent β dw ( t ) defined by Eq. (17). The imbalance that was analyzed above was defined as a measure of an “anti ff erromagneticorder” in the real space. Now we analyze the anomalous dynamics by using an analogous ob-servable defined in energy space. To this end, we diagonalize the non-interacting part of theHamiltonian for a given disorder realization and prepare the initial state where every second ofthe resulting single-particle eigenstates (sorted by energy) is occupied. Our initial state is thus anexact eigenstate of the non-interacting Hamiltonian. We time-evolve these states and calculatethe imbalance between occupation of even and odd states within the above energy ordering.In the left panel of Fig. 6, we show the time dependence of the energy-space imbalance fordisorder from W = W =
8. The corresponding exponent defined via Eq. (15) (e.g., the slopeof the curves in the left panel) is displayed in the right panel. Comparison of both panels ofFig. 6 with the corresponding (i.e., upper) panels of Fig. 4 shows a clear similarity between thedynamics of the real-space and energy-space imbalance. The saturation values are rather closefor both types of imbalances, which makes it plausible that the long-time t − β asymptotics forboth of them is characterized by the same value of the exponent β .Our results demonstrate that the anomalous slow dynamics is rather generic and holds forvery di ff erent types of initial conditions. In this context, we mention that it was proposed inRef. [22] that an exact eigenstate of the non-interacting Hamiltonian would not experience ther-malization within the TDHF approximation. Our findings indicate that this is not the case: sucha state thermalizes in a way rather similar to the thermalization of an initial state with “antifer-romagnetic” site occupation as considered in Sec. 3.2. In fact, since each single-particle stateis localized around a certain site, and since nearby-in-energy single-particle states are far awayin real space, there is some qualitative similarity between our initial state with maximal energy-space imbalance and initial states with random site occupation investigated in Ref. [22]. Indeed,the corresponding imbalance traces I ( t ) look quite similar.11 time t i m b a l a n ce I ( t ) time t s l op e ( t ) W = 2 W = 3 W = 4 W = 5 W = 6 W = 7 W = 8 Figure 6: Energy-space imbalance dynamics in a 1D random system of L =
100 sites with PBC, averaged over ∼
100 disorder realizations. The initial state is obtained by occupying every second single-particle eigenstate (in energyordering) of the non-interacting problem. The imbalance is defined between the initially occupied and unoccupied states.Left and right panels show the time dependence of the imbalance I ( t ) and of the exponent (slope) β ( t ) [defined accordingto Eq. (15)], respectively.
4. One-dimensional quasi-periodic system
In Sec. 3 we have studied the many-body delocalization dynamics in 1D random systemsby self-consistent TDHF approach. The results obtained from the analysis of the real-spaceimbalance were supported by the investigation of the energy-space imbalance and domain-wallbroadening. For all of them, we have numerically determined the running exponent β ( t ). All themethods give consistently strong evidence of the slow, subdi ff usive transport with β < . W ≥ t ∼ studied in our work). A plausible explanation ofthe subdi ff usive dynamics is based on Gri ffi ths physics associated with rare events.In order to shed more light on these results, it is important to compare them with thoseobtained for other types of systems, where the rare-event physics is expected to have di ff erentmanifestations or not to be operative at all. This is done in the present section, where we consider1D systems with quasi-periodic potentials, and in Sec. (5), where 2D random systems are studied.In this section we consider the real-space imbalance I ( t ) in 1D chains described by the inter-acting Aubry-Andr´e model. This model was defined in Sec. 2.1; our choice of parameters hereis J = . U = . Φ = ( √ − /
2. The system length was L =
50 for most of thesesimulations; we also studied L =
100 systems to check the role of finite-size e ff ects, as specifiedbelow. We first compare the TDHF results at relatively short times to those obtained by TDVPand then proceed with the analysis of dynamics at much longer times (inaccessible to TDVP). In Fig. 7, we compare the exponent β of the imbalance I ( t ) obtained from TDHF simulationsin the time interval t ∈ [50 , W = W =
6, includes the MBL transition point W c = . ± . t ∼ β when determined on such time scales. In particular, as a result of such12scillations, β provided by TDVP increases a little from W = W = . W =
5. The TDHF approximation somewhat smears these oscillations. W s l op e ( W ) TDVPTDHF
Figure 7: Comparison of imbalance exponents β in the time interval t ∈ [50 , L =
50 sites with OBC. The TDVP and TDHF imbalances are averaged over ≈
400 and ≈ β =
0. Error bars are two sigma, as in Ref. [31]. time t i m b a l a n ce I ( t ) time t s l op e ( t ) W = 4.0 W = 4.5 W = 5.0 W = 5.5 W = 6.0 time t i m b a l a n ce I ( t ) time t s l op e ( t ) L = 50 L = 100 Figure 8: Imbalance dynamics in a quasi-periodic system with OBC. The data is averaged over ≈
300 disorder real-izations.
Top panels:
Imbalance I ( t ) (left) and the running exponent β ( t ) (right) for L =
50 and di ff erent strengths ofdisorder from W = W = Bottom panels:
Imbalance I ( t ) (left) and the running exponent β ( t ) (right) at disorderstrength W = L =
50 and L =
100 . Dashed lines in the left panels indicate the slope β = .2. Long-time imbalance dynamics In analogy with the model with true randomness, we extend the TDHF analysis up to t = hopping times. The results for disorder between W = W = I ( t ) and the right panel the corresponding runningexponent (slope) β ( t ). As for the random system, we observe an initial increase of the slope withtime. However, there is a crucial di ff erence: for the quasi-periodic disorder β ( t ) increases muchfaster and saturates at a value β ≈
1. After this saturation, β ( t ) decays again down to zero. Thisdecay is related to finite size of the system. To demonstrate this, we compare in the lower panelsof Fig. 8 the data for the same disorder strength, W =
4, and two di ff erent system sizes, L = L = β ≈
1, so that the decay of β ( t ) is shifted to still later times. This confirms the plateau at β = L and the finite-size character of the decay of β ( t ) after this plateau. Theseresults for quasi-periodic systems confirm that the slow, subdi ff usive transport (with β < . ffi ths e ff ects.The strong di ff erence between the dynamics in random and quasi-periodic 1D systems is alsoillustrated in Fig. 9 where we show the real part of the Fourier transform ˜ I ( ω ) of the imbalance I ( t ). The presented data correspond to the weakest disorder strength that we have studied: W = W = I ( ω ) is fitted very well by a straight line on the log-log scale, which means a power-law dependence ˜ I ( ω ) ∼ ω β − , with β ≈ .
25. This value of β is in full correspondence with theone obtained from the analysis of the running exponent in time representation, see upper rightpanel of Fig. 4. At the same time, for the quasi-periodic case, the ˜ I ( ω ) dependence shows a strongcurvature on the log-log scale and becomes nearly flat for the smallest ω , in correspondence withthe 1 / t long-time behavior of I ( t ). I () I () Figure 9: The real part of the Fourier transform ˜ I ( ω ) of the imbalance I ( t ) for random (upper panel; W =
2) and quasi-periodic (lower panel; W =
4) 1D systems of size L = − . I ( ω ) ∼ ω − . ] is shown. β = q but rather at q = π/ a , where a is the lattice constant. However, in analogy with thecase of random systems (see the comment [34]), we expect that mode coupling will induce a 1 / t decay of the imbalance.
5. Two-dimensional random system
In this section, we investigate the long-time delocalization dynamics in 2D systems deter-mined by the Hamiltonian (1) with random potential on square lattices. We set the parameters J = U = ff erent values of disorder strength W . The boundary conditionsare chosen periodic in one direction and open in the other direction, as in Ref. [35], so that thesystem can be viewed as surface of a cylinder. Following Ref. [35], we use as an observable the“columnar imbalance” between rings of sites. In the initial state, every second ring is occupied,so that the columnar imbalance has its maximal value (unity). time t i m b a l a n ce I ( t ) W = 5 W = 7 W = 10 W = 12 W = 15 W = 20 W = 25 W = 30 W = 35 W = 60 time t s l op e ( t ) time t i m b a l a n ce I ( t ) W = 5 W = 7 W = 10 time t s l op e ( t ) Figure 10: Dynamics of columnar imbalance in a 2D random system. The boundary conditions are periodic along onedirection and open along the other.
Upper panels: system size 10 ×
10 sites; disorder from W = W = Lowerpanels: system size 20 ×
20 sites; disorder W = W =
7, and W =
10. Left panels show the imbalance I ( t ), and rightpanels the corresponding exponent β ( t ). Averaging is performed over ≈
300 samples. I ( t ) and for the corresponding running ex-ponent β ( t ) are shown in upper panels of Fig. 10 for systems of 10 ×
10 sites in the disorder rangefrom W = W =
60. The exponent β ( t ) increases with time, reaching values considerablylarger than in the 1D random case. In particular, for disorder strengths W = W = β ≈ . β ≈ .
7, respectively. After reaching these maximal values, theexponent β ( t ) again decreases. As we have already discussed in Sec. 4, this decay is an evidenceof finite-size e ff ects. To explicitly demonstrate this, we show in the low panels of Fig. 10 thedata for disorder W = W =
7, and W =
10 in larger systems—of size 20 ×
20 sites. We seethat now β ( t ) keeps increasing towards β =
1. We argue that in the limit of large L and long t there will be an extended plateau at the 2D di ff usive value β = β =
1, one can ask a question aboutan increase of β ( t ) before this asymptotics is reached. Assuming that this increase is determinedby rare events, i.e., by generalizing the arguments that lead to a power-law behavior in the 1Dcase, one finds [36, 35] f ( t ) ∼ t − γ ln t = e − γ ln t . (18)for the columnar imbalance in two dimensional systems. As shown in Fig. 11, our data for β ( t )are fitted very well by Eq. (18). time t i m b a l a n ce I ( t ) W = 5 W = 7 W = 10 W = 12 W = 15 W = 20 W = 25 W = 30 W = 35 W = 40 Figure 11: Fitted columnar imbalance in a 2D random system of 10 ×
10 sites (the data from the upper left panel ofFig. 10). The figure shows the columnar imbalance (solid, in color) and fits to Eq. (18) (dashed, black). All imbalancecurves are fitted in the time interval [5 , ], with the exception of W = W = , ] toexclude strong finite-size e ff ects at t (cid:38) .
6. TDHF and MBL: analytical considerations
In our numerical analysis of 1D random systems in Sec. 3, we have pointed out that, withinthe TDHF approximation, delocalization takes place at any disorder strength. This is seen, e.g., inFig. 4: the exponent β ( t ) increases with time up to the longest times for disorder strengths that arewell above the actual MBL transition. This means that the TDHF approximation produces finitedephasing, leading to finite quasiparticle life time, relaxation, and therefore delocalization, evenfor strong disorder, where exact solution would yield MBL (i.e., no dephasing and relaxation).Therefore, there should be delicate cancellations in the exact equations of motion for the Green’sfunctions that are not fully respected by the TDHF approximation. In this section, we provide16he corresponding analytical considerations based on an expansion of the equations of motion inpowers of interaction. First we show that, within the self-consistent Hartree-Fock approximation,such dephasing terms indeed emerge in the second order in interaction V . Second, we shed lighton the mechanism of cancelation between the Hartree-Fock terms and those arising from thesecond-order self-energy in the exact solution in MBL phase. Technical details of the analysis inthis section are relegated to Appendix B.In the TDHF approximation, the self-energy Σ depends linearly on the Green’s function: Σ HF αβ ( t ) = V M αβγδ G γδ ( t , t ) . (19)Greek indices here and in the following are used to denote components in the eigenbasis of thenoninteracting Hamiltonian, distinguishing them from position-basis components for which weuse latin indices as in Sec. 2. We separate the interaction strength V from the remaining part ofthe interaction matrix element M αβγδ in order to make the perturbative expansion in powers of V more clear: G = G (0) + VG (1) + V G (2) + O ( V ) , Σ = V Σ (1) + V Σ (2) + O ( V ) . (20)We tackle the self-consistent equation by solving for G ( i ) and then computing Σ ( i + using Eq.(19)order by order in V , with the time evolution of all expressions starting at t =
0. We refer the readerto Appendix B for details of this perturbative solution scheme.Having obtained the Green’s function from Eq. (20), we then compute its Wigner-transform˜ G <αα ( t , (cid:15) ) up to second order in V . Most of resulting terms only renormalize the quasiparticleweight or energy. However, there is also a second-order contribution [to be denoted ˜ G <, (2) b αα ( t , (cid:15) )]that introduces broadening of the level ω α :˜ G <, (2) b αα ( t , (cid:15) ) = π i (cid:88) γµν M αγνν M αγµµ ρ µ ρ ν ( ω α − ω γ ) × (cid:20) ω α − ω γ ) t ] δ (cid:18) (cid:15) − ω α + ω γ (cid:19) ( ρ γ − ρ α ) − δ ( (cid:15) − ω γ ) ρ γ (cid:21) , (21)where ρ α are initial populations of levels α . This shows that the δ -functions δ ( (cid:15) − ω α ) corre-sponding to the non-interacting levels ω α are broadened. This is the source of the Hartree-Fockmany-body delocalization observed in the simulations. Since the interaction is short-ranged, anessential contribution to the sum in Eq. (21) is given by states γ, µ, ν that are located close (withinthe distance of order of single-particle localization length) to the state α .At weak interaction, the 1D random system is in the MBL phase, so that there should be nodephasing if the problem is solved exactly. Thus, the TDHF broadening (21) should be compen-sated by another contribution. Since such a compensation is only possible between terms of thesame order in interaction V , we should inspect other terms (not included in TDHF approxima-tion) in the Green function that are of the second order in V .To this end, we include the second order self-energy (discarded by TDHF approximation) Σ ( D ) ,< i j ( t , t (cid:48) ) = (cid:88) k , l V il V jk G > kl ( t (cid:48) , t ) (cid:104) G < lk ( t , t (cid:48) ) G < i j ( t , t (cid:48) ) − G < l j ( t , t (cid:48) ) G < ik ( t , t (cid:48) ) (cid:105) (22)into the perturbative scheme where we used standard notations, with the superscripts > and < denoting greater and lesser Green’s functions. Substituting the second in the equation of motion17or the Green’s function, we find the following second-order contribution to the Green’s function(see Appendix B for detail):˜ G <, (2 (cid:48) ) b αα ( t , (cid:15) ) = π i (cid:88) µνγ M ναµγ M αµγν ( ω α − ω γ + ω ν + ω µ ) (cid:110) ( ρ γ − ρ ν ρ µ δ ( ω γ − ω ν − ω µ − (cid:15) ) − [ ρ α ρ γ ( ρ ν + ρ µ − − ρ α ρ ν ρ µ + ( ρ γ − ρ ν ρ µ ] e − i t ( ω α − ω γ + ω ν + ω µ ) δ (cid:18) (cid:15) − ω α + ω γ − ω ν − ω µ (cid:19)(cid:27) . (23)The superscript “(2 (cid:48) )” in the l.h.s. of Eq. (23) serves to distinguish it from the TDHF contribution(21).The contributions (23) and (21) have rather similar general structure. Moreover, separatingin Eq. (23) the contribution of matrix elements with only three distinct labels, like those in Eq.(21), one gets an expression whose form is almost identical to Eq. (21), see Appendix B. Thishelps to understand how a cancelation required to restore the MBL may work. Only when takentogether, the TDHF and non-TDHF terms describe fully the quantum-coherent dynamics, andthus localization at strong disorder. Discarding some of the terms leads to decoherence and, asa result, to delocalization. Of course, such compensation should be operative to all orders in V to ensure the MBL. We thus expect that if the self-energy in the equation of motion is calculatednot to the order V (as in the TDHF) but to the order V (or any finite order V k ), the system willstill experience delocalization for arbitrarily strong disorder at su ffi ciently long times.Finally, we note that, since the TDHF approximation destabilizes localization by inducing afinite exponent β of anomalous di ff usion in the MBL phase, W > W c , it leads, by continuity, toan enhancement of β also in the ergodic phase, W < W c —at least, in some vicinity of W c . Thisis indeed what is observed in the lower panel of Fig. 3 (where W c ≈ . β appears tobe rather good for moderately strong disorder. For example, for W =
2, the exact (TDVP) resultfor β ( t ) as determined in the time range [50 , β ≈ .
18, while the TDHF yields in the sametime interval β ≈ .
20. It remains to be understood whether there is a regime in which the TDHFframework would provide a parametrically controlled approximation for the quantum dynamicsof many-body delocalization.
7. Summary
In this paper, we have studied the dynamics of many-body delocalization within the TDHFapproximation. Our study included 1D random, 1D quasi-periodic, and 2D random systems.We have accessed large system sizes and very long times t (up to 100 sites and t = for 1Dsystems and up to 20 ×
20 sites and t = for 2D systems). We have analyzed the dynamicsusing di ff erent initial conditions: in addition to the setting with initial charge-density wave that isparticularly popular in experiments and numerical studies, we have also investigated the domainwall melting as well as the dynamics starting from an initial state with alternating occupation ofnon-interacting states in the energy space. Our key results are as follows.1. For all the settings, we have characterized the dynamical process by a running exponent β ( t ). We have demonstrated that this is a very useful characterization of the quantumdynamics that allows one to distinguish sensitively between distinct regimes of many-bodydelocalization dynamics. 18. While the TDHF method is a priori approximate, it does capture various key propertiesof the dynamics in highly excited states of interacting disordered (or quasi-periodic) sys-tems. A comparison with the (essentially exact) results by TDVP on large systems (and atrelatively short times accessible to TDVP) shows that TDHF performs rather well on theergodic side of the MBL transition, yielding values of β ( t ) reasonably close to exact ones,see Fig. 3.3. For random 1D systems (Sec. 3.2, 3.3, and 3.4), we have found slow, subdi ff usive dy-namics of many-body delocalization, with the exponent β ( t ) saturating at long times at adisorder-dependent value β corresponding to power-law relaxation ∼ t − β . In the wholerange of considered disorder strength, the resulting values of this exponent satisfy β < . ff usive value 0.5. We have found consistent values of β for all three types of initial conditions, for which we investigated the time evolution ofthe real-space imbalance, the domain-wall width, and the energy-space imbalance, respec-tively. These results provide clear evidence of slow, subdi ff usive transport on the ergodicside of the MBL transition, thus confirming and extending the corresponding findings ofRef. [22]. The origin of this slow dynamics is attributed to rare-event Gri ffi ths physics,which is further supported by the study of quasi-periodic systems and of 2D systems (seethe next two items in this list).4. For quasi-periodic (Aubry-Andr´e) 1D systems (Sec. 4), we have found that β ( t ) character-izing the imbalance has a very di ff erent behavior: it increases with time up to its saturationat the ballistic value β =
1. This confirms that the subdi ff usive saturation value of β inrandom systems is related to rare localized spots (which are absent in a quasi-periodicsystem).Our findings for the quasi-periodic systems are consistent with previous numerical resultsin Refs. [22, 31] where signatures of a faster-than-power-law decay at intermediate timeswere found. In fact, one of the earliest experimental studies of MBL considered preciselythe interacting Aubry-Andr´e model [37]. Experimentally, a slow decay of the imbalancewas observed. The time range available to the experiment was, however, relatively short;at such times the di ff erence between truly random and quasi-periodic systems is not sopronounced yet.5. For 2D random systems (Sec. 5), the imbalance decay is found to follow, at intermediatetimes, the exp (cid:16) − γ ln t (cid:17) law (i.e., β ( t ) increases as ln t ), as expected within the Gri ffi ths-type picture. At longest times, our results indicate a saturation at β =
1, which is the valueexpected from the memory-e ff ect coupling between the large-wave-vector imbalance modeand the di ff usive 2D mode.Also in the 2D case, experiments find slow transport [38] (rather than di ff usive behavior).Again, this is not surprising, taking into account not so long times that have been probedexperimentally. In such a limited time window (say, t ∼ ), it is di ffi cult to distinguishthe exp (cid:16) − γ ln t (cid:17) decay from a simple power law. This is also the case for TDVP study[35] where the characteristic time scale is comparable to that in the experiment. It is a keyadvantage of the TDHF approach that it provides access to much longer times and thusallows one to see the qualitative di ff erence between 1D and 2D models.6. While the TDHF approximation characterizes successfully the quantum dynamics in thedelocalized phase, there is no true MBL transition in this approximation. Our conclu-sion in this respect is at variance with that of Ref. [22] and is consistent with Ref. [25].We find that β ( t ) characterizing the TDHF data slowly increases with time even for very19trong disorder, implying no localization. We have provided also analytical arguments(using an expansion of equations of motion in interaction strength) that shed light on thecorresponding mechanism (Sec. 6). In particular, we have demonstrated that terms provid-ing dephasing occur, within the TDHF approximation, in the second-order in interaction.These terms are very similar in structure to terms of the same order that are discarded bythe TDHF approximation (those resulting from the second-order self-energy). We arguedthat only when both TDHF and non-TDHF terms are retained, the evolution equations cor-respond to a time-independent Hamiltonian and thus should yield MBL at strong disorder.Keeping only TDHF terms generates dephasing also in the regime of strong disorder wherethe exact solution would give MBL.In conclusion, the TDHF approximation is a very valuable computational tool (complemen-tary to other existing methods) for the analysis of quantum dynamics of disordered and quasi-periodic interacting many-body systems. Crucially, it allows one to explore the wealth of regimesof many-body delocalization dynamics in large systems and, at the same, at very long time scales.
8. Acknowledgments
We acknowledge useful discussions with S. Bera, F. Evers, and M. Knap. We are particu-larly thankful for fruitful discussions with I. Gornyi, especially for his advice on the analyticalapproach. This research was financially supported by the DFG-RFBR Grant [No. MI 658 / Appendix A. Integration step size
Solving the TDHF equations of motion (11) numerically, we have to specify the step size ∆ t for the integration procedure. All of our results were integrated by steps of ∆ t = − . In thisAppendix, we analyze the dependence of the results on the step size. Figure A.12: Imbalance in the random 1D model as a function of time for two di ff erent step sizes, ∆ t = − and ∆ t = . · − , calculated for the same single disorder realization. The system size is L =
50 sites, with open boundaryconditions. The disorder strength is W = W = J = . U = . In Fig. A.12, we show a comparison of imbalances in the random 1D model at W = W =
4, as calculated for a single disorder realization with our standard step size ∆ t = − and20maller steps ∆ t = . · − for a system of L =
50 sites, open boundary conditions, and ourusual choice of hopping and interaction parameters, J = U = . ∆ t = − and ∆ t = . · − , the same realization of disorder is chosen.We observe that the imbalance values obtained for simulations with di ff erent time steps beginto deviate significantly at t ≈ · hopping times. The reason for this are fast fluctuations of I ( t ) for an individual disorder realization superimposed on a smooth decay of I ( t ). Even thougherrors originating from a finite time steps are small, they accumulate by the time t ∼ andstrongly disturb the specific pattern of these fast fluctuations.We stress, however, that our analysis is based upon averaged quantities: we are interestedin the smooth behavior of the imbalance (and other observables), with exact details of fast fluc-tuations in a given realization of disorder being of no importance. In Fig. A.13 we present theimbalance averaged over ≈
500 realizations of disorder; all other parameters are the same as inFig. A.12. We see that, upon averaging, the results are nearly identical for two di ff erent timesteps up to t = , which justifies the analysis of the long-time regime with the chosen step size.Deviations between the traces may serve as a measure of an expected error resulting from a finitestep size. i m b a l a n ce I ( t ) W = 2 time t i m b a l a n ce I ( t ) W = 5 t = 10 t = 0.5 10 Figure A.13: Imbalance in the random 1D model as a function of time for two di ff erent step sizes, ∆ t = − and ∆ t = . · − , averaged over ≈
500 disorder realizations. All other parameters are the same as in Fig. A.12. The systemsize is L =
50 sites, with open boundary conditions. The disorder strength is W = W = J = . U = . ppendix B. Quasiparticle Broadening This Appendix presents technical details behind Sec. 6 on the mechanism of dephasing(quasiparticle broadening) that governs many-body delocalization within the TDHF approxima-tion. In Appendix B.1, we calculate the TDHF Green’s functions and self-energies to the secondorder in the interaction strength V , thus providing the derivation of Eq. (21). At strong disorder,one expects MBL to result from the exact time evolution of the interacting problem. Thus, theterms in the total Green’s function that result from TDHF approximation should combine withthose discarded by the TDHF in such a way that the MBL is restored (i.e., no level broadeningis generated). To understand how this “cancellation” can happen, in Appendix B.2 we derivethe second-order terms in the Green’s function that are beyond the TDHF. The result is given byEq. (23) and shows that TDHF and non-TDHF contributions have indeed a similar structure. Atruncation (like that performed within the TDHF approximation), when only a part of the terms iskept, apparently leads to decoherence and thus delocalization at arbitrarily strong disorder, sincethe quantum dynamics then does not correspond to a time-independent many-body Hamiltonian. Keldysh equations of motion: Notations.
In the numerics described in the main text, we haveconsidered the TDHF equations of motion for the Green’s function (8). As these equations arenonlinear in G , an exact analytical solution within the TDHF scheme turns out to be impossible;however, the locality of the TDHF approximation in time allows for e ffi cient numerical solutions.We start with the exact equations of motion, including memory integrals that are non-local intime:i ∂ t G < ( t , t (cid:48) ) = (cid:104) H + Σ HF ( t ) (cid:105) G < ( t , t (cid:48) ) + (cid:90) t d t (cid:48)(cid:48) Σ R ( t , t (cid:48)(cid:48) ) G < ( t (cid:48)(cid:48) , t (cid:48) ) − (cid:90) t (cid:48) d t (cid:48)(cid:48) Σ < ( t , t (cid:48)(cid:48) ) G A ( t (cid:48)(cid:48) , t (cid:48) ) , (B.1)i ∂ t G A ( t , t (cid:48) ) = (cid:104) H + Σ HF ( t ) (cid:105) G A ( t , t (cid:48) ) − (cid:90) t (cid:48) d t (cid:48)(cid:48) Σ A ( t , t (cid:48)(cid:48) ) G A ( t (cid:48)(cid:48) , t (cid:48) ) . (B.2)Here, we have singled out the Hartree-Fock part of the self-energy, Σ HF i , j ( t ) = − i δ i , j (cid:88) k V i , k G < k , k ( t , t ) + i V i , j G < i , j ( t , t ) . (B.3)We use Latin indices for position space and Greek indices for the eigenbasis of the non-interactingsingle-particle Hamiltonian H . The remaining part of the self-energy Σ is a functional of greaterand lesser Green’s function G > , G < of a generic form. This means that, in general, skeleton dia-grams of any order (with respect to the interaction V ) contribute to Σ , starting with the second-order diagrams: Σ < = Σ ( D ) ,< + Σ ( D ) ,< + . . . . In particular, the second-order non-HF skeletondiagram yields: Σ ( D ) ,< i , j ( t , t (cid:48) ) = (cid:88) k , l V i , l V j , k G > k , l ( t (cid:48) , t ) (cid:104) G < l , k ( t , t (cid:48) ) G < i , j ( t , t (cid:48) ) − G < l , j ( t , t (cid:48) ) G < i , k ( t , t (cid:48) ) (cid:105) . (B.4)The retarded, advanced, and Keldysh Green’s functions are defined as usual: G R ( t , t (cid:48) ) = Θ ( t − t (cid:48) ) (cid:2) G > ( t , t (cid:48) ) − G < ( t , t (cid:48) ) (cid:3) , G A ( t , t (cid:48) ) = − Θ ( t (cid:48) − t ) (cid:2) G > ( t , t (cid:48) ) − G < ( t , t (cid:48) ) (cid:3) , G K ( t , t (cid:48) ) = G > ( t , t (cid:48) ) + G < ( t , t (cid:48) ) . (B.5)22etarded, advanced, and Keldysh self-energies are defined analogously. The Wigner transform˜ G ( t , (cid:15) ) is defined as ˜ G ( t , (cid:15) ) = (cid:90) ∞−∞ d τ e i (cid:15)τ G ( t + τ/ , t − τ/ . (B.6)It is used to determine properties of quasi-particles at time t . The time dependence of G ( t , t ) canbe obtained by integrating ˜ G ( t , (cid:15) ) over all energies (cid:15) .A convenient for the perturbative expansion notation for the interaction matrix elements isobtained by splitting the interaction strength V out of the potential V i , k = V ( δ i , k + + δ i , k − ) thatmediates the density-density interaction. Specifically, we introduce matrix elements M writtenin terms of the exact wavefunctions φ α ( i ) of the noninteracting disordered Hamiltonian H : V · M αβγδ ≡ (cid:88) i , k V i , k B ∗ αδ ( i , k ) B γβ ( i , k ) , (B.7) B αβ ( i , j ) ≡ φ α ( i ) φ β ( j ) − φ α ( j ) φ β ( i ) . (B.8)We choose our wavefunctions (Anderson model) to be real. Appendix B.1. Hartree-Fock Green’s functions
In this section, we analyze the Green’s functions in the TDHF approximation, discardingthe non-Hartree-Fock contributions to the self-energy. The Hartree-Fock time evolution can bemapped to the problem of a non-interacting Hamiltonian H supplemented by a time-dependentexternal field S ( t ) with an additional constraint that S ( t ) = Σ HF [ G ]( t ). In the interaction picture,we write:i ∂ t G I ( t , t (cid:48) ) = Σ I ( t ) G I ( t , t (cid:48) ) , G I ( t , t (cid:48) ) = exp(i H t ) G ( t , t (cid:48) ) exp (cid:0) − i H t (cid:48) (cid:1) , Σ I ( t ) = exp(i H t ) S ( t ) exp( − i H t ) , G I (0 , = G ( ρ ) . (B.9)Here G ( ρ ) encodes the dependence of the initial condition on the initial density matrix ρ . Thesolution is given by a formal power series in terms of the time-dependent field S ( t ): G I ( t , t (cid:48) ) = U ( t , G ( ρ ) U (0 , t (cid:48) ) , U ( t , t (cid:48) ) = T exp (cid:32) − i (cid:90) tt (cid:48) dt (cid:48)(cid:48) Σ I ( t (cid:48)(cid:48) ) (cid:33) .. (B.10)If we choose the external field S ( t ) to be the Hartree-Fock self-energy Σ HF defined in (B.3), weget a complicated self-consistent equation. We then perform a perturbative expansion of G and Σ HF in powers of V : G = G (0) + VG (1) + V G (2) + O ( V ) , Σ HF = V Σ (1) + V Σ (2) + O ( V ) . (B.11)We know G (0) from the free problem (time evolution with H ). With this knowledge, we cancompute Σ (1) using Eq. (B.3), which in turn enables us to find G (1) via Eq.(B.10). Repeatingthese steps we can go, in principle, to any order in V . We carry out now this procedure explicitlyup to the second order. 23t zeroth order in V , we have the Green’s functions of a non-interacting Hamiltonian H (with eigenenergies ω α ), which we write in the basis of exact single-particle eigenstates:˜ G αβ ( t , (cid:15) ) = (cid:90) ∞−∞ d τ e i( (cid:15) − ( ω α + ω β ) / τ G αβ ( t + τ ) , t ) , ˜ G (0) , R αβ ( t , (cid:15) ) = δ αβ (cid:15) − ( ω α + ω β ) / + i0 , (B.12)˜ G (0) ,<αβ ( t , (cid:15) ) = π i ρ ,αβ δ (cid:18) (cid:15) − ω α + ω β (cid:19) e − i( ω α − ω β ) t , (B.13)˜ G (0) ,>αβ ( t , (cid:15) ) = − π i ( δ αβ − ρ ,αβ ) δ (cid:18) (cid:15) − ω α + ω β (cid:19) e − i( ω α − ω β ) t , (B.14)˜ G (0) , K αβ ( t , (cid:15) ) = − π i ( δ αβ − ρ ,αβ ) δ (cid:18) (cid:15) − ω α + ω β (cid:19) e − i( ω α − ω β ) t . (B.15)The zeroth order Green’s is a Hermitean matrix even if ρ does not commute with H . Whenwe start from a diagonal ρ , e.g., a staggered state in energy space, it is obvious that the Green’sfunction does not explicitly depend on t and the spectrum stays the same.At higher orders in V , for simplicity, we will restrict our consideration to the case when ρ is diagonal, ρ ,αβ = δ αβ ρ α , and hence the self-energy Σ (1) αβ involves the time-independent freeGreen’s function G (0) ( t , t ) and is thus itself time-independent: Σ (1) αβ = (cid:88) γ M αβγγ ρ γ . (B.16)Below we analyze G < with the initial condition G αβ ( ρ ) = G <αβ (0 , = i ρ α δ αβ in equations of motion (B.9) (for the other Green’s functions one has to change the initial condi-tion of the equations of motion). Expanding Eq. (B.10), we get:˜ G (1) ,<αβ ( t , (cid:15) ) = (cid:90) d τ e i (cid:15)τ e − i ω α ( t + τ/ e i ω β ( t − τ/ (cid:34)(cid:90) t + τ/ dt (cid:48) Σ I ,αβ ( t (cid:48) ) ρ β − (cid:90) t − τ/ dt (cid:48) ρ α Σ I ,αβ ( t (cid:48) ) (cid:35) = (cid:90) d τ e i (cid:15)τ e − i ω α ( t + τ/ e i ω β ( t − τ/ Σ (1) αβ (cid:34) e i( ω α − ω β )( t + τ/ − ω α − ω β ) ρ β − e i( ω α − ω β )( t − τ/ − ω α − ω β ) ρ α (cid:35) α (cid:44) β = π i Σ (1) αβ ω α − ω β (cid:20) ρ β δ ( (cid:15) − ω β ) − ρ α δ ( (cid:15) − ω α ) + e − i( ω α − ω β ) t δ (cid:18) (cid:15) − ω α + ω β (cid:19) (cid:16) ρ α − ρ β (cid:17)(cid:21) . (B.17)To read o ff the quasiparticle energies and the lifetimes, the full matrix ˜ G < ( t , (cid:15) ) = (cid:80) n V n ˜ G ( n ) ,< ( t , (cid:15) )needs to be diagonalized (at fixed t , (cid:15) ). This causes additional mixing of the individual ordersin the expansion in V : in particular, the second-order correction to the diagonalized Green’sfunction will contain a contribution from the first-order o ff -diagonal terms from Eq. (B.17). Notethat new levels ( ω α + ω β ) / ∝ V .To the first order in V , the diagonal component of G < can be directly found from the first lineof Eq. (B.17): ˜ G (1) ,<αα ( t , (cid:15) ) = (cid:90) d τ · e i (cid:15)τ e − i ω α τ τ Σ (1) αα ρ α = π i δ (cid:48) ( (cid:15) − ω α ) Σ (1) αα ρ α . (B.18)24his is just the first term in the formal expansion of the zero-order Green’s function (B.13) witha shifted dispersion: − π i δ (cid:16) (cid:15) − ω α − V · Σ (1) αα ρ α (cid:17) = ˜ G (0) ,<αα ( t , (cid:15) ) + V · ˜ G (1) ,<αα ( t , (cid:15) ) + . . . (B.19)This means that the levels ω α are shifted by the t -independent term V Σ (1) αα = V (cid:80) γ M ααγγ ρ γ . Inhigher-order contributions to the self-energy, we will also identify terms with the derivative ofthe exact delta-function, δ (cid:48) ( (cid:15) − ω α ), with the level shifts.At second order in V there are several contributions: G (2) ,< I ( t , t (cid:48) ) = G (2 , ,< I ( t , t (cid:48) ) + G (2 , ,< I ( t , t (cid:48) ) , (B.20) G (2 , ,< I ( t , t (cid:48) ) = (cid:90) t d t (cid:48)(cid:48) Σ I [ G (1) ,< ( t (cid:48)(cid:48) , t (cid:48)(cid:48) )] ρ − ρ (cid:90) t (cid:48) d t (cid:48)(cid:48) Σ I [ G (1) ,< ( t (cid:48)(cid:48) , t (cid:48)(cid:48) )] , G (2 , ,< I ( t , t (cid:48) ) = − i (cid:90) t d t (cid:48)(cid:48) Σ I [ G (0) ,< ( t (cid:48)(cid:48) , t (cid:48)(cid:48) )] (cid:90) t (cid:48)(cid:48) d t (cid:48)(cid:48)(cid:48) Σ I [ G (0) ,< ( t (cid:48)(cid:48)(cid:48) , t (cid:48)(cid:48)(cid:48) )] ρ + i (cid:90) t d t (cid:48)(cid:48) Σ I [ G (0) ,< ( t (cid:48)(cid:48) , t (cid:48)(cid:48) )] ρ (cid:90) t (cid:48) d t (cid:48)(cid:48)(cid:48) Σ I [ G (0) ,< ( t (cid:48)(cid:48)(cid:48) , t (cid:48)(cid:48)(cid:48) )] − i ρ (cid:90) t d t (cid:48)(cid:48) Σ I [ G (0) ,< ( t (cid:48)(cid:48) , t (cid:48)(cid:48) )] (cid:90) t (cid:48)(cid:48) d t (cid:48)(cid:48)(cid:48) Σ I [ G (0) ,< ( t (cid:48)(cid:48)(cid:48) , t (cid:48)(cid:48)(cid:48) )] . (B.21)The contribution G (2 , ,< to G (2) ,< is obtained by plugging the first-order correction to the Green’sfunction (B.18) into the self energy Σ I in the first-order expansion of the evolution operator. Theother type of contributions is G (2 , ,> , where Σ I [ G (0) ,> ] enters the second-order term in expansionof the evolution operator.Following this procedure, we first obtain Σ (2) I ,αβ from Eqs. (B.16) and (B.17): Σ (2) I ,αβ ( t ) ≡ Σ I ,αβ [ G (1) ,< ( t , t )]( t ) = (cid:88) µνγ M αβµν M µνγγ ρ γ (cid:16) ρ ν − ρ µ (cid:17) e i( ω α − ω β ) t − e i( ω µ − ω ν ) t ω µ − ω ν . (B.22)For the diagonal element of G <, (2 , αα we get: G <, (2 , αα ( t + τ/ , t − τ/ = e − i ω α τ (cid:90) t + τ/ t − τ/ d t (cid:48)(cid:48) Σ (2) I ,αα ( t (cid:48)(cid:48) ) ρ α (B.23) = i e − i ω α τ (cid:88) µνγ M ααµν M µνγγ ρ α ρ γ ρ µ − ρ ν ω µ − ω ν τ − e i( ω µ − ω ν ) t sin ( ω µ − ω ν ) τ ω µ − ω ν . We see that this diagonal component of the TDHF Green’s function is fully determined by thematrix elements with only three distinct indices.Note that at τ =
0, the contribution to G (2) ,<αα from Eq. (B.23) vanishes, which implies that thecorresponding contribution to the time dependence of the occupation of state α is absent, thusindicating that Eq. (B.23) does not describe decay processes. In order to elucidate the physicalmeaning of this term, we consider its Wigner transform:˜ G (2 , ,<αα ( t , (cid:15) ) = π i (cid:88) µνγ M ααµν M µνγγ ρ α ρ γ ρ µ − ρ ν ω µ − ω ν (B.24) × (cid:40) δ (cid:48) ( (cid:15) − ω α ) + e i( ω µ − ω ν ) t ω µ − ω ν (cid:20) δ (cid:18) (cid:15) − ω α + ω µ − ω ν (cid:19) − δ (cid:18) (cid:15) − ω α − ω µ − ω ν (cid:19)(cid:21)(cid:41) . δ (cid:48) ( (cid:15) − ω α ) shifts the quasi-particle energy and will be omitted in what follows.The other term G (2 , ,< is written in terms of the static self-energy parts Σ (1) : G (2 , ,< I ,αβ ( t + τ/ , t − τ/ = − i (cid:90) t + τ/ dt (cid:48) (cid:90) t (cid:48) dt (cid:48)(cid:48) Σ (1) I ,αγ ( t (cid:48) ) Σ (1) I ,γβ ( t (cid:48)(cid:48) ) ρ β + i (cid:90) t + τ/ dt (cid:48) (cid:90) t − τ/ dt (cid:48)(cid:48) Σ (1) I ,αγ ( t (cid:48) ) ρ γ Σ (1) I ,γβ ( t (cid:48)(cid:48) ) − i (cid:90) t − τ/ dt (cid:48) (cid:90) t (cid:48) dt (cid:48)(cid:48) ρ α Σ (1) I ,αγ ( t (cid:48) ) Σ (1) I ,γβ ( t (cid:48)(cid:48) ) , (B.25)where a summation over the index γ is assumed. For the diagonal elements ( α = β ) we get theWigner transform in the form:˜ G (2 , ,>αα ( t , (cid:15) ) = π i Σ (1) αγ Σ (1) γα ω αγ (cid:26)(cid:20) δ ( (cid:15) − ω α ) − (cid:16) ω αγ t (cid:17) δ (cid:18) (cid:15) − ω α + ω γ (cid:19) − ω αγ δ (cid:48) ( (cid:15) − ω α ) (cid:21) ρ α − (cid:20) δ ( (cid:15) − ω α ) − (cid:16) ω αγ t (cid:17) δ (cid:18) (cid:15) − ω α + ω γ (cid:19) + δ ( (cid:15) − ω γ ) (cid:21) ρ γ (cid:27) , (B.26)where we have introduced a short-hand notation ω αγ = ω α − ω γ . The terms with δ ( (cid:15) − ω α )and δ (cid:48) ( (cid:15) − ω α ) contribute to the quasiparticle weight and the shift of energy, respectively. Thequasiparticle broadening can be expected to arise from the terms that contain other energy levelsover which the summation is performed, like the term with δ ( (cid:15) − ω γ ). We also note that, incontrast to the contribution G (2 , ,<α , the integral over (cid:15) of Eq. (B.26), which yields the Green’sfunction at coinciding time arguments and hence the contribution to the occupation of state α ,does not vanish and hence the contribution of Eq. (B.26) does lead to the TDHF decay. Keepingonly the terms responsible for broadening and with Σ (1) αγ expressed through the matrix elementsaccording to Eq. (B.16), we get Eq. (21) of the main text.It is convenient to define auxiliary functions combining matrix elements with the combina-tions of the density matrices entering the correction to the Green’s function (B.24) and (B.26): f α ( z ) ≡ (cid:88) µνγ M ααµν M µνγγ ρ α ρ γ ρ ν − ρ µ ( ω µ − ω ν ) δ ( ω µ − ω ν − z ) , p α ( z ) ≡ (cid:88) µνγ M αγµµ M γανν ρ α ρ ν ρ µ ω α − ω γ ) δ ( ω α − ω γ − z ) , q α ( z ) ≡ (cid:88) µνγ M αγµµ M γανν ρ γ ρ ν ρ µ ω α − ω γ ) δ ( ω α − ω γ − z ) . (B.27)The second-order TDHF correction G (2 , ,<αα is expressed as an integral in terms of the function f α : ˜ G (2 , ,<α ( t , (cid:15) ) = π i (cid:90) d z f α ( z ) (cid:110) z δ (cid:48) ( (cid:15) − ω α ) − e zt [ δ ( (cid:15) − ω α + z ) − δ ( (cid:15) − ω α − z )] (cid:111) . (B.28)By exchanging µ, ν in the sum in Eq. (B.27), one finds that f α is odd f α ( z ) = − f α ( − z ). Thetime-dependent contribution to ρ α is determined by (cid:90) d z f α ( z ) cos(2 zt ) δ ( (cid:15) − ω α + z ) = − f α ( (cid:15) − ω α ) cos[2( (cid:15) − ω α ) t ] . (B.29)26he integral over (cid:15) , which gives the time decay of the state α is also zero, as was discussed above.The symmetry properties of the term (B.29) are similar to those of the level-shift contributionwith δ (cid:48) ( (cid:15) − ω α ). We conclude that G (2 , ,<α is not responsible to TDHF dephasing at order V .Let us now turn to G (2 , ,< :˜ G (2 , ,<αα ( t , (cid:15) ) = π i (cid:90) d z (cid:8) zt ) (cid:2) p α ( z ) − q α ( z ) (cid:3) δ ( (cid:15) − ω α + z ) + q α ( z ) δ ( (cid:15) − ω α + z ) (cid:9) + . . . , (B.30)where “ . . . ” denote the terms that are not related to the broadening. By definition, functions p α ( z )and q α ( z ) have both even and odd components as functions of z , unlike the odd function f α ( z ).The even components are responsible for the TDHF quasiparticle broadening. Equation (B.30)represents a compact form of Eq. (21) of the main text. Appendix B.2. Second-order self-energy
In this section, we compare the broadening due to TDHF self-energy with the broadeningdue to the second order self-energy that is not included in the TDHF approximation. Whenputting G (0) into the general formula for Σ ( D ) ,< , Eq. (22), one obtains the lowest (second-order)contribution in V : Σ ( D ) ,<αβ ( t , t (cid:48) ) = − (cid:88) i jkl (cid:88) µνγ V il V jk B µν ( l , i ) B ∗ γα ( i , l ) B βγ ( j , k ) B ∗ µν ( k , j )(1 − ρ γ ) ρ µ ρ ν e i( ω γ − ω µ − ω ν )( t − t (cid:48) ) , (B.31)which yields, upon Wigner transformation,˜ Σ ( D ) ,<αβ ( (cid:15), t ) = − π i (cid:88) µνγ M µανγ M βνγµ (1 − ρ γ ) ρ µ ρ ν δ ( (cid:15) + ω γ − ω µ − ω ν ) . (B.32)The exponential ansatz (B.10) is no longer su ffi cient because of the memory integrals. Theimplicit equation to solve isi ∂ t G < ( t , t (cid:48) ) = (cid:16) H + Σ HF ( t ) (cid:17) G < ( t , t (cid:48) ) + Ω ( t , t (cid:48) ) , Ω ( t , t (cid:48) ) = (cid:90) t d t (cid:48)(cid:48) Σ R ( t , t (cid:48)(cid:48) ) G < ( t (cid:48)(cid:48) , t (cid:48) ) − (cid:90) t (cid:48) d t (cid:48)(cid:48) Σ < ( t , t (cid:48)(cid:48) ) G A ( t (cid:48)(cid:48) , t (cid:48) ) . (B.33)We make the ansatz G < ( t , t (cid:48) ) = i U HF ( t , ρ ( t , t (cid:48) ) U HF (0 , t (cid:48) ) , (B.34)where U HF ( t ,
0) is the evolution operator with the TDHF self-energy. With ρ ( t = , t (cid:48) = = ρ , ρ ( t , t (cid:48) ) has to satisfy: i ∂ t ρ ( t , t (cid:48) ) = U HF (0 , t ) Ω ( t , t (cid:48) ) U HF ( t (cid:48) , . (B.35)We are only interested in self-energy terms up to second in V . We therefore plug the zero-orderGF (B.13) into the general expression for the second-order skeleton diagram (22), which yields: Σ ( D ) ,<αβ ( t , t (cid:48) ) = (cid:88) µνγ M αγµν M γβµν (1 − ρ γ ) ρ µ ρ ν e i ( ω γ − ω µ − ω ν )( t − t (cid:48) ) , Σ ( D ) ,<αβ ( t , t (cid:48) ) = (cid:88) γµν M αγµν M γβµν ρ γ (1 − ρ µ )(1 − ρ ν ) e i( ω γ − ω µ − ω ν )( t − t (cid:48) ) . (B.36)27or U HF we can then use the zeroth-order term: U HF , (0) αβ ( t , t (cid:48) ) = e i ω α ( t − t (cid:48) ) δ αβ . (B.37)We thus have to solve the di ff erential equation:i ∂ t ρ (2) ( t , t (cid:48) ) = U HF , (0) (0 , t ) Ω (2) ( t , t (cid:48) ) U HF , (0) ( t (cid:48) , , i ∂ t (cid:48) ρ (2) ( t , t (cid:48) ) = U HF , (0) (0 , t ) Ω (cid:48) , (2) ( t , t (cid:48) ) U HF , (0) ( t (cid:48) , . (B.38)Here, Ω (2) is the term of the second order in V in the expansion of Ω from Eq. (B.33). The t (cid:48) evolution is governed by Ω (cid:48) , (2) : Ω (2) ( t , t (cid:48) ) = (cid:90) t d t (cid:48)(cid:48) Σ R , (2) ( t , t (cid:48)(cid:48) ) G <, (0) ( t (cid:48)(cid:48) , t (cid:48) ) − (cid:90) t (cid:48) d t (cid:48)(cid:48) Σ <, (2) ( t , t (cid:48)(cid:48) ) G A , (0) ( t (cid:48)(cid:48) , t (cid:48) ) Ω (cid:48) , (2) ( t , t (cid:48) ) = − (cid:90) ∞ t (cid:48) d t (cid:48)(cid:48) G <, (0) ( t , t (cid:48)(cid:48) ) Σ R , (2) ( t (cid:48)(cid:48) , t (cid:48) ) + (cid:90) ∞ t d t (cid:48)(cid:48) G A , (0) ( t , t (cid:48)(cid:48) ) Σ <, (2) ( t (cid:48)(cid:48) , t (cid:48) ) , (B.39)where ˜ Ω ( t , t (cid:48) ) = U HF (0 , t ) Ω ( t , t (cid:48) ) U HF ( t (cid:48) ,
0) includes the phase from the Hartree-Fock time evolu-tion. The condition ∂ t (cid:48) ˜ Ω (2) ( t , t (cid:48) ) = ∂ t ˜ Ω (cid:48) , (2) ( t , t (cid:48) ) is satisfied, which means the di ff erential equationis exact and the following integral expression is a solution for ρ (2) ( t , t (cid:48) ): ρ (2) ( t , t (cid:48) ) = (cid:90) t d t ˜ Ω (2) ( t , t (cid:48) ) + (cid:90) t (cid:48) d t (cid:48) ˜ Ω (cid:48) , (2) ( t , t (cid:48) ) − (cid:90) t d t (cid:90) t (cid:48) d t (cid:48) ∂ t (cid:48) ˜ Ω (2) ( t , t (cid:48) ) . (B.40)This leads to the following expression for the diagonal component of the Wigner-transformedGreen’s function ˜ G <, (2 (cid:48) ) αα ( t , (cid:15) ):˜ G <, (2 (cid:48) ) αα ( t , (cid:15) ) = π i (cid:88) µνγ M ναµγ M αµγν ( ω α − ω γ + ω ν + ω µ ) (cid:110) ( ρ γ − ρ ν ρ µ δ ( ω γ − ω ν − ω µ − (cid:15) ) + [ − ρ α ρ γ ( ρ ν + ρ µ − + ρ α ρ ν ρ µ − ( ρ γ − ρ ν ρ µ ] e − i t ( ω α − ω γ + ω ν + ω µ ) δ (cid:18) (cid:15) − ω α + ω γ − ω ν − ω µ (cid:19) − [ − ρ α ρ γ ( ρ ν + ρ µ − + ρ α ρ ν ρ µ ] δ ( (cid:15) − ω α ) − ρ α [ ρ γ ( ρ ν + ρ µ − − ρ ν ρ µ ] δ (cid:48) ( ω α − (cid:15) ) (cid:111) . (B.41)The first two terms here yield the quasiparticle broadening, cf. Eqs. (23), (B.26). Keeping onlythese terms, we get Eq. (23) of the main text. The last two terms, which are proportional to δ ( (cid:15) − ω α ) and δ (cid:48) ( (cid:15) − ω α ), only influence the quasiparticle weight and energies, respectively.To compare this result with the second-order correction to the TDHF Green’s function derivedin Appendix B.1, we retain in Eq. (B.41) only the matrix elements with three distinct indices.By definition the stucture of the matrix elements M νααγ = −M αανγ is the same as in Hartree-Fockcase. We then define auxiliary functions similar to Eq. B.27: r α ( z ) ≡ (cid:88) µν M ααµν M µναα ρ α ρ ν − ρ µ ( ω µ − ω ν ) δ ( ω µ − ω ν − z ) , h α ( z ) ≡ (cid:88) µν M ααµν M µναα ρ α ρ ν ρ µ − ( ρ µ + ρ ν )( ω µ − ω ν ) δ ( ω µ − ω ν − z ) . (B.42)Up to di ff erent prefactors g α behaves like f α from Eq. (B.27) for z >
0. The function h α is aneven function of z and determines the Wigner transform of the second-order self energy in Eq.(B.32). 28part from the structure of density matrices, the resulting expression for ˜ G <, (2) αα ( t , (cid:15) ) is verysimilar to the TDHF self-energy and can also be rewritten with the level-spacing distributionsintroduced in Eq. (B.42):˜ G <, (2 (cid:48) ) αα ( t , (cid:15) ) (cid:39) π i ρ α (cid:90) d z (cid:110) [ h α ( z ) − r α ( z )] δ ( (cid:15) − ω α + z ) − [ h α ( z ) − ρ α r α ( z )] e itz δ ( (cid:15) − ω α + z ) (cid:111) + . . . , (B.43)where “ . . . ” denote the non-broadening terms, as well as the contributions of matrix elementswith four distinct indices. Equation (B.43) is the compact version of Eq. (23) of the main text,where only matrix elements with three distinct indices are retained. The structure of Eq. (B.43)is very similar to that of the TDHF contribution given by Eq. (B.30). We see that, dependingon the initial state, the TDHF contribution to the decay may, in principle, be compensated by thecontribution of the non-Hartree-Fock type. References [1] P. W. Anderson, Absence of di ff usion in certain random lattices, Phys. Rev. 109 (1958) 1492–1505. doi:10.1103/PhysRev.109.1492 .URL https://link.aps.org/doi/10.1103/PhysRev.109.1492 [2] I. V. Gornyi, A. D. Mirlin, D. G. Polyakov, Interacting electrons in disordered wires: Anderson localization andlow- t transport, Phys. Rev. Lett. 95 (2005) 206603. doi:10.1103/PhysRevLett.95.206603 .URL https://link.aps.org/doi/10.1103/PhysRevLett.95.206603 [3] D. Basko, I. Aleiner, B. Altshuler, Metal–insulator transition in a weakly interacting many-electron system withlocalized single-particle states, Annals of Physics 321 (5) (2006) 1126–1205. doi:10.1016/j.aop.2005.11.014 .URL http://dx.doi.org/10.1016/j.aop.2005.11.014 [4] E. Altman, R. Vosk, Universal dynamics and renormalization in many-body-localized systems, Ann. Rev. Cond.Mat. Phys. 6 (1) (2015) 383–409. doi:10.1146/annurev-conmatphys-031214-014701 .URL https://doi.org/10.1146/annurev-conmatphys-031214-014701 [5] R. Nandkishore, D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Ann.Rev. Cond. Mat. Phys. 6 (1) (2015) 15–38. doi:10.1146/annurev-conmatphys-031214-014726 .URL https://doi.org/10.1146/annurev-conmatphys-031214-014726 [6] D. A. Abanin, Z. Papi´c, Recent progress in many-body localization, Ann. Phys. (Berl.) 529 (7) (2017) 1700169. doi:10.1002/andp.201700169 .URL http://dx.doi.org/10.1002/andp.201700169 [7] F. Alet, N. Laflorencie, Many-body localization: An introduction and selected topics, C. R. Phys. 19 (2018) 498. doi:10.1016/j.crhy.2018.03.003 .URL [8] D. A. Abanin, E. Altman, I. Bloch, M. Serbyn, Colloquium: Many-body localization, thermalization, and entan-glement, Reviews of Modern Physics 91 (2) (2019) 021001.[9] D. J. Luitz, Y. B. Lev, The ergodic side of the many-body localization transition, Annalen der Physik 529 (7) (2017)1600350. doi:10.1002/andp.201600350 .URL http://dx.doi.org/10.1002/andp.201600350 [10] S. Gopalakrishnan, S. Parameswaran, Dynamics and transport at the threshold of many-body localization, PhysicsReports (2020).[11] D. J. Luitz, N. Laflorencie, F. Alet, Many-body localization edge in the random-field heisenberg chain, PhysicalReview B 91 (8) (Feb 2015). doi:10.1103/physrevb.91.081103 .URL http://dx.doi.org/10.1103/PhysRevB.91.081103 [12] J. H. Bardarson, F. Pollmann, J. E. Moore, Unbounded growth of entanglement in models of many-body localiza-tion, Phys. Rev. Lett. 109 (2012) 017202. doi:10.1103/PhysRevLett.109.017202 .URL https://link.aps.org/doi/10.1103/PhysRevLett.109.017202 [13] S. P. Lim, D. N. Sheng, Many-body localization and transition by density matrix renormalization group and exactdiagonalization studies, Phys. Rev. B 94 (2016) 045111. doi:10.1103/PhysRevB.94.045111 .URL https://link.aps.org/doi/10.1103/PhysRevB.94.045111
14] B. Kloss, Y. Bar Lev, D. Reichman, Time-dependent variational principle in matrix-product state manifolds: Pitfallsand potential, Phys. Rev. B 97 (2018) 024307. doi:10.1103/PhysRevB.97.024307 .URL https://link.aps.org/doi/10.1103/PhysRevB.97.024307 [15] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, I. V. Gornyi, Many-body localization and delocalization in large quantum chains, Phys. Rev. B 98 (2018) 174202. doi:10.1103/PhysRevB.98.174202 .URL https://link.aps.org/doi/10.1103/PhysRevB.98.174202 [16] E. V. H. Doggen, A. D. Mirlin, Many-body delocalization dynamics in long Aubry-Andr´e quasiperiodic chains,Phys. Rev. B 100 (2019) 104203. doi:10.1103/PhysRevB.100.104203 .URL https://link.aps.org/doi/10.1103/PhysRevB.100.104203 [17] E. V. H. Doggen, I. V. Gornyi, A. D. Mirlin, D. G. Polyakov, Slow many-body delocalization beyond one dimension,Phys. Rev. Lett. 125 (2020) 155701. doi:10.1103/PhysRevLett.125.155701 .URL https://link.aps.org/doi/10.1103/PhysRevLett.125.155701 [18] T. Chanda, P. Sierant, J. Zakrzewski, Time dynamics with matrix product states: Many-body localization transitionof large systems revisited, Phys. Rev. B 101 (2020) 035148. doi:10.1103/PhysRevB.101.035148 .URL https://link.aps.org/doi/10.1103/PhysRevB.101.035148 [19] Y. Bar Lev, D. R. Reichman, Dynamics of many-body localization, Physical Review B 89 (22) (Jun 2014). doi:10.1103/physrevb.89.220201 .URL http://dx.doi.org/10.1103/PhysRevB.89.220201 [20] J. Wurtz, A. Polkovnikov, D. Sels, Cluster truncated wigner approximation in strongly interacting systems, Annalsof Physics 395 (2018) 341–365. doi:10.1016/j.aop.2018.06.001 .URL http://dx.doi.org/10.1016/j.aop.2018.06.001 [21] G. De Tomasi, F. Pollmann, M. Heyl, E ffi ciently solving the dynamics of many-body localized systems at strongdisorder, Phys. Rev. B 99 (2019) 241114(R). doi:10.1103/PhysRevB.99.241114 .URL https://link.aps.org/doi/10.1103/PhysRevB.99.241114 [22] S. A. Weidinger, S. Gopalakrishnan, M. Knap, Self-consistent hartree-fock approach to many-body localization,Phys. Rev. B 98 (2018) 224205. doi:10.1103/PhysRevB.98.224205 .URL https://link.aps.org/doi/10.1103/PhysRevB.98.224205 [23] S. Gopalakrishnan, K. Agarwal, E. A. Demler, D. A. Huse, M. Knap, Gri ffi ths e ff ects and slow dynamics in nearlymany-body localized systems, Phys. Rev. B 93 (2016) 134206. doi:10.1103/PhysRevB.93.134206 .URL https://link.aps.org/doi/10.1103/PhysRevB.93.134206 [24] K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan, D. A. Huse, M. Knap, Rare-region e ff ects and dynamicsnear the many-body localization transition, Annalen der Physik 529 (7) (2017) 1600326. doi:10.1002/andp.201600326 .URL http://dx.doi.org/10.1002/andp.201600326 [25] S. Nandy, F. Evers, S. Bera, Dephasing in strongly disordered interacting quantum wires (2020). arXiv:2010.07919 .[26] T. Kita, Introduction to nonequilibrium statistical mechanics with quantum field theory, Progress of TheoreticalPhysics 123 (4) (2010) 581–658. doi:10.1143/ptp.123.581 .URL http://dx.doi.org/10.1143/PTP.123.581 [27] A. Hindmarsh, L. L. Laboratory, ODEPACK, a Systematized Collection of ODE Solvers, Lawrence LivermoreNational Laboratory, 1982.URL https://books.google.de/books?id=9XWPmwEACAAJ [28] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson,W. Weckesser, et al., Scipy 1.0: fundamental algorithms for scientific computing in python, Nature Methods 17 (3)(2020) 261–272. doi:10.1038/s41592-019-0686-2 .URL http://dx.doi.org/10.1038/s41592-019-0686-2 [29] P. Weinberg, M. Bukov, QuSpin: a Python Package for Dynamics and Exact Diagonalisation of Quantum ManyBody Systems part I: spin chains, SciPost Phys. 2 (2017) 003. doi:10.21468/SciPostPhys.2.1.003 .URL https://scipost.org/10.21468/SciPostPhys.2.1.003 [30] P. Weinberg, M. Bukov, QuSpin: a Python Package for Dynamics and Exact Diagonalisation of QuantumMany Body Systems. Part II: bosons, fermions and higher spins, SciPost Phys. 7 (2019) 20. doi:10.21468/SciPostPhys.7.2.020 .URL https://scipost.org/10.21468/SciPostPhys.7.2.020 [31] E. V. H. Doggen, A. D. Mirlin, Many-body delocalization dynamics in long aubry-andr´e quasiperiodic chains,Physical Review B 100 (10) (Sep 2019). doi:10.1103/physrevb.100.104203 .URL http://dx.doi.org/10.1103/PhysRevB.100.104203 [32] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, F. Verstraete, Unifying time evolution and optimizationwith matrix product states, Physical Review B 94 (16) (Oct 2016). doi:10.1103/physrevb.94.165116 . RL http://dx.doi.org/10.1103/PhysRevB.94.165116 [33] P. P¨opperl, Many-body localization in the self-consistent Hartree-Fock approximation, master Thesis (unpublished)(2020).[34] In fact, if one describes a di ff usive system within the Boltzmann equation, one gets an exponential decay of theimbalance, since it corresponds to a mode with a large wave vector. However, including coupling between the modesoriginating from memory e ff ects associated with impurity scattering, one finds [P. P¨opperl et al, to be published] adi ff usive decay of the imbalance, i.e. t − / in 1D and t − in 2D system. .[35] E. V. H. Doggen, I. V. Gornyi, A. D. Mirlin, D. G. Polyakov, Slow many-body delocalization beyond one dimension,Phys. Rev. Lett. 125 (2020) 155701. doi:10.1103/PhysRevLett.125.155701 .URL https://link.aps.org/doi/10.1103/PhysRevLett.125.155701 [36] S. Gopalakrishnan, D. A. Huse, Instability of many-body localized systems as a phase transition in a nonstandardthermodynamic limit, Phys. Rev. B 99 (2019) 134305. doi:10.1103/PhysRevB.99.134305 .URL https://link.aps.org/doi/10.1103/PhysRevB.99.134305 [37] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. L¨uschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, I. Bloch,Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349 (6250)(2015) 842–845. doi:10.1126/science.aaa7432 .URL http://science.sciencemag.org/content/349/6250/842 [38] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, C. Gross,Exploring the many-body localization transition in two dimensions, Science 352 (6293) (2016) 1547–1552. doi:10.1126/science.aaf8834 .URL http://science.sciencemag.org/content/352/6293/1547http://science.sciencemag.org/content/352/6293/1547