Non-equilibrium correlation dynamics in the one-dimensional Fermi-Hubbard model: A testbed for the two-particle reduced density matrix theory

We explore the non-equilibrium dynamics of a one-dimensional Fermi-Hubbard system as a sensitive testbed for the capabilities of the time-dependent two-particle reduced density matrix (TD2RDM) theory to accurately describe time-dependent correlated systems. We follow the time evolution of the out-of-equilibrium finite-size Fermi-Hubbard model initialized by a quench over extended periods of time. By comparison with exact calculations for small systems and with matrix product state (MPS) calculations for larger systems but limited to short times, we demonstrate that the TD2RDM theory can accurately account for the non-equilibrium dynamics in the regime from weak to moderately strong inter-particle correlations. We find that the quality of the approximate reconstruction of the three-particle cumulant (or correlation) required for the closure of the equations of motion for the reduced density matrix is key to the accuracy of the numerical TD2RDM results. We identify the size of the dynamically induced three-particle correlations and the amplitude of cross correlations between the two- and three-particle cumulants as critical parameters that control the accuracy of the TD2RDM theory when current state-of-the art reconstruction functionals are employed.


I. INTRODUCTION
Accurately describing the correlated out-ofequilibrium dynamics of interacting many-particle systems has remained a great challenge to date.Frequent realizations of such out-of-equilibrium dynamics involve either quenches and relaxation of initially prepared excited states of systems governed by a time-independent Hamiltonian, or systems driven by an explicitly time-dependent Hamiltonian.Such systems are at the forefront of current experimental and theoretical studies (see e.g.[1][2][3][4][5][6][7][8][9][10][11][12][13][14]).Several recent experiments have shown that exotic states of matter can be generated by ultrashort pulses of external fields or energetic ions and that relaxation and decoherence can be strongly influenced by inter-particle correlations [15][16][17][18][19][20][21][22][23].A versatile method to reliably describe the nonequilibrium scenarios of correlated many-body systems, in particular in extended systems and for extended periods of time, is still lacking.Direct many-body wavefunction based methods can be applied only to systems with a moderate number of degrees of freedom and pure states as they eventually face the exponential wall of computational effort when increasing the number of particles and the time interval of propagation [2,9,24,25].Application of the time-dependent density matrix * iva.brezinova@tuwien.ac.at renormalization group (DMRG) theory [2,26] has been shown to yield numerically accurate results, currently, however, limited to one-dimension (1D) systems and short time scales (see e.g.[27][28][29]).Similarly, the closely related time-dependent matrix product state (MPS) method [25,30] invoking the time-dependent variational principle, is also limited to small propagation times for mesoscopic system sizes of a few tens of particles (see e.g.[31]) with the increasing bond dimension as a function of time as the major bottleneck (see e.g.[32]).The complex multi-dimensional information encoded in the quantum many-body wavefunction is, however, often not needed for the extraction of many physical observables.Therefore, an appealing alternative are time-dependent quantum many-body methods that attempt to bypass the use of the many-body wavefunction altogether.Upon successively tracing out more and more degrees of freedom, information and complexity is lost but, in turn, the reduced system is rendered increasingly tractable.A well-known limit of this reduction is the timedependent particle density n(r, t).The corresponding many-body theory, the time-dependent density functional theory (TDDFT) [33,34] with the Kohn-Sham ansatz features a linear scaling with particle number and remains to date the only time-dependent quantum many-body theory applicable to large extended systems with weak to intermediate correlations.Its major drawback, however, is the fundamental lack of knowledge of the exact exchange-correlation (XC) functional.The pathway towards systematic improvements beyond the currently frequently used approximate adiabatic XC functionals is still a widely open question and the applicability of TDDFT to correlated systems is limited.Alternatively, the so-called time-dependent current-density functional theory has been proposed for which, up to now, however only few approximations for the exchange-correlation vector potential have become available [35][36][37].Going up one step of the ladder of reduction the oneparticle reduced density matrix (1RDM) D 1 (r 1 , r 1 , t) allows one to avoid some of the problems of TDDFT [38][39][40] while facing others.
The equation of motion for the 1RDM corresponds to the first equation within the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [41,42] and thus couples the 1RDM to the two-particle reduced density matrix (2RDM) D 12 (r 1 , r 2 , r 1 , r 2 , t).Closing the equations of motion requires representing the 2RDM as a functional of the 1RDM which is challenging in the presence of medium to strong correlations and time-dependent settings.An alternative route to an accurate description of non-equilibrium correlated quantum many-body systems involves non-equilibrium Green's function (NEGF) methods, going back to the pioneering work of Keldysh [43].They have been applied to a wide range of physical systems (see e.g.[44,45] and references therein) but are impeded by a non-linear time scaling, which has only recently been overcome [29,[46][47][48].Moreover, they exhibit a similar hierarchical coupling between different orders of Green's functions, which is subject to closure approximations as in the case of reduced density matrices (see e.g.[44]).The importance of two-particle correlations as imprinted by the pair-wise interaction potentials in most physical systems calls for the use of the 2RDM itself as the fundamental object for representing the many-body system.When only one-and two-body operators are present in the Hamiltonian, the total energy of the system can be exactly expressed in terms of the 2RDM.The fact that the energy is an exactly known functional of the 2RDM has been meanwhile exploited in numerous calculations of groundstate energies in quantum many-body systems [49][50][51][52][53].In this paper we investigate the time-dependent 2RDM.The equation of motion for propagating the 2RDM of an excited system, the second equation of motion the BBGKY hierarchy, requires, the knowledge of the three-particle reduced density matrix (3RDM).Many important works have been devoted in the past to develop reconstruction functionals of the 3RDM in terms of the 2RDM for the quantum many-body ground state problem [54][55][56][57][58][59][60].Incorporating such reconstruction functionals into the time-depending setting within the time-dependent 2RDM method (TD2RDM), we have recently succeeded in calculating the dynamics of multielectron atoms driven by strong laser fields [61,62].Motivated by the stability and remarkable accuracy of this method, it is the aim of the present paper to explore the application of the TD2RDM theory to extended systems, and to systems featuring stronger correlations than typically present in multi-electron atoms.A paradigmatic model system for this endeavor is the Fermi-Hubbard model due to its structural simplicity and the one-parameter tunability from weakly to strongly correlated dynamics.Moreover, this model system can nowadays be realized and accurately probed with ultracold atoms in optical lattices even with single-site resolution (see e.g.[63][64][65][66][67][68][69] and references therein) and is, of course, of conceptual relevance for the study of correlated quantum matter in real solids.Several state-of-the art methods have been tested by application to the Fermi-Hubbard model.They include the NEGF methods [16,28,70], as well as approaches based on Green's functions exploiting the mapping between the Fermi-Hubbard model and an impurity model where the impurity is treated in a fully correlated fashion and is coupled to an external uncorrelated bath.These methods, such as time-dependent dynamical mean-field theory [6] or the explicit sum of a high-order perturbation series in the interaction on the Keldysh contour using quantum Monte-Carlo methods [71][72][73] have the advantage that extended systems can be treated through the coupling of the impurity to an extended bath.However, correlations between distant sites are not well represented.As a prototypical example, we apply the TD2RDM theory to the dynamics of the Fermi-Hubbard model at half filling initialized by a quench, i.e. by suddenly switching off a confining potential that prepares the initial out-of-equilibrium state (Fig. 1).In order to test and to benchmark the TD2RDM we consider in the present work one-dimensional systems with a relatively small number of sites.For these systems a detailed assessment of the accuracy by comparison with numerically exact or highly accurate solutions is still possible allowing us to perform large and systematic parameter scans over many different interaction strengths and excitation energies.We generate (nearly) exact solutions by direct propagation of the Schrödinger equation or using highly accurate matrix product state calculations (MPS) within the time-dependent variational principle [25,30,31].We follow the dynamics over relatively long times (≥ 50 in units of the inverse hopping amplitude) and study the exact build-up of dynamical correlations which can give valuable hints for the applicability of the TD2RDM method as well as for the improvements of reconstruction functionals.We emphasize that our present restriction to 1D systems of moderate size is due to the difficulty of obtaining exact or highly accurate results for comparison, rather than due to the limitations of the TD2RDM theory itself.The latter can be easily extended to larger systems and higher dimensions without encountering major complications.We also compare with time-dependent Hartree-Fock (TDHF) predictions to access the influence of two-particle correlations ne- glected by mean-field theories.We analyze the accuracy of the TD2RDM theory as a function of the strength of the inter-particle interaction as well as the degree of initial excitation.Our focus is on detailed probes of the accuracy of the time-dependent three-particle correlations resulting from different state-of-the art reconstruction functionals.The structure of the paper is as follows: In Sec.II we briefly present the model system under investigation, the one-dimensional Fermi-Hubbard model at half filling.The key ingredients of the TD2RDM theory are reviewed in Sec.III.We numerically analyze the dynamics of two-and three-particle correlations, the so-called cumulants, which are the key ingredient to reconstruction functionals, for small systems by comparison with exact calculations in Sec.IV.Fully self-consistent TD2RDM simulations for the time evolution of the out-of-equilibrium dynamics as monitored by the one-site occupation number are presented in Sec.V, followed by concluding remarks and an outlook to future improvements in Sec.VI.As units we use = m = e = 1 unless otherwise stated.

II. OUT-OF-EQUILIBRIUM FERMI-HUBBARD MODEL
We consider a 1D chain with a number of M s sites (Fig. 1) and impose Dirichlet boundary conditions.The Hamiltonian of the Fermi-Hubbard model in the presence of an external potential initializing the quench is given in second quantization by where i, j denotes nearest-neighbor hopping, J the hopping amplitude, n , the occupation number operators for particles with spin up (down) at site i, U the strength of on-site interaction controlling the correlation energy in the system.V i (t) is the explicitly time dependent potential chosen to be harmonic in the present case, which determines the initial excited state (the ground state of H in the potential for t < 0), and induces the dynamics by a sudden potential quench at t = 0. We consider in the following the spin-symmetric Fermi-Hubbard system at half filling, i.e. particle number N = M s and the number of spin up particles equals to the number of spin down particles (total spin-singlet case).Fig. 1 illustrates the quench-induced dynamics on the level of the one-particle site-occupation number n i corresponding to the diagonal elements of the one-particle reduced density matrix D 1 .The ground state of the interacting many-body system in the potential [Fig. 1 (a)] represents an excited state of the field-free Fermi-Hubbard system and, thus, an out-of-equilibrium state that evolves in time after the quench [Fig. 1 (b)].It would eventually relax,upon coarse graining, to a new equilibrium state.We explore in the following within the framework of TD2RDM theory the importance of interparticle correlations induced by U (Eq. 1) in both the stationary initial state as well the time-dependent correlations induced by the sudden quench.

III. OUTLINE OF TD2RDM THEORY A. Equation of motion
The central object of our method is the 2RDM which is obtained from the exact pure N -body wavefunction |Ψ(t) by tracing out all but two particles.We denote the 2RDM in a basis-independent notation as D 12 and it follows from |Ψ(t) as with N the number of particles, N (N − 1) the normalization related to particles pairs, and Tr 3...N indicating the tracing out of all particles except for the two particles 1 and 2 of interest.More generally, the pRDM is obtained from with normalization factor N !/(N − p)!.
The equation of motion of the 2RDM corresponds to the second equation within the BBGKY hierarchy and reads where the square brackets denote commutators.The Hamiltonian governing Eq. 5 is given (in first quantization) by where h n is the single-particle Hamilton operator, and W nm the two-particle interaction operator.In a basis of spin orbitals {|ψ iσ } Ms i=1 with σ =↑ or σ =↓ localized at a single site i (given e.g. by s-wave orbitals localized at atomic sites in solids or potential minima in optical latices of ultracold atoms) the terms in Eq. 6 yield the explicit matrix representation for the nearest neighbor hopping as and the on-site interaction of particles with different spins as For any initial state (pure or mixed) described by D 12 (t = 0), Eq. 5 allows to propagate the 2RDM without any knowledge of the many-body wavefunction |Ψ(t) .However, since all equations of the BBGKY hierarchy couple to the density matrix of the next higher order, propagation of the 2RDM requires closure, i.e. a sufficiently accurate representation for the 3RDM in terms of the 2RDM.Closure of the equations of motion by reconstruction (denoted by the superscript R in the following) of the 3RDM by the 2RDM, i.e.
poses thus a major challenge for the implementation of the TD2RDM theory as a useful and accurate computational tool.In the spirit of a quantum Boltzmann transport equation [41], we call the term in Eq. 5 containing the D 123 the collision operator (or "collision integral") C, While for the collision operator an approximation to the reconstruction of the 3RDM is required, the timedependent 2RDM, D 12 (t), fully includes all two-particle interactions and correlations without any additional approximation.The solutions of the equations of motion of the 2RDM (Eq.5) feature an important exact relation to Green's functions which opens the door to employ well established diagramatic methods also within the TD2RDM theory.The pRDMs can be identified with the equal-time limits of the p-particle Green's functions G < 1...p (t 1 , . . ., t p , t 1 , . . ., t p ).For the 1RDM and 2RDM, e.g., we get (see e.g.[29,44]) In a given single particle basis D 12 (t) is represented by the matrix Because of the dependence of C on D 123 , the equation of motion of the 2RDM represented in a single-particle basis of dimension M scales as M 7 for a general pairinteraction W .In the present spin-symmetric realization of the Fermi-Hubbard model with equal number of spin-up and spin-down particles, the complexity of the problem can be considerably reduced.The calculation of Eq. 5 can be reduced to that of the spin block D i1↑i2↓ j1↑j2↓ which contains all the information on the entire D 12 (t).
All other spin blocks can be obtained from this particular block either through trivial exchange or spin-flip symmetries, or through the following relation Correspondingly, only the 3RDM block D i1↑i2↑i3↓ j1↑j2↑j3↓ needs to be constructed instead of the entire 3RDM.The equation of motion for D i1↑i2↓ j1↑j2↓ is given in Appendix A. Due to the simple on-site interaction within the Fermi-Hubbard model (Eq.8), the equation of motion for the 2RDM scales as M 4 s .For simplicity of notation, we drop the explicit spin labeling (↑, ↓) unless specifically needed keeping in mind that only the spin blocks identified above need to be calculated.

B. Cumulant expansion
The pRDM describes, in general, the correlated dynamics of a p-tuple of particles embedded in a larger system, in particular in the pure state |Ψ(t) of an N -particle system.In the absence of inter-particle interactions, the pRDM reduces to the independent-particle limit where only Pauli exchange correlations via anti-symmetrization are present.Accordingly, the pRDM can be expanded in term of correlators, in this context conventionally referred to as cumulants [74], of increasing order in the number of particles within the tuple to be correlated with each other.For D 12 the cumulant expansion reads with the two-particle cumulant (or correlator) ∆ 12 and Â the anti-symmetrization operator acting on the two one-particle density matrices D 1 and D 2 .In the singleparticle site representation and Eq. 15 reads The cumulant expansion of D 12 (Eq.15) can be diagrammatically visualized (Fig. 2).The key feature to be noted is that the cumulant expansion [Fig. 2 (a)] does not invoke any ingredients from perturbation theory.The double lines represent the equal-time limit of the full oneparticle propagator.The cumulant represents the sum over all connected diagrams between two one-particle propagators.For illustrative purposes and to connect to other theories we also indicate in Fig. 2 (b) and (c) the corresponding perturbative diagrammatic expansion of the constituents of Fig. 2 (a), the one-particle propagator [Fig. 2 (b)] and the two-particle cumulant [Fig. 2 (c)].We emphasize that within the TD2RDM theory the full 1RDM as well as the full 2RDM are included such that the use of the perturbation series [Fig. 2 (b), (c)] can be avoided.However, these diagrammatic interrelations provide a helpful guidance for developing reconstruction functionals on the three-particle level.
The cumulant expansion of the 3RDM follows as diagrammatically visualized in Fig. 3 (a).The first term in Fig. 3 (a) represents three uncorrelated particles, the second the contribution of two-particle correlations in the presence of a third uncorrelated particle, and the last the true three-particle correlation or three-particle cumulant ∆ 123 containing all connected three-particle diagrams.
For illustrative purposes we show also in Fig. 3 (b) the first few low-order diagrams of a perturbative expansion of ∆ 123 in terms of Hartree-Fock propagators and pair interactions.We note again that the present TD2RDM theory does not make direct use of perturbation theory but we invoke the structure of these diagrams in the following to motivate the approximations of ∆ 123 in terms of one-particle propagators and two-particle cumulants.

C. Three-particle cumulant reconstruction
The challenge to render the TD2RDM theory operational is the closure of the equations of motion (Eq.5) by developing a reconstruction functional for the threeparticle density matrix D R 123 [D 12 ] (Eq.9).The success of the TD2RDM method in describing the many-body dynamics relies on a sufficiently accurate approximation of this functional as has been shown for multi-electron atoms [61,62].While for the non-degenerate ground state the existence of such a reconstruction is assured through Rosina's theorem [75,76], it is presently unknown, whether such an exact reconstruction also exists in a time dependent setting.As Rosina's theorem is an existence theorem, it does not lend itself to aid in the development of new functionals.The cumulant expansion of D 123 (Eq.18) reduces the task of finding a reconstruction functional to that of reconstructing the cumulant ∆ 123 = ∆  [54-56, 59, 60].They provide the starting point of our analysis of the capability of the TD2RDM theory to capture non-equilibrium dynamics in correlated systems.The simplest approximation attributed to Valdemoro (V) and coworkers [54] amounts to neglecting ∆ 123 altogether.Accordingly, the reconstruction functional becomes A similar approximation has been earlier investigated by Wang and Cassing [77].Reconstruction functionals that include contributions from ∆ 123 and benchmarked in this paper have been derived from different perspectives but all rely on approximating ∆ 123 to second order in ∆ 12 .
Nakatsuji and Yasuda (NY) [55] used diagrammatic techniques to arrive at where the intermediate single-particle projector P i is given by with I i the identity matrix, Γ i a diagonal matrix in the eigenrepresentation of the 1RDM with eigenvalues 1 for the lowest N natural orbitals and zero otherwise.Γ i is frequently (in ground state calculations) referred to as the Hartree-Fock reference matrix.We note, however, that in the present context Γ i refers to the natural orbitals of the non-perturbative 1RDM rather than to mean-field states.It has been shown [56] that this projector can be substantially simplified through an expansion in Γ i − D i , the zeroth order of which yields In practice, the summation over the index 2 in Eq. 20 is performed in the basis of natural orbitals with a matrix as the projector containing −1 for unoccupied and 1 for occupied orbitals.We have checked that in the regime where the NY approximation is applicable (see Sec. V below), both Eq.21 and Eq.22 yield very similar results.We, therefore, use the much simpler approximation (Eq.22).Thus, the NY reconstruction functional for D 123 reads (26) It has been shown that in the perturbative limit the reconstructions (Eqs.23,26,25) agree with each other to second-order in the inter-particle interaction [58].
None of the reconstruction functionals presented above preserves, however, important symmetries of the equations of motion (Eq.5), most importantly the contraction consistency (CC).At each instant of time CC requires to hold.We have recently shown [61,62] that the lack of CC seriously impedes the stability as well as the accuracy of the solutions of the equation of motion of the TD2RDM.This deficiency, however, can be cured for any reconstruction functional [61,62] with the important consequence that the orthogonal component of D 123 as well as of ∆ 123 become now exactly known functionals of the 2RDM.This exact functional for three-particle hermitian matrices has been first given in [61] (see also [78] and [29] for a more detailed description).With this decomposition we can now reconstruct parts of the missing components for the above reconstruction functionals through where the defective part of the 2RDM, D d 12 , corresponds to the contraction error in the two-particle space By construction, D R+CC 123 is now contraction consistent, i.e.
Equally importantly, the CC correction to D 123 , , provides a correction to the approximate three-particle cumulant In Eq. 35 we have used the fact that the first two terms of the cumulant expansion (Eq.18) are already exact functionals of D 12 .One remarkable consequence of restoring parts of ∆ 123 by the CC correction is that even the Valdemoro approximation whose bare version (Eq.19) neglects ∆ 123 entirely contains now in its contraction consistent (V+CC) version a three-particle correlation contribution ∆ 123;⊥ .The residual error for the reconstruction functionals considered can thus be traced to the kernel of the three-particle cumulant ∆ 123;K , either completely missing as in the V+CC approximation or only incompletely reconstructed by the NY+CC, TS+CC, or M+CC approximation.In the following, we refer to functionals without the CC correction as the bare functionals.

IV. PROBING THE DYNAMICS OF THE CUMULANTS
The proposed approximate reconstruction functionals for D 123 (t) or, more specifically, for the three-particle cumulant ∆ 123 (t) (Eqs.20,24,26) are at most quadratic functionals in ∆ 12 (t) and local in time.Higher-order terms in ∆ 12 as well as any memory effects are neglected from the outset.As this simple analytic structure of the approximate reconstruction functionals implies strong temporal correlations between ∆ 12 (t) and ∆ 123 (t), it is instructive to probe for the temporal correlations between the time evolution of ∆ 123 (t) and ∆ 12 (t) in the non-equilibrium few-site Fermi-Hubbard model.Only when such time-correlated dynamics is present within the exact solution, the reconstruction by the time-local reconstruction functionals used here can be expected to be accurate.For the Fermi-Hubbard model we explore the coupling between ∆ 12 (t) and ∆ 123 (t) by following the quench dynamics for varying strength of interparticle correlations (Hubbard parameter U ) and strength of the initial outof-equilibrium excitation (controlled by the confining potential parameter V ).We extract from the exact propagation the time evolution of ∆ 12 (t) and ∆ 123 (t) for the quench dynamics without invoking any reconstruction The Frobenius norm provides an upper bound of the largest eigenvalue of M .The square of the Frobenius norm has been used in previous time-dependent studies as a size-extensive measure of correlations [79].
In Fig. 4 we show a typical example for the nonequilibrium dynamics of cumulants at V = 1J and U = 3.1J.All cumulants start with non-zero values ∆ 12 (t = 0) and ∆ 123 (t = 0) of the initial out-ofequilibrium state.They significantly increase immediately following the potential quench signifying the buildup of dynamical correlations in non-equilibrium dynamics.Direct visual inspection reveals that the variations of ∆ ↑↓ 12 , ∆ ↑↑ 12 and ∆ ↑↑↓ 123 are correlated in time with each other.To quantify this time correlations, we calculate the equal-time limit C f g (τ = 0) of the normalized crosscorrelation function with T the total time interval considered, the standard deviation and the mean f = 1 T T t0 dtf (t) (similarly for g).C f g (τ = 0) is also referred to as the Pearson correlation coefficient [80].With this normalization −1 ≤ C f g (τ = 0) ≤ 1 where C f g (τ = 0) = 1(−1) corresponds to perfect (anti-)correlation and C f g (τ = 0) = 0 to absence of correlation in time.The behavior of C f g (τ = 0) for different cumulant pairs in the U -V plane is displayed in Fig. 5 and Fig. 6 for different system sizes (Fig. 5 for M s = 6 sites, Fig. 6 for M s = 8 sites).For the whole parameter scan we use T = 50J −1 , and t 0 = 10J −1 .We use a finite t 0 in Eq. 37 (instead of evaluating the correlation starting with t = 0) because the initial rise of the cumulants is always correlated and its inclusion could lead to an overestimate of the correlation coefficient.We focus here on the long-time average over the fluctuations after the initial build-up.We use t 0 = 10J −1 found to be large enough to separate the initial build-up from the fluctuations around the mean for most parameters in the U -V plane.To delimit and identify structures in the U -V landscape we also display the distribution of two characteristic variables in the U -V plane.One is the ratio of the time-averaged correlation energy Ēcor to the initial degree of excitation parameterized by E pot (0) [Figs.5 (c) and 6 (c)].The latter is given by while which in case of the Fermi-Hubbard model reduces to The other variable measures the build-up of dynamical three-particle correlations during time evolution relative to the three-particle correlations already present in the initial state at t = 0 prior to the quench [Figs.5, 6 (d)], and 6 (a) and (b) marking quite accurately the borderline between strong and weak time correlation (or anticorrelation) between ∆ 123 and ∆ 12 .We also display the borderline between high and low relative correlation energy by plotting the contour line Ēcor = −0.1Epot (0) in Figs. 5 and 6 (a), (b), and (c) which accurately delimits the region of strong time correlation (i.e.C ∆12,∆123 1) in the cumulant dynamics.Obviously, distinct parameter regimes exist for which ∆ 123 (t) and ∆ 12 (t) are strongly correlated with each other: one region pertains to small U (U 0.1J) and a wide range of excitation energies (0 ≤ V 2J).In this region, the cumulants build up over the whole time interval investigated of T = 50J −1 and have not reached saturation for most V .This build-up is naturally strongly correlated over the whole time interval.The other region of positive correlations can be associated with negative relative correlation energies Ēcor −0.1E pot (0) present for the whole interval of U tested and moderate levels of excitation (V J).Furthermore, we find for the larger system (M s = 8) also a region of time-correlation between ∆ 12 (t) and ∆ 123 (t) for positive correlation energy 0.5 1.0 of Ēcor 0.05E pot (0).In other regions the dynamics of ∆ 12 (t) and ∆ 123 (t) is either uncorrelated or even anticorrelated.
In view of the quadratic dependence of the approximate reconstruction functionals of ∆ 123 (t) on ∆ 12 (t) (Eqs.20, 24, 26) the time-correlation maps (Figs. 5 and 6) determined here from exact calculations, allow predictions for the anticipated accuracy of the TD2RDM theory.The time evolution of the many-body system should be captured quite well with the present set of reconstruction functionals in those parameter regions in the U -V plane where the time-correlation between ∆ 12 (t) and ∆ 123 (t) is strong.As will be shown below, the approximate reconstruction functionals are reasonably accurate as long as the build-up of three-particle correlations over time (Eq.42) remains moderate.
To assess the accuracy of the reconstruction functionals locally in time and without the accumulation of errors during time evolution, we also perform exact calculations of both D 12 (t) and D 123 (t) and compare the latter with the reconstructed D R 123 (t) using the exact D 12 (t) as input for the reconstruction, Taking into account the cumulant expansion of D 123 (see Eq. 18) this error coincides with the error in the threeparticle cumulant δ[D 123 ] = δ[∆ 123 ] as only the latter is subject to reconstruction errors.In Fig. 7 we present exemplary results for δ[D 123 (t)] for the ∆ ↑↑↓ 123 block and for the parameters V = 0.8J and U = 1.9J localized in the region of strong temporal correlations [Fig.5 (a), (b)] as well as moderate build-up of three-particle correlations over time [Fig.5 (d)].As the bare Valdemoro approximation neglects ∆ 123 entirely, its error is largest and corresponds to the exact value of ∆ 123 itself.The bare NY, M, and TS perform better with the NY-reconstruction performing best.The difference between NY and TS is very small indicating that the normalization N does not play a significant role in this case.Inclusion of the CC corrections improves the performance of all reconstruction functionals (Fig. 7).As expected, the changes are largest for V+CC for which the CC correction given by the orthogonal component of the cumulant, ∆ 123;⊥ (Eq.35), represents the only contribution to ∆ 123 .For the TS and NY functionals, on the other hand, the corrections due to CC are small in this particular case.
To further probe the accuracy of the reconstruction functionals within the TD2RDM theory locally in time in more detail we now take into account that only a fraction of the elements of the full 3RDM enters the equations of motion of the 2RDM via the collision operator C[D 123 ] (see Eq. 5).We therefore determine the corresponding relative error in the collision operator using the exact input from D exact 12 in D R 123 .The error in the collision operator (shown in Fig. 8 for the ↑↑↓block) mirrors closely that of ∆ 123 (Fig. 7).It is largest for the V functional and smallest for the NY+CC functional.In the following benchmark calculations of the non-equilibrium dynamics of the Fermi-Hubbard model for different pairs of (U , V ) we will restrict ourselves to these two functionals which provide a clear indication of the bandwidth of the expected accuracy.It is furthermore instructive to directly compare the time-local reconstruction error δ[D 123 (t)] (Eq.43) of the V+CC and NY+CC reconstruction functionals for the cumulants with the norm of the cumulants themselves (Fig. 9).Note that the norm ||∆ 123 (t)|| [Fig.7 (a)] coincides with the error of the bare V reconstruction func-  tional in which the three-particle is neglected.In turn, the difference to the V+CC functional [Fig.9 (b)] directly indicates the size of ∆ 123;⊥ included by enforcing contraction consistency.This correction amounts in the present system to an approximate scaling factor of ≈ 0.85 ± 0.05.The time-local reconstruction functional NY+CC improves the reconstruction substantially compared to V+CC.We find that NY+CC performs better in regions where ∆ ↑↑↓ 123 has local minima but performs similarly as the V+CC in regions of local maxima.This gives an indication of current limitations of the reconstruction accuracy and also useful hints for directions of future improvements.We now analyze the time-averaged reconstruction error in the collision operator (Eq.44) in the U -V plane (Fig. 10).The time-averaged error closely mirrors the behavior of the equal-time cross-correlation between the two-particle and three-particle cumulants (Eq.37, Fig. 5).For a Fermi-Hubbard system with weak interparticle interactions U 0.1J the reconstruction error in the collision operator is very small for both the V+CC and the NY+CC reconstruction.For much larger U of up to U ≈ 3J and moderately strong initial excitation (V J) the NY+CC reconstruction performs markedly better.Remarkably, this region is quite faithfully delimited by the region where the build-up over time of the correlations (Eq.42) is moderate, δ∆ ↑↑↓ 123 0.65.Nevertheless, it should be pointed out that the accuracy also of this reconstruction is limited for larger U 3J. Interestingly, in the U -V region of time-anti-correlated or uncorrelated dynamics of the cumulants (Figs. 5 and 6) the time-local reconstruction within NY+CC can cause even larger errors than the V+CC.This is not surprising in view of the fact that the V+CC approximation neglects (apart from the CC correction) ∆ 123 entirely and does not enforce time correlations as the NY+CC reconstruction does through the quadratic dependence of ∆ NY 123 on ∆ 12 (see Eq. 20).Also this observation may point to avenues for further improvements of reconstruction functionals.

V. SELF-CONSISTENT PROPAGATION OF THE 2RDM
We present now examples of the fully self-consistent solution of the equations of motion of the 2RDM (Eq.5) starting from the pure excited state, the many-body ground state in the potential V i (t) prior to the quench at t = 0. We present examples of these simulations for different parameters in the U -V plane marked in Fig. 10.As observable for the quench dynamics we chose the occupation n 1 (t) = D 1 1 (t) of the first site.We compare the results of the TD2RDM theory for a given reconstruction functional with the corresponding exact calculations .System as in Fig. 10.Shown is the time evolution of the occupation of site 1, n1(t) (see Fig. 1), as predicted by TD2RDM using the V+CC and NY+CC functionals and compared with the exact results for different parameter combinations of (U ,V ) as indicated in the frame and marked in Fig. 10 .
(Fig. 11).As a figure of merit we use the time-integrated deviation sensitively probing the amplitude, frequency, and phase of the quench-induced density fluctuations (Fig. 11).We note that the convergent and accurate propagation of the 2RDM requires, in addition to an accurate three-particle cumulant reconstruction functional, also the preservation of N -representability which is a priori not guaranteed when errors due to the approximate reconstruction functionals pile up.N -representability is approximately restored during propagation by purification "on the fly" (see [61,62]).The specific purification algorithm employed in the present simulation is summarized in Appendix C. The time evolution of n 1 (t) for selected values of U and V marked by roman numbers in Fig. 10  Figure 12.System as in Fig. 10.Shown is the time-integrated error δn1 (Eq.45) of the occupation number n1(t) predicted by TD2RDM relative to the exact result in the U -V plane.
The black line denotes the δ∆ ↑↑↓ 123 = 0.65 contour of the buildup of three-particle correlations over time (Eq.42) (same as in Fig. 5).The gray areas mark regions where the method did not meet the convergence criteria discussed in Appendix C.
for t ≤ 50J −1 in Fig. 11.For weak on-site interaction U 0.1J both the V+CC and the NY+CC reconstruction functionals yield excellent agreement with the exact results over a wide range of out-of-equilibrium excitations 0 ≤ V ≤ 1.5J (see also Fig. 12).For stronger U and intermediate quenches with V up to V = 0.8J [Fig.11 (a)-(c)] we observe excellent agreement for the NY+CC reconstruction which performs better than the V+CC reconstruction.For larger V [e.g.V = 1.1J and U = 1.9J,Fig. 11 (d)] deviations for both functionals from the exact result are larger with the V+CC reconstruction functional performing slightly better.To survey the accuracy of the site occupation n 1 (t) we display in Fig. 12 the time-integrated deviations δn 1 from the exact result (Eq.45) for the 6-site Fermi-Hubbard model in the U -V plane.This distribution closely resembles the equal-time correlation between ∆ 12 and ∆ 123 (see Fig. 5, Eq. 37).Obviously, the limitation to moderate build-up of three-particle correlations over time (δ∆ ↑↑↓ 123 0.65) is one reliable predictor for accurate longterm simulations of the correlated non-equilibrium dynamics.We emphasize that the present figure of merit (δn 1 0.1, i.e. the area shaded green or blue in Fig. 10) puts the theory to a fairly stringent test.Even when δn 1 0.3 the agreement with the exact calculation is qualitatively and even semi-quantitatively satisfactory, capturing key features of the fluctuations even though not in all details (see, e.g., Fig. 11 IV).Likewise, when the purification protocol does not fully converge relative to the criteria imposed (see Appendix C), the results for the time-dependent occupation numbers still contain qualitatively correct information on the mean occupation and dominant frequencies.We now turn to larger systems (M s = 18, 20) and time scales for which exact or highly accurate wavefunction based methods (such as MPS) are presently still a challenge.We demonstrate the straightforward applicability of the TD2RDM theory for such systems.For a meaningful comparison with mean-field methods such as TDHF we restrict ourselves to a weakly correlated Fermi-Hubbard model with U = 0.1J but a high degree of excitation.We start with an initial state where all M s /2 sites around the center of the system are doubly occupied amounting to a quench in the limit V → ∞.For the system size M s = 18 [Fig.13 (a)] we can still compare with exact propagation for the time interval 0 ≤ t ≤ 80J −1 .
For M s = 20 [Fig.13 (b)] we can compare to MPS calculations.The latter are, however, limited to short times t 10J −1 .For M s = 18 we find excellent agreement with the exact results up to t ≈ 60J −1 and still reasonable agreement for longer times.For M s = 20 we find excellent agreement with the MPS results for the short time interval for which the MPS data could be gener-ated.We emphasize that increasing the system sizes here from M s = 18 to M s = 20 does not pose any major challenge for TD2RDM theory.The extension towards larger systems approaching extended periodic systems thus appears feasible.We also present in Fig. 13 a comparison with a time-dependent Hartree-Fock (TDHF) simulation as a representative mean-field description within which larger systems are accessible.Even though the system is only weakly correlated, TDHF fails after a short time interval (t 35J −1 ) and substantially overestimates the oscillation amplitude of n 1 (t) (Fig. 13).By contrast the TD2RDM method is able to correctly capture the dynamics in this system for extended periods of time.

VI. CONCLUSIONS AND OUTLOOK
In this paper, we have applied the time-dependent twoparticle reduced density matrix (TD2RDM) theory to the non-equilibrium dynamics of the finite-size Fermi-Hubbard model at half-filling in 1D for a wide range of number of sites and interaction strengths U and initial out-of-equilibrium excitations controlled by the potential strength V of the quench.The Fermi-Hubbard model serves here as a benchmark model to demonstrate the applicability and performance of the theory to extended systems with non-negligible correlations relevant for current research in condensed matter physics and ultra-cold atoms.The TD2RDM theory fully incorporates two-particle correlations and includes approximate three-particle correlations via reconstruction functionals.Key to an accurate description of the dynamics within the TD2RDM theory is the reconstruction of the 3RDM, D 123 , by means of the 2RDM, D 12 , to close the equations of motion, and application of contraction consistency.The underlying assumption of the closure is the existence of a sufficiently accurate reconstruction functional of ∆ 123 .Currently used functionals assume ∆ 123 to be local in time, i.e. they do not take into account possible memory effects.The existence of such a reconstruction functional is guaranteed for the ground state via Rosina's theorem [75,76] but its extension to time-dependent settings is currently unknown.By comparing with exact results for small system sizes we have analyzed the dynamics of both two-and three-particle cumulants.Over a wide range of U and V we could identify parameter regimes in which the dynamics of the three-particle and two-particle cumulants are indeed strongly correlated in time with each other, a key prerequisite for the applicability of current state-of-the art time-local reconstruction functionals.For this particular model system we could show the applicability and accuracy of the TD2RDM theory well into the regime of moderately strong correlations (U 3J), of moderately strong out-of-equilibrium excitations (V J), and for long propagation times (close to hundred time units J −1 ).As an approximate parameter controlling the applicability of TD2RDM theory with the present function-als we could identify the difference δ∆ 123 between the dynamically built up and the initially present (ground state) three-particle correlations.The present observation of the key role of temporal correlations between the two-and three-particle cumulants as well as of the buildup of three-particle correlations over time will provide us with directions for further improvements of the reconstruction functionals.Moreover, they may serve as a guidance for the applicability of TD2RDM theory for systems where exact benchmarks are not available.Application to larger systems indicates that TD2RDM theory is still capable of providing accurate results and may outperform wavefunction based methods.We have showcased an example in the regime of weak interactions and a high degree of excitation in a system with 18 sites, where a numerically exact solution of the full Schrödinger equation is still possible, and with 20 sites where a MPS solution can be generated, however only for a limited time span (≤ 10J −1 ).Based on the excellent agreement with exact results we conclude that the TD2RDM theory has the potential to develop into a versatile tool to study the correlated and strongly driven dynamics of extended models relevant to ultracold atoms and solid state physics.Extensions of the TD2RDM theory to extended systems in two-and three-dimensional systems are planned.steps applied in the iterative process as long the iteration convergences and the smallest geminal occupation number is close to zero (but can still be slightly negative).Since calculating the geminal occupation numbers through exact diagonalization is numerically costly [scaling as O(M 6 s ) with the number of sites] we restrict ourselves to applying Eq.C5 only once each time the smallest geminal occupation number drops below zero for the largest systems in the present paper (with M s = 18 and M s = 20) to save computational time.We have checked that the obtained results are converged with respect to the time step dt of the propagation (i.e. the number of time steps within the whole time interval [0, T ]).This also means that applying different purification schemes (with respect to the threshold on the smallest geminal occupation number and the number of iterative steps) will give the same results on the level of accuracy set by the threshold.Reaching numerical convergence as a function of the number of time steps of the propagation when purification is applied poses a challenge in case of large reconstruction errors.We have used an overall global time step dt, however, within each time step we use a time-adaptive propagation (using a Runge-Kutta-Fehlberg propagator of 4 th and 5 th order) to split each time step into sub-steps within which the prescribed tolerance of the local error is reached.Whenever the smallest geminal occupation number falls below a certain threshold, purification is ap-plied after the global time step.For the scan in Fig. 12 we have applied a threshold of −10 −4 and the maximal allowed number of iterations to reach this threshold is set to 100.We observe that in regions of large errors in the reconstruction the iterative purification often does not reach the threshold (this happens mostly in the lower right corner of Fig. 12).This does not lead necessarily to instabilities.In fact, we did not observe any instabilities for all values of U and V scanned in Fig. 12.However, frequent applications of purification and large numbers of iterations required to reach the threshold or not reaching the threshold at all may induce undesirable numerical noise.This may prevent convergence as a function of the size of the global time step dt (i.e. as a function of the number of time steps used in the propagation) in cases when the uncontrolled noise accumulates (Fig. 12 gray regions).These areas are determined by discarding all results with a local error in the physical observable of > 5 × 10 −3 .Even if the simulation does not meet strict convergence criteria some of the observables are reasonably well represented, e.g. the mean (i.e.time-averaged) occupation number.

Figure 1 .
Figure 1.Fermi-Hubbard model with Ms sites and external harmonic potential (t < 0).(a) The one-particle siteoccupation number ni of the initial ground state in the potential representing an excitation of the potential-free Fermi-Hubbard model after the quench.(b) Snapshots of the time evolved ni(t) at t = 2J −1 and t = 4J−1.The parameters used are V = J and U = J.

Figure 2 .
Figure 2. (a) Diagrammatic representation of the 2RDM.The first term corresponds to ÂD1D2 (the anti-symmetrized contribution is not explicitly shown for brevity).The second term represents the cumulant ∆12 which contains all connected diagrams between two single-particle propagators.Note that no perturbative expansion in the inter-particle interaction W12 is involved in (a).(b) Diagrammatic perturbative expansion of D1, i.e. the equal-time limit of the full single-particle propagator.Two second order diagrams in W12 are shown as illustrative examples.The double line represents the full D1, while each single line in (b) and (c) stands for a Hartree-Fock propagator.(c) Diagramatic perturbative expansion of the cumulant ∆12 with three prototypical diagrams to first and second order in W12.

Figure 3 .
Figure 3. (a) Diagrammatic representation of the 3RDM.The first term corresponds to ÂD1D2D3, the second corresponds to Â∆12D3, and the last to the three-particle cumulant ∆123.The anti-symmetrized contributions of each term are not shown for brevity.(b) Diagrammatic representation of the perturbative expansion of the three-particle cumulant, with several prototypical second/order connected diagrams (∝ W 2 12 ) shown.

Figure 4 .
Figure 4. Time evolution of the exact cumulants (Frobenius norm) for the six-site Fermi-Hubbard model at half filling after the potential quench with V = 1J and U = 3.1J.

Figure 5 .
Figure 5. Six-site Fermi-Hubbard model at half filling (Ms = 6) with interaction U starting from a ground state at varying potential of strength V .Equal-time limit of the normalized cross correlation function between the three-particle cumulant ∆ ↑↑↓ 123 and the two-particle cumulants (a) ∆ ↑↑ 12 and (b) ∆ ↑↓ 12 in the U -V plane.In (c) we depict the mean (i.e.time averaged) correlation energy Ēcor relative to the initial excitation energy Epot(0).The colorbar is cut below and above ±0.1 (i.e.larger and smaller values than shown on the color bar are present.)The black dashed line indicates the −0.1 contour line.In (d) we show the dynamical build-up of three-particle correlations relative to the initial correlations at t = 0, δ∆ ↑↑↓ 123 (Eq.42).The solid black line denotes the δ∆ ↑↑↓ 123 = 0.65 contour.

Figure 6 .
Figure 6.As in Fig. 5, but for the Fermi-Hubbard model with Ms = 8 sites.

Figure 7 .
Figure 7. Six-site Fermi-Hubbard model at half filling (Ms = 6) starting from the many-body ground state in the potential with V = 0.8J defining the out-of-equilibrium initial state in the potential-free Fermi-Hubbard model after the quench.The interaction parameter is U = 1.9J.Error in the three-particle reconstruction δ[D123(t)] (Eq.43) for different reconstruction functionals (a) Valdemoro (V) (Eq.19) with and without contraction consistency (CC), (b) Mazziotti (M) (Eq.26) with and without CC, and (c) Tohyama-Schuck (TS) (Eq.25) and Nakatsuji-Yasuda (NY) (Eq.23) with and without CC.The inset in (c) shows a zoom into a region where the difference between TS+CC and NY+CC becomes visible.

Figure 8 .
Figure 8.As in Fig. 7 however for the relative error of the collision operator δ[C] (Eq.44).Note the different scale of the vertical axis in (c) as compared to (a) and (b).

Figure 10 .
Figure 10.Six-site Fermi-Hubbard model at half filling (Ms = 6) and varying pair interaction U initially (t = 0) confined by the harmonic potential of strength V and released for t > 0. Shown is the reconstruction error in the U -V plane measured by the time-averaged relative error of the collision operator δ[C] (Eq.44) using (a) the V+CC reconstruction functional and (b) the NY+CC reconstruction functional.For the points in the U -V plane marked by (I)-(IV) the time evolution of the occupation number n1(t) is displayed in Fig.11and their error in Fig.12.The black solid line corresponds to the same contour line as in Fig.5.

Figure 11
Figure 11.System as in Fig.10.Shown is the time evolution of the occupation of site 1, n1(t) (see Fig.1), as predicted by TD2RDM using the V+CC and NY+CC functionals and compared with the exact results for different parameter combinations of (U ,V ) as indicated in the frame and marked in Fig.10.

Figure 13 .
Figure13.The Fermi-Hubbard model with (a) Ms = 18 sites and (b) Ms = 20 sites at half filling and U = 0.1J.All particles are initially located at the center of the system amounting to a potential quench of V → ∞.The exact solution of the Schrödinger equation in (a) has been obtained using a Trotter decomposition of the sparse time-evolution operator (time step dt = 0.025J −1 ), and via the one-site timedependent variational principle with matrix product state bond dimension of χ = 512 (dt = 0.025J −1 ) in (b).Convergence with bond-dimension and time step has been verified to be within an absolute precision of 10 −3 for the expectation values shown.We compare with the TD2RDM prediction using the V+CC and NY+CC reconstruction functionals as well as to the mean-field solution within time-dependent Hartree-Fock (TDHF) theory.
R 123 [D 12 ] as the other terms contributing to D 123 are already known functionals of D 12 (and D 1 ).Several approximate functionals D R 123 [D 12 ] or ∆ R 123 [D 12 ] have been recently proposed