Enforcing conservation laws in nonequilibrium cluster perturbation theory
EEnforcing conservation laws in nonequilibrium cluster perturbation theory
Christian Gramsch
1, 2 and Michael Potthoff
1, 2 I. Institute for Theoretical Physics, University of Hamburg, Jungiusstraße 9, 20355 Hamburg, Germany The Hamburg Centre for Ultrafast Imaging, Luruper Chaussee 149, 22761 Hamburg, Germany
Using the recently introduced time-local formulation of the nonequilibrium cluster perturbation theory (CPT),we construct a generalization of the approach such that macroscopic conservation laws are respected. This isachieved by exploiting the freedom for the choice of the starting point of the all-order perturbation theoryin the inter-cluster hopping. The proposed conserving CPT is a self-consistent propagation scheme whichrespects the conservation of energy, particle number and spin, which treats short-range correlations exactly upto the linear scale of the cluster, and which represents a mean-field-like approach on length scales beyond thecluster size. Using Green’s functions, conservation laws are formulated as local constraints on the local spin-dependent particle and the doublon density. We consider them as conditional equations to self-consistentlyfix the time-dependent intra-cluster one-particle parameters. Thanks to the intrinsic causality of the CPT, thiscan be set up as a step-by-step time propagation scheme with a computational effort scaling linearly with themaximum propagation time and exponentially in the cluster size. As a proof of concept, we consider thedynamics of the two-dimensional, particle-hole-symmetric Hubbard model following a weak interaction quenchby simply employing two-site clusters only. Conservation laws are satisfied by construction. We demonstratethat enforcing them has strong impact on the dynamics. While the doublon density is strongly oscillating withinplain CPT, a monotonic relaxation is observed within the conserving CPT.
PACS numbers:
I. INTRODUCTION
A major challenge of the theory of strongly correlatedlattice-fermion models is to predict the real-time dynam-ics of local observables on a long time scale.
Usingexact-diagonalization techniques, i.e., full diagonalization orKrylov-space methods, only lattices with a small number ofsites can be addressed such that artificial boundary effects startto dominate the dynamics after a few elementary hopping pro-cesses. Much larger systems are in principle accessible bymeans of quantum Monte-Carlo methods, at least in thermalequilibrium. As concerns real-time dynamics, however, thesign (or complex phase) problem still prevents a computa-tionally efficient simulation, even for impurity-type modelswhich are typically sign-problem-free at thermal equilibrium,and despite substantial progress in the recent past.
Forimpurity and for one-dimensional systems, recent extensionsof the numerical renormalization group and of the density-matrix renormalization group to the time domain havebeen shown to be highly efficient and accurate.For lattice models in two or higher dimensions, on the otherhand, one has to resort to approximations, e.g., to the time-dependent variational principle evaluated with Gutzwiller or with Jastrow-like variational wave functions. Using aGreen’s-function-based approach, on the other hand, one mayalso treat the problem within weak-coupling perturbation the-ory. Naive perturbative techniques, however, usually violatethe macroscopic conservation laws that result from the con-tinuous symmetries of the lattice-fermion model. As has beenshown by Baym and Kadanoff, “conserving approxima-tions” can be constructed diagrammatically by deriving theself-energy from a (truncated) Luttinger-Ward functional involving, e.g., certain infinite re-summations of diagramclasses, and by calculating the single-particle Green’s functionself-consistently. Due to the necessary approximation of the Φ functional, however, certain low-order diagrams are neglectedwhich implies that, strictly speaking, such conserving approx-imations are usually restricted to the weak-coupling limit. Nonperturbative conserving approximations can eitherbe constructed with the help of the many-body wavefunction, or, using Green’s functions, within the frame-work of the nonequilibrium generalization of self-energy-functional theory (SFT). Here, the Green’s function isself-consistently obtained from an optimal self-energy whichmakes the grand potential of the initial thermal state, ex-pressed as a functional of the nonequilibrium self-energy, sta-tionary. The equilibrium SFT comprises different approxima-tions, such as the variational cluster approximation andthe dynamical impurity approximation. These techniqueshave been extended to real-time dynamics and have been ap-plied recently to study the dynamical Mott transition in theHubbard model and a variant of the periodic Andersonmodel. Another nonperturbative conserving approach, which canbe derived within the SFT framework but has actually beenproposed much earlier, is the (nonequilibrium) dynamicalmean-field theory.
Being the exact theory in the limitof infinite spatial dimensions, conservation laws are in princi-ple naturally satisfied in this case. In practice, however, thisrequires the exact solution of a highly nontrivial quantum-impurity model out of equilibrium. First cluster extensionsof the DMFT have been reported as well.
Those combinethe mean-field concept with an improved description of spatialcorrelations.The nonequilibrium extension of cluster-perturbationtheory (CPT) is a strongly simplified variant of a cluster-embedding approach. Still, the numerical solution of the basicCPT equation is complicated by the presence of memory ef-fects which are encoded in real-time Green’s functions withinthe Keldysh formalism. This is very similar to the nonequilib- a r X i v : . [ c ond - m a t . s t r- e l ] M a y rium Dyson or Kadanoff-Baym equations in other diagram-matic approaches. As has been shown recently, how-ever, the problem can be mapped exactly onto a noninteractingproblem with additional auxiliary degrees of freedom. Adopt-ing this idea, we could demonstrate that the CPT real-timedynamics can be understood as a simple Markovian dynamicsof a system of noninteracting fermions but in a much largertime-dependent bath of virtual degrees of freedom. Using thisreformulation of the CPT, it has been possible to formallystudy the real-time dynamics of an inhomogeneous setup inthe two-dimensional Hubbard model consisting of 10 × where the inverse nearest-neighbor hopping serves as the time unit.Those plain CPT calculations, however, suffer from a cou-ple of conceptual problems. The drawback of any mean-fieldtheory is the missing feedback of certain correlations on thedynamics of the observables of interest, such as, e.g., themissing feedback of nonlocal spatial correlations on the lo-cal self-energy in the case of the DMFT. In the case of plainCPT, the situation is even worse as there is no feedback atall. In particular, plain CPT calculations cannot be expectedto respect the macroscopic conservation laws emerging fromthe symmetries of the underlying Hamiltonian. This can betraced back to the fact that the plain CPT does not containany element of self-consistency. Therefore, it is not surpris-ing that a violation of, e.g., total-energy conservation has beenobserved. With the present study we give a proof of principle that thisdrawback can be overcome. We make use of the fact that theCPT can be viewed as an all-order perturbation theory inthe inter-cluster hopping around a system of decoupled clus-ters, where the starting point, i.e., the intra-cluster Hamil-tonian, is not at all predetermined. The idea is to formu-late the macroscopic conservation laws as local constraintson the spin-dependent particle and doublon density. Theseequations are then used to fix the intra-cluster one-particleparameters and thereby to optimize the starting point for thecluster-perturbation expansion. This defines a novel “conserv-ing cluster-perturbation theory.” The theory is conserving byconstruction, it is nonperturbative, and in principle controlledby the inverse cluster size as a small parameter. In practice,however, the accessible cluster size is limited by the expo-nential growth of the cluster Hilbert space. Hence, conservingCPT must be seen as a typical cluster mean-field theory whichcorrectly accounts for nonlocal correlations up to the linearscale of the cluster. Opposed to standard mean-field theories,the “mean-field” or the renormalization of the one-particle pa-rameters is determined by imposing local constraints express-ing conservation laws, i.e., it is finally the symmetries of thelattice model which dictates the time-dependent cluster em-bedding. As the theory relies on local self-consistency or con-ditional equations, it can easily be extended to inhomogeneousmodels or inhomogeneous initial states.While the underlying idea is conceptually simple, its prac-tical realization requires a couple of new theoretical conceptswhich are discussed here in detail. In particular, the imple-mentation of a causal time-stepping algorithm requires a care-ful analysis to which order the renormalization of the intra- cluster parameters at a certain time slice enters the conditionalequations. We are able to demonstrate that an efficient numer-ical implementation of the theory is possible and discuss firstresults for weak interaction quenches in a two-dimensionalHubbard model. The algorithm scales linearly with the prop-agation time and exponentially in the cluster size. Conser-vation laws are satisfied with numerical accuracy. Yet, longtime scales cannot be achieved with the present implementa-tion due to singular points which are found to evolve duringthe time propagation.The next section briefly states the model and the necessaryelements of the Keldysh formalism. Section III introduces theCPT and discusses the formulation of the local constraints.The mapping onto a noninteracting auxiliary problem is de-scribed in Sec. IV. The main theoretical work addresses thesolution of the local constraints for the optimal starting pointof the all-order perturbation theory. This is presented in Sec.V. Numerical results are discussed in Sec. VI, and the conclu-sions are summarized in Sec. VII. II. MODEL AND NONEQUILIBRIUM FORMALISM
We consider the single-band, fermionic Hubbard model onan arbitrary lattice with a time-dependent hopping matrix T ( t ) and interaction strength U ( t ) . The hopping is assumed asspin-diagonal for simplicity. The Hamiltonian reads H T , U ( t ) = ∑ i j σ ( T i j σ ( t ) − δ i j µ ) c † i σ c j σ + U ( t ) ∑ i n i ↑ n i ↓ , (1)where the operators c † i σ ( c j σ ) create (annihilate) a fermionwith spin σ ∈ {↑ , ↓} at site i ( j ), and where n i σ = c † i σ c i σ de-notes the spin-dependent local density operator. At time t = β and chemical potential µ . Nonequilib-rium real-time dynamics for t > or experimentswith ultracold gases in optical lattices. Our central quantity of interest is the one-particle Green’sfunction which is defined as [ G T , U ] i j σ ( t , t (cid:48) ) = − i (cid:104)T C ˆ c i σ ( t ) ˆ c † j σ ( t (cid:48) ) (cid:105) H T , U (2) ≡ − iZ tr (cid:16) exp ( − β H T , U ( )) (cid:104) T C ˆ c i σ ( t ) ˆ c † j σ ( t (cid:48) ) (cid:105)(cid:17) . Here “tr ( . . . ) ” traces over the Fock space, i.e., wetake averages using the grand-canonical ensemble. Z = tr ( exp ( − β H T , U ( ))) defines the grand-canonicalpartition function and T C the time-ordering operator onthe Keldysh-Matsubara contour C . The time variables t and t (cid:48) are thus understood as contour times. An in-depthintroduction to the Keldysh formalism can be found inRefs. 47, 48. Throughout the text we use the convention thatoperators with a hat carry a time dependence according tothe Heisenberg picture, i.e., ˆ c i ( t ) = U † ( t , ) c i U ( t , ) , where U ( t , t (cid:48) ) = T exp (cid:0) − i (cid:82) tt (cid:48) H T , U ( t ) d t (cid:1) , for t > t (cid:48) , is the system’stime-evolution operator and T the time-ordering operator.The dependence of the Green’s function on T ( t ) and U ( t ) is made explicit in the notation using subscripts, i.e., G T , U ,where convenient. A similar notation is also used for otherquantities.Through Dyson’s equation, the Green’s function is linkedto the self-energy G T , U = G T , + G T , ◦ Σ T , U ◦ G T , U , (3)where G T , is the noninteracting propagator. Its (contour) in-verse is [ G − T , ] i j σ ( t , t (cid:48) ) = [ δ i j ( i ∂ t + µ ) − T i j σ ( t ))] δ C ( t , t (cid:48) ) . (4)In Eq. (3) we made use of the shorthand notation “ ◦ ” for theconvolution of contour matrices. In particular [ Σ T , U ◦ G T , U ] i j σ ( t , t (cid:48) ) = (cid:90) C d t ∑ l [ Σ T , U ] il σ ( t , t )[ G T , U ] l j σ ( t , t (cid:48) ) . (5)Note that in this context we implicitly assume a contour Diracdelta function δ C ( t , t (cid:48) ) present in case of time-local quanti-ties. For example, T ( t ) should be replaced by T ( t ) δ C ( t , t (cid:48) ) ina contour convolution, so that [ T ◦ G T , U ] i j σ ( t , t (cid:48) ) = (cid:90) C d t ∑ l T il σ ( t ) δ C ( t , t )[ G T , U ] l j σ ( t , t (cid:48) )= ∑ l T il σ ( t )[ G T , U ] l j σ ( t , t (cid:48) ) . (6)Combining Dyson’s equation with the equation of motionfor the one-particle Green’s function, we get [ Σ T , U ◦ G T , U ] i j σ ( t , t (cid:48) ) = − iU ( t )[ G ( l ) T , U ] i j σ ( t , t (cid:48) ) , (7)where G ( l ) is the two-particle Green’s function [ G ( l ) T , U ] i j σ ( t , t (cid:48) ) = (cid:104)T C ˆ n i ¯ σ ( t ) ˆ c i σ ( t ) ˆ c † j σ ( t (cid:48) ) (cid:105) , (8)and where ¯ σ indicates a flip of the spin index σ , i.e., ¯ ↑ = ↓ andvice versa. Analogously to Eq. (8), we furthermore have − i [ G ( r ) T , U ] i j σ ( t , t (cid:48) ) U ( t (cid:48) ) = [ G T , U ◦ Σ T , U ] i j σ ( t , t (cid:48) ) , (9)where G ( r ) is defined as [ G ( r ) T , U ] i j σ ( t , t (cid:48) ) = (cid:104)T C ˆ c i σ ( t ) ˆ n j ¯ σ ( t (cid:48) ) ˆ c † j σ ( t (cid:48) ) (cid:105) . (10)The local doublon density d i ( t ) can be expressed through thetwo-particle Green’s functions as d i ( t ) ≡ (cid:104) ˆ n i ↑ ( t ) ˆ n i ↓ ( t ) (cid:105) H T , U (11) = − [ G ( l ) T , U ] ii σ ( t , t + ) = − [ G ( r ) T , U ] ii σ (cid:48) ( t , t + ) , with t + being infinitesimally “later” than t in the sense of timeordering on the Keldysh-Matsubara contour. III. PREPARATIONS
To enforce conservation laws within cluster perturbationtheory (CPT), we proceed in two steps. First, we work outthat the CPT approach is not unique and that there are freeparameters at one’s disposal. Second, to fix these parameters,we suggest to employ local constraints expressing the con-servation laws that result from continuous symmetries of theHubbard model. We start with a discussion of the main ideaof the CPT and of the local constraints on the spin-dependentparticle and doublon density.
A. Conventional cluster perturbation theory
The idea of the CPT is to partition the lattice into clus-ters small enough to be treated exactly, e.g., using Krylov-space methods or full diagonalization, and to subsequently in-clude the connections between the clusters perturbatively. Onthe level of the Hamiltonian one starts by partitioning the fullhopping matrix T into the intra-cluster hopping T (cid:48) and theinter-cluster hopping V so that T (cid:48) only contains terms whichconnect lattice sites within the individual clusters, while V contains all remaining terms such that T = T (cid:48) + V , see Fig. 1.Corresponding to the intra-cluster hopping, we define a clusterHamiltonian H T (cid:48) , U ( t ) which describes the system of isolatedclusters, also referred to as the reference system. Its Green’sfunction and self-energy are denoted as G T (cid:48) , U and Σ T (cid:48) , U , re-spectively.For the equilibrium as well as for the nonequilibriumcase, the CPT can be seen as an all-order perturbationtheory in the inter-cluster hopping V which provides the one-particle Green’s function of the original system by expandingaround the cluster Green’s function: G CPT = G T (cid:48) , U + G T (cid:48) , U ◦ V ◦ G T (cid:48) , U + · · · = G − T (cid:48) , U − V . (12)We also have: G CPT = G − T , − Σ T (cid:48) , U . (13)In the noninteracting case, this is exact since Σ T (cid:48) , =
0. For fi-nite U ( t ) , however, the CPT Green’s function G CPT representsan approximation of the exact Green’s function G T , U .A closer look reveals that the CPT is not unique since onemay consider a different starting point for the all-order pertur-bation theory in V . To this end, consider a starting point witha renormalized intra-cluster hopping, T (cid:48) → T (cid:48) − λ , resultingin a renormalized cluster Green’s function G T (cid:48) − λ , U and self-energy Σ T (cid:48) − λ , U . Correspondingly, also the inter-cluster hop-ping V must be renormalized as V → V + λ . Summation ofthe geometrical series yields G CPT [ λ ] = G − T (cid:48) − λ , U − ( V + λ ) = G − T , − Σ T (cid:48) − λ , U , (14)where we emphasized the special role of the renormaliza-tion parameter λ by square brackets. For U ( t ) =
0, we have
FIG. 1: Sketch of plain CPT (a) and conserving CPT (b).
PlainCPT: the hopping matrix T is decomposed as T = T (cid:48) + V into theintra-cluster ( T (cid:48) ) and the inter-cluster hopping V . The problem de-fined by T (cid:48) (and the local Hubbard-type interaction) is solved ex-actly. V is treated by all-order perturbation theory (neglecting ver-tex corrections), see Eq. (12). Conserving CPT: the same as plainCPT but with “renormalized” intra- ( T (cid:48) − λ ) and inter-cluster hop-ping V + λ , where the time-dependent renormalization λ (indicatedin red) is used to enforce conservation laws. Note that λ may alsocomprise the on-site energies. G CPT [ λ ] = G T , for any λ . For an interacting system, how-ever, the choice for λ is crucial, i.e., the resulting CPT Green’sfunction does depend on the starting point of the all-order per-turbation theory in the inter-cluster hopping.This ambiguity in the definition of the CPT seems to beproblematic on first sight, yet it can be turned into an advan-tage by interpreting the renormalization λ as an optimizationparameter. This has been worked out systematically in thecontext of the (nonequilibrium) self-energy functional theory(SFT), where the optimal λ is derived from a varia-tional principle based on the self-energy. Here, we will take adifferent route and use the freedom in λ to enforce the localconstraints on spin-dependent particle and doublon density.Physically, the parameter set λ must be interpreted as a non-local mean-field and the resulting conserving CPT as a clustermean-field theory. B. Formulation of the conservation laws as local constraints
While conservation laws like particle-number or energyconservation are naturally fulfilled if one is able to treat aphysical problem exactly, this is not necessarily the case whenworking with approximate methods. For Green’s-function-based methods it was shown by Baym and Kadanoff that respecting certain symmetry relations for the two-particleGreen’s function is sufficient to ensure that an approximationis conserving.Here, we build on an equivalent formulation of the macro-scopic conservation laws for the particle number, spin and en-ergy and reformulate them as local constraints for the spin-dependent particle density and the doublon density, respec-tively. This is in the spirit of expressing conservation laws ofa classical field theory as continuity equations and follows thework of Baym and Kadanoff.
One should note, however, that in our case the local constraints cannot be written in thestandard form of continuity equations, as here we aim at anapproach for a discrete lattice model.To discuss the local constraints, we first consider the ex-act time evolution of a system described by the HubbardHamiltonian H T , U ( t ) . We write G ≡ G T , U , G ( l ) ≡ G ( l ) T , U and G ( r ) ≡ G ( r ) T , U in this subsection to keep the notation simple.The exact time evolution of the system will preserve the totalparticle number and the z -component of the total spin as canbe expressed by the following local constraint for the spin-dependent density:0 = ∂ t (cid:104) ˆ n i σ ( t ) (cid:105) H T , U − [ G ◦ T − T ◦ G ] ii σ ( t , t + ) , ⇔ F i σ ( t ) ≡ G ( l ) ii σ ( t , t + ) − G ( r ) ii σ ( t , t + ) = , (15)as can be verified directly using Eq. (11).The first line of Eq. (15) constitutes the discrete-lattice ana-log of the continuity equation for the spin-dependent parti-cle density. Opposed to a continuum theory, however, the di-vergence of the spin-dependent particle-current density is re-placed by the commutator. The second line of Eq. (15) is anequivalent formulation of the same constraint as has originallybeen mentioned by Baym and Kadanoff. Next, we consider the following local constraint for thedoublon density [cf. Eq. (11)]: C i σ ( t ) ≡ i ∂ t (cid:104) G ( l ) ii σ ( t , t + ) + G ( r ) ii σ ( t , t + ) (cid:105) (16) − ∑ j σ (cid:104) T i j σ ( t ) G ( r ) ji σ ( t , t + ) − G ( l ) i j σ ( t , t + ) T ji σ ( t ) (cid:105) = . In the exact theory, this constraint together with the aboveconstraint F i σ ( t ) = G ( l ) or G ( r ) and for each spin component σ in Eq. (11).More important, in case of a time-independent Hamilto-nian, i.e., if H T , U ( t ) = const . for t > t , Eq. (16) implies total-energy conservation. This is explicitly shown in the AppendixA where, for completeness, also a formal derivation of Eq.(16) is carried out.While in the exact theory the equations F i σ ( t ) = C i σ ( t ) = G ( l ) and G ( r ) , can be expressed in terms of the single-particleGreen’s function and the self-energy and thus both equations F i σ ( t ) = C i σ ( t ) = F i σ ( t ) = C i σ ( t ) = λ = λ opt . If λ opt can be found, this automatically en-sures the conservation of particle number, spin and energy. IV. HAMILTONIAN-BASED FORMULATION
In the last section we have introduced the CPT in its usualform, i.e., based on the self-energy Σ T (cid:48) − λ , U of the referencesystem and Dyson’s equation. A major drawback of this ap-proach is its limitation for the maximum propagation time thatcan be reached in a practical numerical calculation. This isdue to the fact that the CPT Green’s function and the self-energy of the reference system are nonlocal in time throughtheir dependence on two contour times. The necessary stor-age for these quantities scales quadratically with the maxi-mum propagation time, the effort for solving Dyson’s equa-tion scales cubically. This intrinsic limitation can be over-come if a so-called Lehmann representation of the self-energyis available. This allows us to solve the Dyson equation bymeans of a Markovian propagation scheme which permits toreach much longer time scales. In the following, we con-sider this Lehmann representation of the self-energy as given.Its existence for an arbitrary, fermionic lattice system out ofequilibrium has been shown in Ref. 39 where also a construc-tive numerical scheme has been presented. It can be used incase of small clusters accessible to exact-diagonalization tech-niques. In the following we briefly recall the main results andthen discuss the application to conservation laws and the re-spective local constraints.
A. Convolution-free definition of G ( l ) and G ( r ) The nonequilibrium self-energy Σ ≡ Σ T , U of any lattice-fermion model has a unique Lehmann representation: Σ i j σ ( t , t (cid:48) ) = δ C ( t , t (cid:48) ) Σ HF i j σ ( t ) + ∑ s σ h is σ ( t ) g ( h ss σ ; t , t (cid:48) ) h ∗ js σ ( t (cid:48) ) . (17)Here, Σ HF i j σ ( t ) is the time-local Hartree-Fock term. The secondterm has a hybridization-function-like structure where h is σ ( t ) denotes the hopping between a physical site i andan additional virtual site labeled by the index s . The time-independent on-site energy of the virtual site is given by h ss .Furthermore, g ( ε ; t , t (cid:48) ) is the noninteracting Green’s functionof an isolated one-particle mode ( h mode = ε c † c ) with excita-tion energy ε : g ( ε ; t , t (cid:48) ) = i [ f ( ε ) − Θ C ( t , t (cid:48) )] e − i ε ( t − t (cid:48) ) . (18)Here, f ( ε ) = ( e βε + ) − denotes the Fermi-function, and Θ C ( t , t (cid:48) ) refers to the contour variant of the Heaviside stepfunction, i.e., Θ C ( t , t (cid:48) ) = t ≥ C t (cid:48) , and Θ C ( t , t (cid:48) ) = H eff ( t ) = ∑ i j σ (cid:0) T i j σ ( t ) − δ i j µ + Σ HF i j σ ( t ) (cid:1) c † i σ c j σ (19) + ∑ is σ (cid:16) h is σ ( t ) c † i σ c s σ + h . c . (cid:17) + ∑ s h ss σ c † s σ c s σ , FIG. 2: The effective, one-particle Hamiltonian, (21), has three dis-tinct kind of elements. The hybridization elements h is σ ( t ) = h ∗ si σ ( t ) ,the elements h i j σ ( t ) of the physical sector and the time-independentelements h ss (cid:48) σ ∝ δ ss (cid:48) of the virtual sector. which reproduces the one-particle Green’s function G ≡ G T , U exactly when evaluated at the physical sites i , j . We empha-size that H eff ( t ) includes all correlation effects through thehybridization strengths h is σ ( t ) and on-site energies h ss σ of thevirtual sites. Furthermore, the previously mentioned interpre-tation of the s -degrees of freedom as additional virtual sitesbecomes obvious in Eq. (19).For a given arbitrary self-energy one would typically haveto consider a continuum of virtual sites. Here, however, thesituation is much simpler since Σ i j σ ( t , t (cid:48) ) is the CPT self-energy, i.e., the self-energy of our reference model consistingof a system of decoupled clusters. In this case the total numberof physical and of necessary virtual sites equals the number ofsingle-particle excitations with nonzero weight. As the lat-ter grows exponentially with the size of the individual cluster,the exact mapping and the exact numerical construction of theeffective Hamiltonian is limited to clusters small enough to al-low for an exact numerical diagonalization. The computationof the parameters of the effective Hamiltonian (19) is non-trivial but straightforward and numerically completely stable.Details are described in Ref. 39.For clarity, we use the following convention throughout thepaper: physical sites: i , j , k , l , (20)virtual sites: s , s (cid:48) , physical or virtual sites: x , y . Defining h i j σ ( t ) ≡ T i j σ ( t ) − δ i j µ + Σ HF i j σ ( t ) , the effectiveHamiltonian can be written as (cf. Fig. 2) H eff ( t ) = ∑ xy σ h xy σ ( t ) c † x σ c y σ . (21)The corresponding Green’s function is given on the physicalbut also on the virtual sites, G xy σ ( t , t (cid:48) ) = − i (cid:104)T C ˆ c x σ ( t ) ˆ c † y σ ( t (cid:48) ) (cid:105) H eff , (22)so that the original, physical Green’s function G i j σ ( t , t (cid:48) ) is obtained if we restrict x , y to physical sites only, i.e., ( x , y ) = ( i , j ) .Many-particle correlation functions, e.g., spin-spin correla-tions, are in general not accessible from the effective Hamil-tonian. There is, however, one important exception. Namely,the two-particle Green’s functions G ( l ) and G ( r ) can be ex-pressed as contour convolutions of the system’s self-energywith the Green’s function, cf. Eqs. (7) and (9). In theHamiltonian-based formalism this convolution is greatly sim-plified and becomes a straightforward matrix multiplication.By comparing Dyson’s equation (3) with the equation of mo-tion that follows from Eq. (22), one readily finds the identity [ Σ ◦ G ] i j σ ( t , t (cid:48) ) = ∑ l [ h il σ ( t ) − T il σ ( t ) + δ il µ ] G l j σ ( t , t (cid:48) )+ ∑ s h is σ ( t ) G s j σ ( t , t (cid:48) ) . (23)An analogous relation can be derived for G ◦ Σ .The result can be written in a more compact form by defin-ing a new quantity η xy σ ( t ) via h i j σ ( t ) = U ( t ) η i j σ ( t ) + T i j σ ( t ) − δ i j µ , (24) h is σ ( t ) = U ( t ) η is σ ( t ) , and η ss (cid:48) σ ( t ) ≡
0. This is consistent with the alternative defi-nition given in Appendix B which also holds for U ( t ) =
0. Inthe physical sector it implies η i j σ ( t ) = δ i j (cid:104) ˆ n i ¯ σ ( t ) (cid:105) H T , U , (25)as follows from U ( t ) η i j σ ( t ) = Σ HF i j σ ( t ) [cf. Eqs. (19), (21)and (24)]. With this definition for η and with the relations Σ ◦ G = − iU ◦ G ( l ) and G ◦ Σ = − iG ( r ) ◦ U we get G ( l ) i j σ ( t , t (cid:48) ) = i ∑ x η ix σ ( t ) G x j σ ( t , t (cid:48) ) , (26) G ( r ) i j σ ( t , t (cid:48) ) = i ∑ x G ix σ ( t , t (cid:48) ) η ∗ ix σ ( t (cid:48) ) , (27)Recall at this point that quantities like h ( t ) = h T , U ( t ) or η ( t ) = η T , U ( t ) (as well as H eff ( t ) , G ( l ) , etc.) are function-als of T and U . B. Hamiltonian-based formulation of the CPT
Let us now discuss how the CPT Green’s function can beobtained from an effective one-particle Hamiltonian and howto set up a Markovian time-propagation scheme. As dis-cussed in Sec. III A, we have T = T (cid:48) + V where T (cid:48) − λ is therenormalized intra-cluster and V + λ the renormalized inter-cluster hopping. For each set of parameters λ , an effectiveone-particle CPT Hamiltonian can be defined by adding theinter-cluster hopping to the effective Hamiltonian (19) of thereference system: H CPT [ λ ]( t ) = H eff T (cid:48) − λ , U ( t ) + ∑ i j σ [ V i j σ ( t ) + λ i j σ ( t )] c † i σ c j σ ≡ ∑ xy σ h CPT xy σ ( t ) c † x σ c y σ . (28)The CPT Green’s function, as computed from H CPT [ λ ]( t ) , G CPT [ λ ] xy σ ( t , t (cid:48) ) = − i (cid:104)T C ˆ c x σ ( t ) ˆ c † y σ ( t (cid:48) ) (cid:105) H CPT [ λ ] (29) then coincides with the original definition in Eq. (14) if onlythe physical sector is considered, i.e., ( x , y ) = ( i , j ) . This canbe verified easily by integrating out the virtual, s degrees offreedom from H CPT . Eq. (28) reflects the freedom we havein the CPT construction as the λ -terms cancel in the phys-ical sector. λ only enters implicitly through the hybridiza-tion strengths h (cid:48) is σ ( t ) , through the on-site energies h (cid:48) ss σ (where h (cid:48) ≡ h T (cid:48) − λ , U ) and through the Hartree-Fock term Σ HF T (cid:48) − λ , U ofthe reference system’s Hamiltonian H T (cid:48) − λ , U .For each set of parameters λ , the two-particle correlationfunction G ( l ) is approximated within the context of the CPTas G ( l ) [ λ ] i j σ ( t , t (cid:48) ) = i ∑ x η (cid:48) [ λ ] ix σ ( t ) G CPT [ λ ] x j σ ( t , t (cid:48) ) , (30)where we have defined η (cid:48) [ λ ] ≡ η T (cid:48) − λ , U . (31)Eq. (30) corresponds to the exact expression given by Eq. (26). G ( r ) [ λ ] is defined analogously, and thus the symmetry rela-tion G ( r ) [ λ ] ji σ ( t , t + ) = (cid:104) G ( l ) [ λ ] i j σ ( t , t + ) (cid:105) ∗ (32)is ensured within the CPT independently of λ . Note that thissymmetry is not sufficient to allow for an unambiguous defi-nition of the doublon density based on G ( l ) and G ( r ) [cf. Eq.(11)]. Instead, it requires both constraints to be respected asdiscussed in Sec. III B. In case of an arbitrary, non-conservingset of parameters λ this ambiguity needs to be circumventedby defining the doublon density as an average d i [ λ ]( t ) = − ∑ σ (cid:104) G ( l ) [ λ ] ii σ ( t , t + ) + G ( l ) [ λ ] ii σ ( t , t + ) (cid:105) . (33)For λ = λ opt , however, we have d i [ λ opt ]( t ) = − G ( l ) [ λ opt ] ii σ ( t , t + ) = − G ( r ) [ λ opt ] ii σ (cid:48) ( t , t + ) . (34)The final forms of the conditional equations for λ opt are ob-tained by replacing G ( r ) and G ( l ) by their CPT approxima-tions G ( l ) [ λ ] and G ( r ) [ λ ] in the expressions for F and C given by Eqs. (15) and (16): F [ λ opt ] i σ ( t ) ! = , C [ λ opt ] i σ ( t ) ! = . (35)We note that the number of free parameters λ must be chosento match the number of linear independent constraints definedby Eq. (35) to ensure the existence of a unique solution λ opt . V. SOLVING THE SELF-CONSISTENCY EQUATIONS
Having formulated the self-consistency conditions (35), itremains to explicitly solve these equations for λ opt . An im-portant simplification arises from the fact that the CPT is byconstruction a fully causal theory, i.e., the time-local elements G CPT ( t , t + ) of the CPT Green’s function at time t , for exam-ple, only depend on quantities at earlier times. The same holdsfor G ( l ) ( t , t + ) and for h CPT xy σ ( t ) . This allows us to construct astrategy for the solution of Eq. (35) in the form of a time-propagation algorithm. Let us therefore assume that λ opt isknown for all time points on a discrete time grid and thatonly the parameters λ opt ( t ) at the latest point of time t areunknown.Therewith, the actual task is to solve Eq. (35) for λ opt ( t ) only at the given latest point of time t . To this end we have toanalyze at time t the λ ( t ) dependence of the relevant quanti-ties, i.e., of G ( l ) ( t , t + ) and G ( r ) ( t , t + ) , see Eqs. (15) and (16).First of all, the dependence of G ( l ) ( t , t + ) (and G ( r ) ( t , t + ) )on λ ( t ) at time t is due to the CPT Hamiltonian h CPT xy σ ( t ) [seeEq. (28) and see Eqs. (29) and (30)]. The λ ( t ) -dependenceof the latter is exclusively due to the time-evolution operator U (cid:48) [ λ ] ≡ U T (cid:48) − λ , U of the reference system. The detailed con-struction of h CPT xy σ ( t ) is not important here, and we refer to Ref.39 for a comprehensive discussion. Finally, the functional de-pendence of U (cid:48) [ λ ]( t , ) on λ is through an integration over alltimes between 0 and t . With this information at hand, we arein fact able to characterize the dependence on λ ( t ) at time t ofthe quantities G ( l ) ( t , t + ) and G ( r ) ( t , t + ) which enter the localconstraints (35) that serve to enforce the conservation laws.The most important point for the following discussion isthe fact that, in the limit of vanishing time step ∆ t →
0, theparameter set λ ( t ) at the latest point of time enters basicallyall central quantities as a null set only: Consider, for exam-ple, G ( l ) ( t , t + ) . Its first-order response due to a variation of λ ( t ) at time t vanishes (as shown below). On the one hand,this missing sensitivity implies a complication of the theorysince one has to account for this mathematical property ex-plicitly when setting up a numerical implementation. On theother hand, once one has recognized the property, it actu-ally helps to the solve Eqs. (35). Consider a given arbitrarycausal functional M [ λ ]( t ) . The main trick is to enhance thesensitivity of M [ λ ]( t ) to variations of λ ( t ) at time t by tak-ing its time derivative. Typically, if the first-order responseof M [ λ ]( t ) vanishes, ∂ t M [ λ ]( t ) is a linear function of λ ( t ) at time t . Clearly, this is the key to solve an equation like M [ λ ]( t ) = λ opt ( t ) .In the following subsections V A – V D the above-sketchedideas are worked out on a more technical level. Finally, thesection V E addresses the initial state at time t = A. Time-local variations
Assume that we have found the optimal renormalization λ opt ( t ) for t ≤ t n ≡ n ∆ t . We introduce a variation δ n loc whichaffects the current (the n -th) time step only: δ n loc λ i j σ ( t ) = δ λ i j σ ( t ) Θ n loc ( t ) , Θ n loc ( t ) = (cid:26) t ∈ [ t n , t n + ] , . (36)For simplicity, we require the variations to be symmetric, i.e., δ λ i j σ ( t ) = δ λ ji σ ( t ) . This implies a restriction to symmetricsolutions λ opt . Consider now an arbitrary, causal functional M [ λ ]( t ) , i.e., a functional that at time t only depends on λ ( t (cid:48) ) with t (cid:48) ≤ t . For such an object, the variational operator δ n loc isrelated to the conventional functional derivative through δ n loc M [ λ ]( t ) = ∑ σ ∑ i ≥ j (cid:90) tt n d t (cid:48) δ M [ λ ]( t ) δ λ i j σ ( t (cid:48) ) δ λ i j σ ( t (cid:48) ) , (37)where the restriction i ≥ j is necessary because of the symme-try requirement λ i j σ = λ ji σ .We now take the combined limit n → ∞ , ∆ t → t ∈ [ t n , t n + ] to define the time-local variation δ loc in the continuum limit δ loc M [ λ ]( t ) = lim ∆ t → n → ∞ δ n loc M [ λ ]( t ) , (38)with the corresponding variational quotient δ loc M [ λ ]( t ) δ loc λ i j σ ( t ) ≡ lim ∆ t → n → ∞ (cid:90) tt n d t (cid:48) δ M [ λ ]( t ) δ λ i j σ ( t (cid:48) ) . (39)This variational quotient describes the linear response of M [ λ ]( t ) when varying the parameters at the latest time step: δ loc M [ λ ]( t ) = ∑ σ ∑ i ≥ j δ loc M [ λ ]( t ) δ loc λ i j σ ( t ) δ λ i j σ ( t ) . (40) B. Integrated quantities in λ With the appropriate variation for our purposes at hand,we can study the effect of the variation on the main quanti-ties within the CPT framework. We first consider the time-evolution operator (“propagator”) of the reference system U (cid:48) [ λ ] ≡ U T (cid:48) − λ , U . It is instructive to study the effect of theoperator δ n loc first, i.e., the effect of a time-local variation withfinite time step ∆ t . Keeping only terms of the order O ( ∆ t ) onefinds δ n loc U (cid:48) [ λ ]( t , ) = − i (cid:34) ∑ i j σ (cid:90) tt n δ λ i j σ ( t (cid:48) ) ˆ c † i σ ˆ c j σ d t (cid:48) (cid:35) U (cid:48) [ λ ]( t n , )+ O ( ∆ t ) . (41)In lowest order we thus have δ n loc U T (cid:48) − λ , U ( t , ) ∝ ∆ t δ λ ( t ) .This means that the linear response vanishes identically in thelimit ∆ t →
0. This property originates from the fact that λ ( t ) is integrated over time within the propagator U T (cid:48) − λ , U ( t , ) ,and that the contribution of a single time step, t ∈ [ t n , t n + ] ,to this integral is of zero measure in the limit ∆ t → δ loc [ i ∂ t U (cid:48) [ λ ]( t , )] δ loc λ i j σ ( t ) = − (cid:2) c † j σ c i σ + c † i σ c j σ − δ i j c † i σ c i σ (cid:3) U (cid:48) [ λ ]( t , ) . (42)Multiplying this equation with λ i j σ ( t ) , summing over i , j , σ and comparing with the standard equation of motion i ∂ t U (cid:48) [ λ ]( t , ) = H T (cid:48) − λ , U ( t ) U (cid:48) [ λ ]( t , ) , shows that the timederivative of the propagator is of the general form i ∂ t U (cid:48) [ λ ]( t , ) = ∑ σ ∑ i ≥ j δ loc [ i ∂ t U (cid:48) [ λ ]( t , )] δ loc λ i j σ ( t ) λ i j σ ( t ) + ξ U (cid:48) [ λ ]( t ) , (43)where ξ U (cid:48) [ λ ]( t ) = H T (cid:48) , U ( t ) U (cid:48) [ λ ]( t , ) . Note that the depen-dence on λ i j σ ( t ) at time t is strictly linear in the limit ∆ t → ξ U (cid:48) [ λ ]( t ) on the right-hand side ofEq. (43) depend on λ ( t ) only through an integration over timewithin the propagator U (cid:48) [ λ ]( t , ) . We will call such quantities integrated quantities in λ . Integrated quantities in λ inherit animportant property from the cluster propagator U (cid:48) [ λ ] , see Eq.(41): Their time-local variation vanishes in the limit ∆ t → M [ λ ] thatis integrated in λ , i.e., the time derivative of a functional of theform M [ λ ]( t ) = M ( U (cid:48) [ λ ]( t , )) , can be brought into a formanalogous to Eq. (43). This follows immediately from thechain rule in calculus as i ∂ t M [ λ ]( t ) = i ∂ M ( U (cid:48) ) ∂ U (cid:48) ∂ U (cid:48) [ λ ]( t , ) ∂ t . Ex-plicitly this result reads i ∂ t M [ λ ]( t ) = ∑ σ ∑ i ≥ j δ loc [ i ∂ t M [ λ ]( t )] δ loc λ i j σ ( t ) λ i j σ ( t ) + ξ M [ λ ]( t ) , (44)where δ loc [ i ∂ t M [ λ ]( t )] δ loc λ ij σ ( t ) and ξ M [ λ ] are again integrated quantitiesin λ . We furthermore conclude that a time-local variation ofthe time derivative of an integrated quantity in λ is non-zeroin general.The main idea in the following is to combine the conditionalequations (35) into a single equation Γ [ λ opt ]( t ) ! = Γ [ λ ] is of the form Γ [ λ ]( t ) = J [ λ ]( t ) λ ( t ) + ξ Γ [ λ ]( t ) where J [ λ ] and ξ Γ [ λ ] are integrated quantities in λ . This is formallyeasily solved for λ opt ( t ) by matrix inversion and allows to de-rive an efficient propagation scheme for numerical purposes. C. λ -dependence of G ( l ) and G ( r ) The main building blocks of the local constraints on thespin-dependent density, Eq. (15), and the doublon density, Eq.(16), are given by the two-particle correlation functions G ( l ) and G ( r ) . Within the CPT approximation they are definedthrough Eq. (30). We therefore have to understand the λ de-pendence of η (cid:48) [ λ ] ≡ η T (cid:48) − λ , U and G CPT [ λ ] .One can easily see that η (cid:48) [ λ ] is an integrated quantity in λ .Consider, for example, the physical sector. From Eq. (25) wehave η (cid:48) i j σ [ λ ]( t ) = δ i j (cid:104) ˆ n i ¯ σ ( t ) (cid:105) H T (cid:48)− λ , U . The only λ -dependenceof this expression indeed stems from the propagator U (cid:48) [ λ ] .To obtain a non-vanishing time-local variation we thus haveto consider the first derivative with respect to time. This isworked out in Appendix C: δ loc [ i ∂ t η (cid:48) i j σ ( t )] = η (cid:48) ii σ ( t ) δ λ i j σ ( t ) − ∑ l σ (cid:48) [ δ λ il σ (cid:48) ( t )] γ l σ (cid:48) i j σ ( t ) , δ loc [ i ∂ t η (cid:48) is σ ( t )] = − ∑ l σ (cid:48) [ δ λ il σ (cid:48) ( t )] γ l σ (cid:48) is σ ( t ) , (45) where the newly introduced tensor γ [ λ ] l σ (cid:48) is σ ( t ) is cluster-diagonal, i.e., γ [ λ ] l σ (cid:48) is σ ( t ) (cid:54) = i and l refer tolattice sites within the same cluster. It furthermore followsthat i ∂ t η (cid:48) [ λ ]( t ) can be brought into the form specified by Eq.(44), where the variational derivative δ loc [ i ∂ t η (cid:48) [ λ ] ix σ ( t )] δλ jl σ ( t ) , as givenby Eq. (45), and ξ η (cid:48) [ λ ] ix σ ( t ) are integrated quantities in λ .An explicit expression for the latter is not needed for our pur-poses.Let us now take a look at the CPT Green’s func-tion. It depends on λ through the Hamiltonian H CPT [ λ ] ,which in turn depends on λ through the hybridizationstrengths h (cid:48) [ λ ] is σ ( t ) = U ( t ) η (cid:48) [ λ ] is σ ( t ) and the Hartree-Fock term Σ (cid:48) [ λ ] HF i j σ ( t ) = U ( t ) η (cid:48) [ λ ] i j σ ( t ) . The Hamilto-nian H CPT ( t ) is therefore an integrated quantity in λ . Asthe propagator U CPT [ λ ]( t , ) = T exp (cid:0) − i (cid:82) t d t (cid:48) H CPT [ λ ]( t (cid:48) ) (cid:1) involves a second integral over time, we conclude that δ n loc G CPT ( t , t + ) ∝ ∆ t δ λ ( t ) . In this sense, G CPT must be seenas an integrated quantity in λ of second order. Consequently,the time-local variation of its first derivative with respect totime vanishes: δ loc [ i ∂ t G CPT [ λ ]( t , t + )] = . (46)We note that the time derivative involves the product of thematrix elements of the CPT Hamiltonian, Eq. (28), with G CPT [ λ ]( t , t + ) , i.e., the product of an integrated quantity in λ with an integrated quantity in λ of second order, respec-tively. Obviously, the product scales like an ordinary inte-grated quantity in λ when a time-local variation is applied,i.e., δ n loc h CPT ( t ) G CPT ( t , t + ) ∝ ∆ t δ λ ( t ) in lowest order in ∆ t .Concluding, to get a non-vanishing time-local variation,one must consider the first time derivative of the two-particleGreen’s functions G ( l ) and G ( r ) . We find δ loc (cid:104) i ∂ t G ( l ) i j σ ( t , t + ) (cid:105) = ∑ x (cid:0) δ loc [ i ∂ t η (cid:48) ix σ ( t )] (cid:1) G CPT x j σ ( t , t + ) (47)and an analogous expression for G ( r ) . Only the η (cid:48) -term con-tributes to the variation, cf. Eq. (45), while the variation of theCPT Green’s function vanishes, cf. Eq. (46). We also note thatEq. (45) may be used at this point and that i ∂ t G ( l ) ( t , t + ) [andanalogously i ∂ t G ( r ) ( t , t + ) ] is of the form i ∂ t G ( l ) i j σ ( t , t + ) = ∑ σ (cid:48) ∑ k ≥ l δ loc G ( l ) i j σ ( t , t + ) δ loc λ kl σ (cid:48) ( t ) λ kl σ (cid:48) ( t ) + [ ξ G ( l ) ] i j σ ( t ) , (48)where both, the variational derivative δ loc G ( l ) ( t , t + ) δ loc λ ( t ) , as givenby Eq. (47), and the quantity ξ G ( l ) ( t ) , which is not neededin explicit form for our purposes, scale like integrated quan-tities in λ under time-local variations. This follows from thefact that i ∂ t G CPT ( t , t + ) scales like an integrated quantity in λ under time-local variations and the related discussion above. D. λ -dependence of the local constraints The local constraint on the spin-dependent density, Eq.(15), is formulated in terms of the difference between G ( l ) and G ( r ) . Therefore its first derivative with respect to timemust be considered to obtain a non-vanishing time-local vari-ation: δ loc [ i ∂ t F [ λ ] i σ ( t )] = δ loc (cid:104) i ∂ t G ( l ) [ λ ] ii σ ( t , t + ) (49) − i ∂ t G ( r ) [ λ ] ii σ ( t , t + ) (cid:105) , For the time-local variation of the local constraint on the dou-blon density, Eq. (16), on the other hand, one finds δ loc C [ λ ] i σ ( t ) = δ loc (cid:104) i ∂ t G ( l ) [ λ ] ii σ ( t , t + ) (50) + i ∂ t G ( r ) [ λ ] ii σ ( t , t + ) (cid:105) , since δ loc [ T ◦ G ( r ) − G ( l ) ◦ T ]( t , t + ) =
0, where we made useof the fact that T = T (cid:48) + V is the hopping of the original sys-tem and thus independent of λ .To treat both constraints in a combined formal frame, wedefine the functional Γ [ λ ] : Γ [ λ ] a ( t ) = (cid:26) i ∂ t F [ λ ] i σ ( t ) if 0 ≤ a < L , C [ λ ] i σ ( t ) if 2 L ≤ a < L , (51)where L is the number of lattice sites. With this, theconditional equation for the optimal renormalization reads Γ [ λ opt ] ! =
0. From the previous discussion and Eq. (48) it fol-lows that Γ [ λ ] a ( t ) is of the form Γ [ λ ] a ( t ) = ∑ b J [ λ ] ab ( t ) λ b ( t ) + ξ Γ [ λ ] a ( t ) , (52)where we introduced the super-index b which labels the set offree parameters: λ b = λ i j σ , i ≥ j . Both J [ λ ] and ξ Γ [ λ ] scalelike integrated quantities in λ under time-local variations. TheJacobian matrix J is defined as J [ λ ] ab ( t ) ≡ δ loc Γ a [ λ ]( t ) δ loc λ b ( t ) . (53)The matrix J [ λ ]( t ) is quadratic if the number of free parame-ters λ b is chosen such that it equals the number of conditionalequations [see Eq. (51)]. Assuming that J [ λ ]( t ) is regular,one can formally solve the conditional equation for the opti-mal renormalization: Γ [ λ opt ]( t ) ! = ⇔ λ opt ( t ) = − (cid:2) J [ λ opt ]( t ) (cid:3) − ξ Γ [ λ opt ]( t ) , (54)see Eq. (52). This completes our derivation.Let us emphasize that the single point λ opt ( t ) represents anull set with respect to the time-integrations in J [ λ opt ]( t ) and ξ Γ [ λ opt ]( t ) . This can be exploited to derive an efficient numer-ical scheme to obtain λ opt ( t ) step by step on the time axis asdetailed in Appendix D. There we also argue why finding anexplicit expression for ξ Γ [ λ ]( t ) can in fact be circumvented.An explicit expression for J [ λ ]( t ) in terms of known quanti-ties, on the other hand, is available via Eqs. (45), (47), (49)and (50). E. The equilibrium initial state
Initially, at time t = λ eq ≡ λ ( ) for the initial state, and thus a nontriv-ial self-consistency condition is not available, unfortunately.This can be seen as follows: Let us assume that the hoppingmatrix T , and consequently T (cid:48) and V , are real and symmetric.Consider G ( l ) at times t = t (cid:48) =
0. Via Eq. (30) this is givenas G ( l ) [ λ ] i j σ ( , + ) = i ∑ x η (cid:48) [ λ ] ix σ ( ) G CPT [ λ ] x j σ ( , + ) . InAppendix B it is shown that at time t = { η (cid:48) [ λ ] ix σ ( ) } vanishes, independently of λ eq . Hence,Eq. (24) implies that H CPT [ λ ]( ) is real and symmetric,and therefore G CPT [ λ ] xy σ ( , + ) = i (cid:104) ˆ c † y σ ( ) ˆ c x σ ( ) (cid:105) H CPT [ λ ] ispurely imaginary. Consequently, G ( l ) [ λ ]( , + ) is real. Fi-nally, we conclude with Eq. (32) that G ( l ) [ λ ] i j σ ( , + ) = G ( r ) [ λ ] ji σ ( , + ) . (55)This directly proves that F [ λ ] i σ ( ) =
0, and furthermore C [ λ ] i σ ( ) = ∑ j σ T i j σ ( ) (cid:104) G ( r ) [ λ ] ji σ ( , + ) (56) − G ( l ) [ λ ] i j σ ( , + ) (cid:105) = λ eq , i.e., both constraints hold trivially.For the concrete numerical calculations we therefore cir-cumvent this issue and consider a noninteracting initial state.In this case the CPT is exact, independent of the choice of λ .The initial value λ eq is then fixed by requiring λ to be contin-uous so that λ eq = λ ( + ) . VI. NUMERICAL RESULTS
The conserving CPT has been implemented numerically.First results are discussed for the two-dimensional Hubbardmodel on an L = ×
10 square lattice with periodic boundaryconditions. As these results shall serve as a proof of conceptonly, we restrict ourselves to the most simple approximation,i.e., to the smallest meaningful cluster as the building blockof the reference system, namely a cluster consisting of 2 × µ =
0. Note that the CPT de-scription of this initial state is exact (and independent of therenormalization). The hopping of the original system is re-stricted to nearest neighbors, and we set the nearest-neighborhopping T = FIG. 3: Partitioning of the Hubbard model with nearest-neighborhopping on the two-dimensional square lattice used for the numeri-cal calculations. Original system: L = ×
10 lattice with periodicboundary conditions. Reference system: 50 clusters of size 2 × × T and optimization parameter λ ( t ) are indicated by black and red lines. The time-dependent renor-malization λ ( t ) is employed to enforce the conservation of energy inthe real-time dynamics following an interaction quench. where the Hubbard- U is suddenly, at time t =
0, switched onto a finite value U fin : U ( t ) = Θ ( t ) U fin . (57)Here, Θ ( t ) denotes the Heaviside step function. For times t > µ = µ = U / For a spin-independentparameter quench, as considered here, the CPT also triviallyrespects the conservation of the total spin. Total-energy con-servation, on the other hand, is violated in a conventionalCPT approach as has been explicitly demonstrated recently. For the present setup we will therefore employ the nearest-neighbor hopping within the 2 × λ ( t ) (see Fig. 3).We note that the computational effort to self-consistentlyevaluate the presented theory numerically is essentially deter-mined by the underlying solver for the conventional nonequi-librium CPT with little overhead. Here, we use the time-local,Hamiltonian-based solver developed in Ref. 39 which con-structs the effective Hamiltonian of each cluster by exact di-agonalization. For the 2 × ∆ t = .
001 to propagate the sys-tem up to 26 ,
500 steps up to a maximum propagation time of t max = .
5. For each such time step t n , the scheme developedin Appendix D has been employed with n p =
1, i.e., we havecalculated the Taylor coefficients λ n , and λ n , .While the required computational resources are very mod-erate, accessing longer time scales has turned out to be hin-dered by mathematical complications. As is obvious from Eq. t . . . . . . d ( t ) CPT (plain)CPT (conserving)CPT (plain)CPT (conserving)
FIG. 4: Time evolution of the local doublon density after an in-teraction quench U = → U fin . Grey lines: plain CPT. Blue lines:conserving CPT. Results for different U fin ranging from U fin = . U fin = . ∆ U fin = . (54), an inversion of the Jacobian matrix J ( t ) is necessary toobtain λ opt ( t ) at each time step. However, with increasing U fin this matrix exhibits singular points of non-invertibilityat earlier and earlier times. In fact, one finds numericallythat also the starting point t = + is singular, namely the Ja-cobian matrix vanishes: J ( + ) =
0. Fortunately, one alsohas ξ Γ ( + ) =
0, such that this problem is fixed by applyingL’Hˆopital’s rule. At time t = + , the defining equation for λ opt ( + ) becomes λ opt ( + ) = − (cid:2) ∂ t J [ λ opt ]( + ) (cid:3) − [ ∂ t ξ Γ [ λ opt ]( + )] . (58)While this solves the problem at time t = + , finding a sys-tematic and convenient way to propagate beyond the singularpoints of the Jacobian matrix at finite times remains topic forfuture investigations.Apart from this technical problem, the suggested schemeworks as expected. Results for the time evolution of thedoublon density are shown in Fig. 4. It is evident that therenormalization λ has a strong influence on the dynamicsand leads to qualitatively different results when comparingthe plain unoptimized CPT calculation with the novel con-serving CPT. While the dynamics is characterized by ongo-ing oscillations when using plain CPT, there is a monotonicdecay of the doublon density in case of the conserving CPT.The longest maximum propagation time is achieved for thequench U = → .
5. Here, the first singular point of the Ja-cobian shows up at t max ≈ .
5. On this time scale, the dou-blon density seems to relax to a stationary state with little tono oscillations.The qualitatively different time dependence of the doublondensity reflects the qualitatively different behavior found forthe total energy in the plain and the conserving CPT: This isshown in the inset of Fig. 5. For the conserving CPT, the to-tal energy is perfectly conserved within numerical accuracy –by construction of the approach. In the plain CPT calculation,however, the total energy shows unphysical oscillations. Here,maxima and minima of E tot ( t ) nicely correspond to maximaand minima in the plain-CPT doublon density seen in Fig. 4.1 t − . − . . . . λ ( t ) . . . . . . . . . . . . t − . − . − . − . − . E t o t ( t ) U = → . CPT (plain)CPT (conserving)CPT (plain)CPT (conserving)
FIG. 5: Time evolution of the optimal renormalization parameter λ opt ( t ) for different U fin as indicated and corresponding to Fig. 4. In-set: time dependence of the total energy (plain and conserving CPT). It must be concluded that those are artifacts of the plain CPTapproach. We also note that small unphysical oscillations ofthe total energy density E tot ( t ) / L (with L = .
01 lead to much stronger oscillations in thedoublon density with amplitudes of about 0.04.The main part of Fig. 5 displays the results for the time evo-lution of the renormalization parameter λ ( t ) . Its dependenceon U fin turns out to be rather weak on a time scale of a fewinverse hoppings. Irrespective to the final interaction strength U fin , the initial equilibrium value is found as λ eq ≈ − .
86. For t > U fin , the renormalization parameter rapidly in-creases to λ ≈ . t ≈ .
7. This corre-sponds to the rapid initial drop of the doublon density (cf. Fig.4). Results for longer times are again only available for thequench U = → .
5. On the time scale up to t max ≈ .
5, weobserve a subsequent slow relaxation of λ ( t ) toward an aver-age final value λ ∞ ≈ . U fin .One should note that λ ∞ = T (cid:48) − λ ∞ =
0, i.e.,a vanishing renormalized intra-cluster hopping. Apart fromthe remaining oscillations of the renormalization parameteraround λ ∞ =
1, this means that the system “chooses” theatomic limit of the Hubbard model as the optimal startingpoint for the all-order perturbation theory in the inter-clusterhopping for long times. This may be interpreted as follows:First of all, the remaining oscillations are understood as beingnecessary to keep the total energy constant within the con-serving CPT. Disregarding the oscillations, the value λ ∞ = λ ( t ) = λ ∞ = using the density-matrix renormalization group (DMRG) and for the model in infinite dimensions usingthe dynamical mean-field theory (DMFT). In both cases, avery fast relaxation of the doublon density on a time scale ofone inverse hopping has been found in fact. Typically, how-ever, the doublon density first develops a minimum before itsaturates to an almost constant value. This minimum is absentin the conserving CPT calculations (see Fig. 4). Note, that forweak quenches and on the intermediate time scale discussedhere and in the DMRG and DMFT studies, the doublon den-sity does not relax to its thermal value due to kinematic con-straints becoming active after the ultrashort initial relaxationstep.
Indeed, one expects a subsequent relaxation on amuch longer time scale. Let us emphasize that while our datain Fig. 4 are compatible with these expectations, serious pre-dictions using the conserving CPT are not yet possible. Thiswould require a much more systematic study involving differ-ent and in particular larger clusters, an analysis of the depen-dence on the cluster shape and also a systematic discussion ofthe different possibilities to choose renormalization parame-ters for the self-consistent procedure.
VII. CONCLUSIONS
Cluster-perturbation theory, as proposed originally, repre-sents the most simple way to construct a mean-field theorywhich incorporates to some extent the effects of short-rangecorrelations. We have emphasized that the starting point ofthe perturbational expansion in the inter-cluster hopping is byno means predetermined and that the according freedom inthe choice of the intra-cluster hopping parameters can be ex-ploited to “optimize” the mean-field theory. There are dif-ferent conceivable optimization schemes. One way is to addan additional self-consistency condition as, for example, aself-consistent renormalization of the on-site energies whichwould be very much in the spirit of the Hubbard-I approxima-tion. The disadvantage of such ideas is their ad hoc charac-ter. An optimization following a general variational principleis much more satisfying and physically appealing. This is theroute that is followed up by self-energy-functional theory. Un-fortunately, total-energy conservation is not straightforwardlyimplemented within the SFT context. An appealing idea isthus to use the above-mentioned freedom to enforce energyconservation, and actually any conservation law dictated bythe symmetries of the problem at hand. This leads to the con-serving CPT proposed in the present paper.As we have argued (see Sec. V E), this idea can exclusivelybe used to constrain the CPT real-time dynamics while otherconcepts must be invoked for the initial thermal state. On theother hand, there is an urgent need for numerical approaches,even for comparatively simple cluster mean-field concepts,which are able to address the real-time dynamics of stronglycorrelated lattice fermion models beyond the more simple ex-treme limits of one and infinite lattice dimension.With the present paper we could give a proof of principlethat a nonequilibrium conserving cluster-perturbation theoryis possible and can be evaluated numerically in practice. Anhighly attractive feature of this approach is the linear scaling2with the propagation time, while the exponential scaling withthe cluster size is the typical bottleneck of any cluster mean-field theory.The mapping of the original nonequilibrium CPT ontoan effective auxiliary problem specified by a noninteractingHamiltonian with additional virtual (“bath”) degrees of free-dom is crucial for the practical implementation of the ap-proach. One should note that the number of virtual sites is re-lated to the number of one-particle excitations and thus growsexponentially with the original cluster size. Hence, any practi-cal calculation is limited to a few (say, at most 10) cluster sitesonly. This implies that a systematic finite-size scaling anal-ysis will be problematic if long-range correlations dominatethe essential physical properties – this is the above-mentioneddrawback that is shared with any available cluster-mean-fieldtheory. We therefore expect that the field of applications ofconserving CPT is limited to problems with possibly strongbut short-ranged correlations.Due to its formulation in terms of Green’s functions withtime arguments on the Keldysh contour, the CPT has an in-herently causal structure. With the present paper we could inparticular demonstrate how to exploit this causality for an effi-cient time-stepping algorithm where updates of the parameterrenormalization can be limited to the respective last time sliceduring time propagation. The essential problem that had to besolved here consists in controlling the order (in the sense ofa Taylor series) at which the parameter renormalization on asingle time slice enters other physical quantities, such as thebasic time-evolution operator, Green’s functions, etc. This hasallowed us to set up a highly accurate numerical algorithmwhere conservation laws are respected with machine preci-sion.For convenience, first numerical results have been gener-ated for interaction quenches of the two-dimensional Hub-bard model on a square lattice at half-filling, where particle-number and spin conservation are respected trivially. Energyconservation has been enforced by time-dependent renormal-ization of the intra-cluster hopping in the 2 × L c sites, at most 4 L c parameters are needed. In addition, boththe number of constraints and the optimization parameters de-pend on the spatial symmetries and other symmetries, e.g.,particle-hole symmetry, of the original and the reference sys-tem. If necessary, more degrees of freedom and correspond-ingly more parameters can be generated by coupling uncorre-lated “bath” sites to the physical sites in the reference systemin the spirit of (cluster) dynamical mean-field theory. Sys-tematic studies addressing the mentioned issues are necessarybefore a systematic and quantitative comparison with otherapproaches or with experiments is meaningful.Interestingly, the conditional equations for the renormal-ization parameters feature singular points of non-invertibility.Technically, this currently restricts our investigations toquenches with small final interaction and short propagationtimes. It is not clear at the moment whether or not a physicalmeaning can be attributed to those singular points; also thisrequires further systematic studies. According to our presentexperience, it is well conceivable that, with a suitable regu-larization scheme, time propagation through a singularity ofthe Jacobian is possible and has no apparent impact on thetime dependence of physical observables. Developing such aregularization scheme is the next task for future studies andthe most important issue to make the conserving CPT a pow-erful numerical tool to address, e.g., real-time magnetizationdynamics, even of inhomogeneous models and on long timescales. Acknowledgments
We thank Roman Rausch for providing an exact-diagonalization solver for the Hubbard model and Felix Hof-mann for valuable discussions. This work has been supportedby the Deutsche Forschungsgemeinschaft through the excel-lence cluster “The Hamburg Centre for Ultrafast Imaging -Structure, Dynamics and Control of Matter at the AtomicScale” and through the Sonderforschungsbereich 925 (projectB5). Numerical calculations were performed on the PHYSnetcomputer cluster at the University of Hamburg.
Appendix A: Local constraint on the doublon density
Within this subsection we will use the shorthand notation G ( l ) ≡ G ( l ) T , U and G ( r ) ≡ G ( r ) T , U . To prove the local constrainton the doublon density, we consider − i ∂ t d i ( t ) = (cid:104) [ ˆ H T , U ( t ) , ˆ n i ↑ ( t ) ˆ n i ↓ ( t )] (cid:105) H T , U = ∑ σ (cid:104) (cid:34) ∑ jk T jk σ ( t ) ˆ c † j σ ( t ) ˆ c k σ ( t ) , ˆ n i σ ( t ) (cid:35) ˆ n i ¯ σ ( t ) (cid:105) H T , U , (A1)where we used that the double occupation operator commuteswith the interaction term of the Hamiltonian H T , U . Using fur-3ther that ∑ jk T jk σ ( t ) (cid:104) ˆ c † j σ ( t ) ˆ c k σ ( t ) , ˆ n i σ ( t ) (cid:105) (A2) = ∑ jk T jk σ ( t ) (cid:16) δ ki ˆ c † j σ ( t ) ˆ c i σ ( t ) − δ ji ˆ c † i σ ( t ) ˆ c k σ ( t ) (cid:17) , we find the final form by comparing with Eqs. (8) and (10) andusing the relation d i ( t ) = − G ( l ) ii σ ( t , t + ) = − G ( r ) ii σ (cid:48) ( t , t + ) . Thisimplies − i ∂ t d i ( t ) = i ∂ t (cid:104) G ( l ) ii σ ( t , t + ) + G ( r ) ii σ ( t , t + ) (cid:105) (A3) = ∑ j σ (cid:104) T i j σ ( t ) G ( r ) ji σ ( t , t + ) − G ( l ) i j σ ( t , t + ) T ji σ ( t ) (cid:105) , which completes our derivation of Eq. (16).To prove that Eq. (A3) indeed ensures energy conserva-tion, let us consider a time-independent Hamiltonian with T i j σ ( t ) = T i j σ and U ( t ) = U . We consider the time-derivativeof the kinetic energy first. Since the kinetic part of the Hamil-tonian trivially commutes with itself, one obtains i ∂ t E kin ( t ) = (cid:104) ∑ i j σ T i j σ ˆ c † i σ ˆ c j σ , U ∑ l n l ↑ n l ↓ (cid:105) (A4) = U ∑ i j σ (cid:104) T i j σ G ( r ) ji σ ( t , t + ) − G ( l ) i j σ ( t , t + ) T ji σ (cid:105) , This term cancels with i ∂ t E int ( t ) = U ∑ i i ∂ t d i ( t ) assumingEq. (A3) holds thus proving energy conservation for a time-independent Hamiltonian. Appendix B: Mathematical structure of the effectiveHamiltonian
In Eq. (19) we have stated the effective one-particle Hamil-tonian H eff ( t ) ≡ H eff T , U ( t ) for an interacting lattice system. Itscorresponding matrix elements are given by h ( t ) ≡ h T , U ( t ) ,cf. Eq. (21). In Ref. 39 it is shown that h can be constructedfrom the Lehmann representation of the one-particle Green’sfunction G ≡ G T , U . It is given by h xy σ ( t ) = ∑ m , n [ i ∂ t O x σ ( m , n ) ( t )] O ∗ y σ ( m , n ) ( t ) . (B1)The matrix O ( t ) can explicitly be stated within the physicalsector only (i.e., ( x , y ) = ( i , j ) ). There it takes the form O i σ ( m , n ) ( t ) = (cid:114) e − β E m + e − β E n Z (cid:104) m | ˆ c i σ ( t ) | n (cid:105) , (B2)where | m (cid:105) denotes the m -th eigenstate with correspond-ing eigenenergy E m of the initial Hamiltonian, i.e., H T , U ( ) | m (cid:105) = E m | m (cid:105) . Z denotes the grand-canonical partitionfunction, β the inverse temperature of the equilibrium initialstate. The remaining matrix elements O s σ ( m , n ) ( t ) are defineduniquely, up to rotations in invariant subspaces, by requiring O ( t ) to be a unitary transform, see Ref. 39 for details. Eq. (24)is now easily derived from i ∂ t O i σ ( m , n ) ( t ) = ∑ j ( T i j σ ( t ) − δ i j µ ) O j σ ( m , n ) ( t ) (B3) + U ( t ) R i σ ( m , n ) ( t ) , where R i σ ( m , n ) ( t ) = (cid:114) e − β E m + e − β E n Z (cid:104) m | ˆ n i ¯ σ ( t ) ˆ c i σ ( t ) | n (cid:105) . (B4)With the definition η ix σ ( t ) = ∑ mn R i σ ( m , n ) ( t ) O ∗ x σ ( m , n ) ( t ) (B5)we arrive at Eq. (24). It is instructive to note that the elements η i j σ ( t ) in the physical sector are easily evaluated as η i j σ ( t ) = (cid:104) ˆ n i ¯ σ ( t ) (cid:110) ˆ c i σ ( t ) , ˆ c † j σ ( t ) (cid:111) (cid:105) H T , U = δ i j (cid:104) ˆ n i ¯ σ ( t ) (cid:105) H T , U , (B6)where { A , B } = AB + BA denotes the anti-commutator. Froma numerical point of view we note, that η ( t ) as well as itsderivatives with respect to time of arbitrary order can be easilyevaluated using the exact-diagonalization-based scheme pre-sented in Ref. 39.Regarding our discussion of the initial state in Sec. V E itis now easy to proof that Im { η ix σ ( ) } =
0. At time t = R i σ ( m , n ) ( t ) can be chosen as real, since the Hamilto-nian H T , U ( ) is symmetric and therefore has real eigenvectors | m (cid:105) . Similarly, the physical sector O i σ ( m , n ) ( t ) is real in thiscase and completion to a unitary matrix O ( t ) allows for an or-thogonal, i.e., real O ( t ) . It then follows that Im { η ix σ ( ) } = O ( t ) can bebrought into the form ˜ O ( ) = O ( ) S , where the matrix S con-tains the phase factors and possibly rotations in invariant sub-spaces. However, from Eq. (B3) it follows that ˜ R ( ) = R ( ) S and therefore the S matrix cancels. The same proof holds forthe reference system, i.e., Im { η (cid:48) [ λ ] ix σ ( ) } =
0, assuming λ is real. Appendix C: Calculating the time-local variation of η Let h (cid:48) = h T (cid:48) − λ , U denote the matrix elements of the effectiveHamiltonian of the reference system H eff T (cid:48) − λ , U , cf. Eqs. (19)and (21). Through Eq. (24), or equivalently Eq. (B5), a corre-sponding matrix η (cid:48) ≡ η T (cid:48) − λ , U is defined. We are interested inhow it transforms under time-local variations. Since η (cid:48) is anintegrated quantity in λ , we first calculate its derivative withrespect to time. Eq. (B5) implies i ∂ t η (cid:48) ix σ ( t ) = − ∑ y η (cid:48) iy σ ( t ) h (cid:48) yx σ ( t ) (C1) + ∑ mn (cid:104) i ∂ t R (cid:48) i σ ( m , n ) ( t ) (cid:105) [ O (cid:48) ] † ( m , n ) x σ ( t ) , R (cid:48) ≡ R T (cid:48) − λ , U and O (cid:48) ≡ O T (cid:48) − λ , U . With H (cid:48) = H T (cid:48) − λ , U ,the time-local variation of i ∂ t R (cid:48) i σ ( m , n ) ( t ) is given by δ loc (cid:104) i ∂ t R (cid:48) i σ ( m , n ) ( t ) (cid:105) = z ( m , n ) (cid:104) m | δ loc (cid:2) ˆ n i ¯ σ ( t ) ˆ c i σ ( t ) , ˆ H (cid:48) ( t ) (cid:3) | n (cid:105) = − z ( m , n ) ∑ j (cid:104) m | δ λ i j σ ( t ) ˆ n i ¯ σ ( t ) ˆ c j σ ( t )+ δ λ i j ¯ σ ( t ) (cid:104) ˆ c † i ¯ σ ( t ) ˆ c j ¯ σ ( t ) − ˆ c † j ¯ σ ( t ) ˆ c i ¯ σ ( t ) (cid:105) ˆ c i σ ( t ) | n (cid:105) , (C2)where we further introduced z ( m , n ) = (cid:112) ( e − β E m + e − β E n ) / Z and exploited δ λ i j σ = δ λ ji σ . The time-local variation of O (cid:48) ( t ) , on the other hand, vanishes since it is an integratedquantity in λ . This follows directly from the definition of itsphysical sector, cf. Eq. (B2). We define γ l σ ix σ ( t ) = ∑ mn z ( m , n ) (cid:104) m | ˆ n i ¯ σ ( t ) ˆ c l σ ( t ) | n (cid:105) [ O (cid:48) ] ∗ x σ ( m , n ) ( t ) , (C3) γ l ¯ σ ix σ ( t ) = ∑ mn z ( m , n ) (cid:104) m | (cid:104) ˆ c † i ¯ σ ( t ) ˆ c l ¯ σ ( t ) (C4) − ˆ c † l ¯ σ ( t ) ˆ c i ¯ σ ( t ) (cid:105) ˆ c i σ ( t ) | n (cid:105) [ O (cid:48) ] ∗ x σ ( m , n ) ( t ) , and therewith obtain Eq. (45). Appendix D: High-order time propagation scheme
Finally, we like to set up an efficient numerical schemeto determine λ opt ( t ) . This should be based on a time-propagation algorithm where the error is of high order in thebasic time step ∆ t . Let us assume that for each time step theTaylor expansion of λ opt ( t ) is well defined. For each timeinterval and for arbitrary t ∈ [ t n , t n + ] we then have λ opt ( t ) = λ ( ) ( t ) = n p ∑ p = λ , p p ! t p + O ( ∆ t n p + ) if t ∈ [ , t [ , λ ( ) ( t ) = n p ∑ p = λ , p p ! ( t − t ) p + O ( ∆ t n p + ) if t ∈ [ t , t [ ,. . . (D1) where each λ -term must be considered as a tuple with com-ponents labelled by the super-index b , e.g., λ n , p = ([ λ n , p ] b ) ,where n refers to the n -th time interval, and where p runs from p = n p . Duringthe time propagation, the polynomial approximation must beupdated after each time step. This is done by fixing the co-efficients at each interfacing time t n such that Γ [ λ opt ]( t n ) = t (cid:54) = t n we then have Γ [ λ opt ]( t ) = O ( ∆ t n p + ) . Writ-ing J ( t ) ≡ J [ λ opt ]( t ) and ξ Γ ( t ) ≡ ξ Γ [ λ opt ]( t ) for short andapplying the product rule to J ( t ) λ opt ( t ) = ξ Γ ( t ) , the self-consistency condition (54) is readily rewritten in terms of theTaylor coefficients: λ n , p = J − ( t n ) (cid:32) p − ∑ r = (cid:18) pr (cid:19) [ ∂ p − rt J ( t )] t = t n λ n , r − [ ∂ pt ξ Γ ( t )] t = t n (cid:33) . (D2)Suppose that λ ( q ) ( t ) is known for all q < n , i.e., sup-pose that the propagation has been completed over the in-terval [ , t n [ . The next step is to update the coefficients. Atthis point we can exploit that J ( t ) and ξ Γ ( t ) scale like inte-grated quantities in λ under time-local variations which im-plies δ loc J ( t ) = δ loc ξ Γ ( t ) =
0. Hence, at t = t n , both areindependent of λ opt ( t n ) . We define˜ λ ( t ) = (cid:26) λ opt ( t ) if t < t n , . (D3)Then, ξ Γ ( t n ) = Γ [ ˜ λ ]( t n ) , (D4)and we are now able to solve Eq. (D2) for λ n , . The firstderivatives ∂ t J ( t ) (cid:12)(cid:12) t = t n and ∂ t ξ Γ ( t ) (cid:12)(cid:12) t = t n explicitly depend on λ opt ( t n ) = λ n , , which is now known, but are integrated quan-tities in the first derivative ∂ t λ opt ( t ) , i.e., they are indepen-dent of ∂ t λ opt ( t ) (cid:12)(cid:12) t = t n = λ n , . Therefore, the same idea canbe applied and in fact be repeated again and again until fi-nally λ ( n ) ( t ) is known up to the desired order. We em-phasize that the presented algorithm gives a fully converged λ opt ( t ) = λ ( n ) ( t ) + O ( ∆ t n p + ) for t ∈ [ t n , t n + [ within a singleiteration. A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev.Mod. Phys. , 863 (2011). H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner,Rev. Mod. Phys. , 779 (2014). M. Hochbruck and C. Lubich, SIAM J. Numerical Anal. , 1911(1997). R. Blankenbecler, D. Scalapino, and R. L. Sugar, Phys. Rev. D ,2278 (1981). S. R. White, D. J. Scalapino, R. L. Sugar, and N. E. Bickers, Phys.Rev. Lett. , 1523 (1989). F. F. Assaad and H. G. Evertz, in: Computational Many-ParticlePhysics, Vol. 739 of Lecture Notes in Physics, edited by H.Fehske, R. Schneider, and A. Weiße, pp. 277 (Springer, Berlin,2008). L. M¨uhlbacher and E. Rabani, Phys. Rev. Lett. , 176403(2008). P. Werner, T. Oka, and A. J. Millis, Phys. Rev. B , 035320(2009). M. Schir´o and M. Fabrizio, Phys. Rev. B , 153302 (2009). G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, Phys. Rev. Lett. , 266802 (2015). F. B. Anders and A. Schiller, Phys. Rev. Lett. , 196801 (2005). S. R. White and A. E. Feiguin, Phys. Rev. Lett. , 076401 (2004). G. Vidal, Phys. Rev. Lett. , 40502 (2004). J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, andF. Verstraete, Phys. Rev. B , 165116 (2016). M. Sandri, M. Schir´o, and M. Fabrizio, Phys. Rev. B , 075122(2012). K. Ido, T. Ohgoe, and M. Imada, Phys. Rev. B , 245106 (2015). G. Carleo, F. Becca, M. Schir´o, and M. Fabrizio, Sci. Rep. , 243(2012). G. Baym and L. P. Kadanoff, Phys. Rev. , 287 (1961). G. Baym, Phys. Rev. , 1391 (1962). J. M. Luttinger and J. C. Ward, Phys. Rev. , 1417 (1960). N. E. Bickers, D. J. Scalapino, and S. R. White, Phys. Rev. Lett. , 961 (1989). A. V. Joura, J. K. Freericks, and A. I. Lichtenstein, Phys. Rev. B , 245153 (2015). F. Hofmann, M. Eckstein, E. Arrigoni, and M. Potthoff, Phys.Rev. B , 165124 (2013). M. Potthoff, Euro. Phys. J. B , 429 (2003). M. Potthoff, In:
Strongly Correlated Systems: Theoretical Meth-ods , Ed. by A. Avella and F. Mancini, Springer Series in Solid-State Sciences, Vol. 171, p. 303 (Springer, Berlin, 2012). M. Potthoff, M. Aichhorn, and C. Dahnken, Phys. Rev. Lett. ,206402 (2003). C. Dahnken, M. Aichhorn, W. Hanke, E. Arrigoni, and M. Pot-thoff, Phys. Rev. B , 245110 (2004). M. Potthoff, Euro. Phys. J. B , 335 (2003). M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. Lett. ,056403 (2009). F. Hofmann, M. Eckstein, and M. Potthoff, J. Phys.: Conf. Ser. , 012002 (2016). F. Hofmann, M. Eckstein, and M. Potthoff, Phys. Rev. B ,235104 (2016). F. Hofmann and M. Potthoff, Euro. Phys. J. B , 178 (2016). P. Schmidt and H. Monien, cond-mat/0202046. J. K. Freericks, V. M. Turkowski, and V. Zlati´c, Phys. Rev. Lett. , 266408 (2006). N. Tsuji, P. Barmettler, H. Aoki, and P. Werner, Phys. Rev. B ,075117 (2014). A. J. Herrmann, N. Tsuji, M. Eckstein, and P. Werner, Phys. Rev.B , 245114 (2016). M. Balzer and M. Potthoff, Phys. Rev. B , 195132 (2011). P. Jurgenowski and M. Potthoff, Phys. Rev. B , 205118 (2013). C. Gramsch and M. Potthoff, Phys. Rev. B , 235135 (2015). D. S´en´echal, D. P´erez, and M. Pioro-Ladri`ere, Phys. Rev. Lett. , 522 (2000). C. Gros and R. Valenti, Phys. Rev. B , 418 (1993). C. Gramsch, K. Balzer, M. Eckstein, and M. Kollar, Phys. Rev. B , 235106 (2013). K. Balzer and M. Eckstein, Phys. Rev. B , 035148 (2014). L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen,H. Berger, S. Biermann, P. S. Cornaglia, A. Georges, andM. Wolf, Phys. Rev. Lett. , 067402 (2006). I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885(2008). L. V. Keldysh, J. Exptl. Theoret. Phys. , 1515 (1964) [Sov.Phys. JETP 20, 1018 (1965)]. J. Rammer,
Quantum Field Theory of Non-equilibrium States (Cambridge University Press, Cambridge, UK, 2007). R. van Leeuwen, N. E. Dahlen, G. Stefanucci, C.-O. Almbladh,and U. von Barth,
Introduction to the Keldysh formalism , vol.706 of
Lecture Notes in Physics (Springer, Heidelberg, Germany,2006). R. Rausch and M. Potthoff, Phys. Rev. B , 045152 (2017) M. Moeckel and S. Kehrein, Phys. Rev. Lett. , 175702 (2008). D. M. Kennes, J. C. Pommerening, J. Diekmann, C. Karrasch,and V. Meden, Phys. Rev. B95