Conservation laws in coupled cluster dynamics at finite-temperature

We extend the finite-temperature Keldysh non-equilibrium coupled cluster theory (Keldysh-CC) [{\it J. Chem. Theory Comput.} \textbf{2019}, 15, 6137-6253] to include a time-dependent orbital basis. When chosen to minimize the action, such a basis restores local and global conservation laws (Ehrenfest's theorem) for all one-particle properties, while remaining energy conserving for time-independent Hamiltonians. We present the time-dependent orbital-optimized coupled cluster doubles method (Keldysh-OCCD) in analogy with the formalism for zero-temperature dynamics, extended to finite temperatures through the time-dependent action on the Keldysh contour. To demonstrate the conservation property and understand the numerical performance of the method, we apply it to several problems of non-equilibrium finite-temperature dynamics: a 1D Hubbard model with a time-dependent Peierls phase, laser driving of molecular H$_2$, driven dynamics in warm-dense silicon, and transport in the single impurity Anderson model.


I. INTRODUCTION
The simulation of real-time electronic dynamics in molecules and materials is a challenge due to the manybody nature of the Hamiltonian combined with the need to propagate for many time steps in order to see the relevant phenomena. The problem is further complicated when finite-temperature effects are important, as is the case in many condensed-phase systems.
Here, we are interested in the dynamics of ab initio electronic Hamiltonians at finite-temperature in both condensed-phase and molecular systems. We include finite-temperature from the outset to address certain classes of problems. For instance, thermal effects are already important at the spin exchange energy scale in correlated materials with low-temperature electronic phase transitions 30,31 . Temperature also plays a role in the coupling of electrons with lattice vibrations [32][33][34] . Finally, the study of matter under extreme conditions such as a) R. Peng and A. F. White contributed equally to this work those found in planetary cores 35 also necessitates finitetemperature methods.
Various wavefunction methods traditionally used for ab initio zero-temperature quantum chemistry have been extended into the finite temperature regime. For example, the equilibrium finite temperature coupled cluster equations have been previously derived via the thermal cluster cumulant (TCC) method [36][37][38][39] and more recently through a time-dependent diagrammatic approach 40,41 . CI 42 and CC 43,44 methods have also been extended into the finite temperature regime by means of the thermofield dynamics language. In a similar manner, there are also several formalisms to generalize equilibrium finite temperature theories to non-equilibrium dynamics. Within a time-dependent diagrammatic language, non-equilibrium dynamics corresponds to moving time integration from the imaginary axis to the Keldysh (or Kadanoff-Baym) contour in the complex time plane 45,46 . In the context of coupled cluster theory, this diagrammatic extension leads to the Keldysh coupled cluster (Keldysh-CC) theory 47 .
In traditional Feynman diagrammatic approximations, a set of sufficient conditions, known as the "conserving" conditions 48,49 , are known. When diagrammatic approximations for the self-energy are constructed to satisfy these conditions, the resulting approximate singleparticle Green's function dynamics satisfies the physical conservation laws that, for example, relate the timederivative of local single-particle observables such as the density and momentum density to their corresponding currents. However, the Keldysh-CC method is not constructed as a conserving approximation, and this can lead to unphysical results when propagated for long times. It is this issue which we aim to correct.
It is known from time-dependent wavefunction theories that the use of time-dependent orbitals that obey an equation of motion derived from a time-dependent variational principle (TDVP) [50][51][52] leads to the conservation of 1-particle properties in the sense that Ehrenfest's theorem is satisfied for the 1-particle reduced density matrix (1RDM). At zero-temperature, this type of orbital dynamics has been used for time evolution with various ab initio methods, including time-dependent CI 53 , CASSCF 2 , DMET 13 , and CC 17,18,54,55 . Starting from this same idea, we can extend the Keldysh-CC method to include the variational orbital dynamics via a finitetemperature TDVP. This restores the consistency and local conservation laws for all 1-particle properties, and further maintains the global conservation law for the energy present in the original Keldysh-CC theory.
In Section II we review ground-state and finitetemperature coupled cluster theory for systems in and out of equilibrium. We discuss some of the deficiencies of the finite-temperature, non-equilibrium Keldysh-CC theory presented in Ref. 47 and we show how we can remedy some of these problems with orbital dynamics in analogy with ground state theories. This leads to the derivation of the Keldysh orbital-optimized coupled cluster doubles (Keldysh-OCCD) method. In Section IV, we apply the method to the model problem of a 2-site Hubbard model with a Peierls phase as well as an ab initio model of laser driven warm-dense silicon discussed in Ref. 47. In both of these, we demonstrate the exact conserving behaviour of the theory. We then further assess the numerical behaviour of the method in an application to laser driven molecular H 2 , as well as non-equilibrium transport in the 1D single impurity Anderson model (SIAM). The performance of Keldysh-OCCD on SIAM is compared to a numerically accurate benchmark result from finite-temperature DMRG.

II. THEORY
A. Coupled cluster dynamics at zero temperature The ground-state coupled cluster method is defined by an exponential wavefunction ansatz, where Φ 0 is a reference Slater determinant. The T operator is defined in a space of excitation operators indexed by ν, whereν † excites the reference determinant such that This space of excitations is usually truncated based on excitation level. For example, constraining T to the space of single and double excitations from the reference yields the commonly-used coupled cluster singles and doubles (CCSD) method. The coupled cluster energy and ampli-tudes are determined from a projected Schrodinger equation: Here, we have written the similarity-transformed Hamiltonian asH ≡ e −T He T .
Properties are computed from response theory which is, in practice, accomplished by solving for Lagrange multipliers, λ ν associated with the solution conditions (6). These appear in the Lagrangian, The Lagrange multipliers may be associated with elements of a de-excitation operator, such that Φ 0 | (1 + Λ) is the left eigenstate of the projected Schrodinger equation. The coupled cluster model can be extended to treat dynamics by allowing T and Λ to become functions of time. The appropriate equations of motion may be obtained from an action, by making it stationary with respect to variations in the T and Λ amplitudes. The resulting equations of motion are given by Here ν † is an excitation operator and ν is the corresponding de-excitation operator. The properties of such coupled cluster dynamics have been discussed in many works over the years [56][57][58][59][60][61] . In particular, we note that conservation laws are not generally obeyed and Ehrenfest's theorem (Equation 1) is not generally satisfied 55,59,60,62 . However, there are some special cases in which Ehrenfest's theorem will be satisfied. Energy will be conserved for a time-independent Hamiltonian, and particle number will be conserved also.

B. Coupled cluster dynamics with time-dependent orbitals
It has been shown that making the above action (Equation 10) stationary with respect to time-dependent orbitals leads to dynamics that satisfy Ehrenfest's theorem for all 1-electron operators 2,17,18,53,55 . This can be important because, as pointed out by Pedersen et al 55 , it leads to local conservation and gauge invariance in 1-electron properties.
Orbital dynamics can be included by adding a unitary, orbital-rotation operator to the wavefunction ansatz: Here, κ is an anti-Hermitian, time-dependent 1-electron operator. A Lagrangian of the form will lead to a biorthogonal approach termed OATDCC by Kvaal 54 . Alternatively, one may use a Lagrangian given by which has the advantage of being real and leads to orthonormal orbital equations. This is the approach taken by Pedersen et al and later by Sato et al 17 . We take the analogous approach in this work.
To derive the equations of motion for such a theory, it is convenient to define a 1-electron operator, Up to this point, we have assumed a fixed set of reference orbitals, but it is easier to represent the equations using a picture where the orbitals of the reference, and corresponding representation of operators, are timedependent. In this representation, the matrix elements of X are related to the time derivative of the molecular orbitals C so that, We find equations for the amplitudes, and a linear equation for the orbital rotation parameters (X), where α is a single excitation operator. The occupiedoccupied and virtual-virtual rotations are found to be arbitrary, so that we only need to consider occupied-virtual rotations 17,63 .
In practice, including coupled cluster singles amplitudes in the ansatz introduces a degree of redundancy that can lead to numerical problems 54 . For this reason, we will now specialize our discussion to the case of only doubles amplitudes. This is the TD-OCCD method of Sato et al 17 . With this simplification, the terms involving X in Equations 18 and 19 vanish leading to amplitude equations that are the same as those of of TD-CCD, but with time-dependent orbitals. We may write the orbital equation as a simple linear equation, where i, j label occupied orbitals and a, b virtual orbitals, and we have defined a Hermitian operator, The A-matrix has a simple form in terms of the symmetrized 1RDM computed from the T and Λ amplitudes, and the right-hand-side of the equation is given by, where i † a denotes the operator for the de-excitation a → i.
As mentioned previously, a consequence of including the stationary action orbital dynamics is that Ehrenfest's theorem is satisfied for all 1-electron operators. In Appendix A, we show that this is true for a very general class of wavefunction Ansätze. Kadanoff and Baym introduced the notion of conserving approximations for many-body theories based on Green's functions. These have the property that the self-energies satisfy the selfconsistent Dyson equation, or equivalently, are obtained as the derivative of approximate Luttinger-Ward functionals 64 . In Appendix B, we briefly review such approximations. A consequence of this construction is that the resulting Green's function dynamics satisfies important single-particle conservation laws, such as for the particle density and momentum density, as well as conserve energy. Although the truncated coupled cluster methods are not constructed as conserving approximations, the particle density and momentum density are 1-particle properties. Thus satisfying Ehrenfest's theorem via orbital dynamics makes the zero-temperature coupled cluster dynamics "conserving" for these properties.

C. Finite-temperature coupled cluster
Here we review the finite temperature coupled cluster theory presented in Refs. 40 and 41 which forms the basis of this work. This theory can be viewed as a realization of the TCC theory and is somewhat different from the thermal coupled cluster method described in Ref. 43. We intend to discuss the precise relationship between various formulations of finite-temperature CC in a future work.
The theory is most compactly expressed using the imaginary time Lagrangian (Equation 5 of Ref. 41): We note that this quantity is analogous to the action defined in Equation 10 via its structure as a time integral. The E and S kernels resemble the ground state energy and amplitude equations respectively and are given in Appendix A of Ref. 41. The amplitude equations are given by while the equations for the Λ amplitudes are given by The L kernel is also given in Appendix A of Ref. 41. Definingλ we may rewrite the Lagrangian in the form The boundary conditions of s andλ are clear from the integral equations (Equations 26 and 28), In practice, this means that the s(τ ) amplitudes are propagated with the differential equation, starting from τ = 0 to τ = β along the imaginary time axis. Once the s(τ ) are known in the interval [0, β], theλ amplitudes are computed according to, from τ = β back to τ = 0. Using these differential relations, it is insightful to re-express the Lagrangian as (33) where we have integrated by parts and used the boundary conditions in Equation 30, and the latter expression is analogous to the form of action in Equation 10, withλ ν playing the role of the Λ amplitudes in that equation.
As with the zero-temperature theory, finitetemperature coupled cluster is not "conserving" in the sense of Kadanoff and Baym 48,49 . The equilibrium theory nonetheless provides a framework to compute properties as derivatives of the grand potential. Details relating to this response formulation of properties are discussed in Ref. 41.

D. Coupled cluster dynamics at finite temperature
FT-CC theory can be extended to treat out-ofequilibrium systems using the Keldysh formalism as discussed in Ref. 47. The key idea is to analytically continue the imaginary time formalism of Section II C onto the real axis using a contour like the one shown in Figure 1. Computation of the response density matrices along this contour provide observables of a thermal system driven out of equilibrium as discussed in Appendix A of Ref. 47. Extending the equilibrium Lagrangian onto the Keldysh contour yields where the integral over t is understood to proceed along the contour. This Lagrangian leads to the same differen-tial equations, now expressed along the real axis, These equations are the analytic continuation of Equations 31 and 32. When ν is restricted to singles and doubles, these equations define the Keldysh-CCSD method.
Using the differential relations, we can similarly rewrite the Lagrangian as In practical calculations, there is a degree of freedom in choosing the position of the real branch of the contour relative to the imaginary-time integration limits. This corresponds to the choice of σ in Figure 2. Furthermore, once s andλ are known along the whole of the imaginary branch of the contour, they can be propagated simultaneously in real time, starting from t = −iσ. In effect, s is propagated forward along the forward branch of the contour andλ is propagated backwards along the reverse branch of the contour so that at each time the density matrices for the non-equilibrium system can be constructed.
The Keldysh-CC method has two deficiencies which appear when propagating for longer times. (i) As is the case for zero-temperature CC dynamics, Keldysh-CC truncated to singles and doubles does not in general satisfy Ehrenfest's theorem. (ii) The time-dependent properties are not stationary under propagation with the same time-independent Hamiltonian that generates the equilibrium density matrix. (However, the energy remains stationary under propagation by a time-independent Hamiltonian, for the same reasons that it is conserved The contour used in this work (red) with propagation directions for s andλ shown in blue and green respectively. Note that some branches of the contour are drawn slightly away from the actual contour so that all branches can be shown clearly.
in zero-temperature coupled cluster dynamics). (i) was discussed already in Ref. 47 which, in particular, showed analytically and numerically that approximate Keldysh-CC dynamics violates global particle number conservation. As suggested by previous arguments in Section II B, (i) can be addressed by introducing stationary orbital dynamics. (ii) is a different problem, and we will discuss the "stationarity of equilibrium approximations" in Section II F.

E. Orbital-dependent FT-CC dynamics
We now extend the formulation of Section II D to a time-dependent orbital basis with the goal of satisfying Ehrenfest's theorem for all 1-particle observables at finite temperature. As in the zero-temperature case, we reexpress the time-dependence of the orbitals via the timedependence of the matrix elements of operators in the orbital basis. For example, the representation of the 1particle part of the Hamiltonian undergoes dynamics as where here r and s refer to the starting orbitals. In analogy with the zero-temperature theory, we define an X operator (see Equation 17) for non-redundant pairs, X ai (t) = −X ia (t) * . In the finite-temperature case, recall that i, j, . . . and a, b, . . . are used to refer to orbitals with hole and particle character respectively, but these orbitals will span the full space in general. (In principle, we can introduce the matrix elements in the pure-hole and pure-particle sector, X ij (t), X ab (t), but as shown in Appendix A, by the choice of Lagrangian below, such elements vanish). As in the zero-temperature case, the introduction of orbital rotations eliminates the need for singles amplitudes. We will refer to the resulting doubles theory as Keldysh-OCCD.
To derive the equations, we define a symmetrized Lagrangian, The Lagrangian is symmetrized to ensure that the X matrix is anti-Hermitian, and the second and third terms correspond to the thermal Hartree-Fock contribution. L has a similar form to Equation 37 (or equivalently Equation 34) but we choose to remove the term containing ∆ ν (further discussion below), and similarly for h ia . As in Ref. 41, our notation includes factors of the square root of the occupation numbers in the definition of the tensors. For example, (iii) the Fock matrix becomes time-dependent and is computed from the time-dependent Hamiltonian tensors, and unlike in Keldysh-CC, the diagonal is not subtracted (further discussed below) where in the Keldysh-CC formulation, ε a , ε i are orbital energies, and n i ,n a = 1 − n a are orbital occupancies.
In the Keldysh-CC equations, the ∆ ν term in the Lagrangian and the corresponding subtraction of occupancy weighted eigenvalues from the Fock operator both originate from the choice of zeroth order Hamiltonian, and together ensure that the diagrams of Keldysh-CC correspond to a well-defined time-dependent perturbation theory. The mathematical effect of these terms is to introduce eigenvalue time-dependent phases on the indices of the amplitudes during the propagation, as seen from the amplitude equations 35, 36. In the orbital-optimized Keldysh-CC, such time-dependent phases are fully determined by stationarity of the Lagrangian with respect to the rotation elements X ij , X ab . Thus we can remove both terms involving the zeroth order eigenvalues discussed above, while obtaining the same result at stationarity. As shown explicitly in Appendix A this choice leads to X ij , X ab = 0.
To obtain the amplitude equations, we set the variations of Q with respect to the amplitudes to zero, which is equivalent to setting the amplitude variations of L to zero. Since all terms containing singles amplitudes vanish, f ia and f ai do not enter the kernels directly. Therefore, as in the analogous zero-temperature theory, the appearance of X ai in h ai (t) does not affect the Keldysh-OCCD amplitude equations. This leads to Equations (35) and (36) but with the S and L kernels containing time-dependent matrix elements, and the ∆ ν contribution missingṡ Precise equations for the kernels are given in Appendix D.
To obtain an equation for X, we vary Q with respect to the orbital parameters and set the resulting expression to zero. As in the zero-temperature case, this yields a linear equation of the form where we again use the Hermitian R defined in Equation 22. The elements of these tensors are found to be and Here we have used d pq to represent elements of the symmetrized reduced density matrices and where d ps uv is an element of the 2-particle reduced density matrix (2RDM). As in the zero temperature case, the inclusion of such orbital dynamics leads to the satisfaction of Ehrenfest's theorem for all 1-particle properties, as shown in Appendix A.

F. Stationarity of equilibrium approximations
In general, there is no guarantee that an approximate equilibrium density matrix will be stationary when propagated in real-time with the equilibrium Hamiltonian. This property holds for conserving approximations because the approximate Luttinger-Ward functional is expressed in terms of Feynman diagrams, where each diagram implicitly includes all time-orderings of the interactions. For a more general class of theories, this suggests that a necessary condition for stationarity is that all time-ordering counterparts of a particular time-ordered diagrammatic contribution to the density matrix should be included in the theory. This is clearly violated in approximate coupled cluster theory. Perturbation theory at finite order includes all time-orderings at a particular order which means that the density defined strictly as the derivative of the action with respect to an external potential has this property (see Appendix C for a more detailed demonstration). For example, consider perturbation theory at 2nd order (PT2) with a 1-particle perturbation. In Figure 3 (a) we show the diagrammatic contributions to the PT2 one-particle density matrix obtained by differentiating the energy expression with respect to the applied potential. The density matrix defined in this way includes both time orderings of the relevant diagram and has the stationary property. However, one often includes a contribution from the response of the reference via terms like those in Figure 3 (b). In this case, not all time-orderings of the same diagrammatic contribution are included and the density is no longer stationary when propagated in real time.
In the case of perturbation theory, this property can be restored by including some higher order terms. For coupled cluster theory with limited excitations (such as truncation to singles and doubles) this is not generally possible, and, as we have stated previously, none of the methods discussed in this work have this stationarity property. For short or moderate propagation times this can be corrected by subtracting the anomalous dynamics. To be precise, given a Hamiltonian of the form we can define a corrected density matrix where ρ(t) is the density propagated with the full H(t) and ρ 0 (t) is the density propagated with the timeindependent Hamiltonian H 0 .

III. IMPLEMENTATION
In the following Section IV, the Keldysh-CCSD results are obtained using the implementation as described in Ref. 47. The Keldysh-OCCD results are obtained as follows: (i) the equilibrium amplitudes s(τ ),λ(τ ) are first computed as described in Section II C. In particular, the differential form of the s-amplitude equation (Equation 31) is first propagated from τ = 0 to τ = β along the imaginary time contour. The s(τ ) amplitudes are then used in Equation 32 to propagateλ(t) from τ = β back to τ = 0. These equilibrium amplitudes are computed within the coupled cluster doubles approximation, using fixed orbitals, and without any truncation of the occupied or virtual space. The 4th order Runge-Kutta scheme is used to propagate the differential equations. (ii) The equilibrium amplitudes at τ = β/2 are taken as the initial t = 0 amplitudes for the subsequent dynamics. This corresponds to a choice of contour where the real part extends from −iβ/2 in Figure 2. (iii) The dynamical s(t) andλ(t) amplitudes are computed as described in Section II D. In particular, Equation 35, 36 are simultaneously propagated from t = 0 to t = t f on the real contour where the kernels S, L are modified to include orbital dynamics as described in Section II E. (iv) At each update of s(t) andλ(t), the orbital Equation 47 is solved to obtain R, which is used to update the integrals according to e.g. equation 38. We use the 4th order Runge-Kutta (RK4) scheme to update both the amplitudes and the orbitals in a coupled manner as described below. The amplitude Equations 35, 36 and orbital coefficient Equation 17 are of the form dy dt = f (t, y(t), h(C(t))) dC dt = CX(t, d(y(t)), h(C(t))) where y ∈ {s(t),λ(t)}, h(C(t)) denotes the Hamiltonian integral matrix elements in the time-dependent orbital basis C(t), and d(y(t)) denotes the reduced density matrices computed from the amplitudes at time t. In RK4, the time-step from t → t + δt is assembled from intermediate amplitudes y i and orbitals C i as well as their respective finite difference δy i and where t i = t + a i δt y i = y(t) + a i δtδy i−1 where a i are standard 4th order Runge-Kutta weights {0, 1 2 , 1 2 , 1}. Note that all quantities for i = 1 correspond to their values at time t.
After all intermediate quantities are obtained, the amplitudes and orbitals are updated as y(t + δt) = y(t) + 1 6 δt(δy 1 + 2δy 2 + 2δy 3 + δy 4 ) X(t + δt) = 1 6 (X 1 + 2X 2 + 2X 3 + X 4 ) The present Keldysh-OCCD theory has the same scaling as the previous Keldysh-CCSD theory in terms of computational cost and memory. The amplitude update cost scales as N 6 for each real or imaginary time step, where N is the size of 1-particle basis. Furthermore for each real time step, Keldysh-OCCD additionally involves solving a linear set of equations for N 2 variables in the orbital equation 47, and an integral update which scales as N 5 .

A. Numerical demonstration of Ehrenfest's theorem
We first demonstrate the satisfaction of Ehrenfest's theorem for the Keldysh-OCCD method for the 2-site time-dependent Peierls-Hubbard model. This was previously studied with Keldysh-CCSD in Ref. 47. The timedependent Hamiltonian is given by where the first term is the Peierls driving term, which mimics the effect of the coupling of an underlying nuclear lattice to an external laser pulse. Here, the driving takes the form We use U = 1.0, for which a chemical potential of µ = 0.5 gives half filling at equilibrium, a temperature T = 1.0, and pulse parameters of σ = 0.8, t 0 = 2, ω = 6.8. We show the Keldysh-CCSD and Keldysh-OCCD results for the population difference between the 2 sites (L, R) defined as n L − n R , where n L = n L,↑ + n L,↓ , along with the exact result computed by propagating the density matrix in the full Liouville space. As shown in Figure 4, the deviation from the exact result increases for both Keldysh-CCSD and Keldysh-OCCD with increasing pulse amplitude A 0 . However, the orbital dynamics in Keldysh-OCCD ensures that it is much closer to the exact result for all pulse parameters. As discussed in detail in Ref. 47, Keldysh-CCSD does not in general conserve global symmetries, such as the total particle number, whereas such 1-particle quantities are conserved both locally and globally by Keldysh-OCCD. In Figure 5 the 2-site Hubbard model for µ = U/2 where the particlehole symmetry is no longer present. As expected, we see that as µ is decreased from half-filling, the Keldysh-CCSD total particle number begins to deviate from the equilibrium value at longer times, whereas the total particle number remains conserved for Keldysh-OCCD for each µ.
In Figure 6, we demonstrate the stronger condition of conservation of local particle number for the population of the left site. From Ehrenfest's theorem, where the expression on the right hand side is the contribution from the right-site flux. The red curve in the upper panel shows the local population change computed with the flux. The green curve gives the timederivative of n L , computed by the finite-difference expression ( n L (t + δt) − n L (t))/δt, with time step 5×10 −3 . The two curves are on top of each other, thus we also plot the difference between the r.h.s. and the l.h.s.
quantities. This stays close to 0 at all times, demonstrating the local conservation law. In the lower panel, we also plot the difference for a smaller time-step of 2.5 × 10 −3 where the quantity is reduced by a factor of 2. This illustrates that Ehrenfest's theorem will be fully satisfied for an infinitesimal time step.

B. Comparison with Keldysh-CCSD: warm-dense Si
We next compare the performance of Keldysh-CCSD and Keldysh-OCCD for the field-driven Si system in Ref. 47. This treats a single primitive cell of Si in a minimal basis (SZV 66 with GTH-Pade pseudopotentials 67,68 ) and with the ions frozen at the experimental lattice constant (3.567Å). The matrix elements were obtained from the PySCF software package using plane-wave density fitting [69][70][71] . We use the dipole approximation in the velocity gauge where the coupling to an external field is of the form where p is the momentum operator, and we choose a pulse shape described by We choose parameters of σ = 2.0, t 0 = 15.0, ω = 0.97529, A 0 = 0.6752 at a temperature T = 0.2. This corresponds to a maximum laser intensity of 2.36 × 10 16 W/cm 2 . Figure 7 again demonstrates that, as expected, Keldysh-OCCD conserves total particle number in contrast to Keldysh-CCSD. In Figure 8  at short times where we would expect both approximations to be valid. However at longer times, different physics is predicted: Keldysh-CCSD yields a sharp drop in the valence to conduction population transfer at the point where there is a strong violation of total particle number conservation, while Keldysh-OCCD predicts that the valence to conduction population transfer continues at a reduced rate after the pulse ends.

C. Molecular H2 in a laser field
We further apply the Keldysh-OCCD theory to a fielddriven molecular system. The purpose of this calculation is to provide a simple benchmark so that future implementations of Keldysh-OCCD and similar theories can be easily tested. The total molecular Hamiltonian is de- H(0) is the time-independent molecular Hamiltonian, and is the molecular dipole operator for N electrons and N A nuclei. We use a field of the same form as in Equation 61.
In Figure 9, we show the x-component of the dipole moment for H 2 in the STO-3G basis 72 with the molecule aligned along the x-axis. We use field parameters t 0 = 15.0, σ = 2.0, ω = 1.0, z = (1.0, 1.0, 1.0) and T = 1.0, µ = 0. The raw data for µ x (t) are provided in the supporting information. . µx from Keldysh-OCCD as a function of time for a field-driven molecular H2 benchmark. The shape of the electric field as a function of time is also plotted in the lower panel.

D. Single Impurity Anderson Model
Finally we apply the Keldysh-OCCD method to transport in the single impurity Anderson model, with a central impurity ("dot") coupled to two one-dimensional leads. We use a Hamiltonian of the form whereĤ dot = V g n d + U n d↑ n d↓ (65) as used previously in Ref. 13. Here the impurity is associated with fermion operators a ( †) d and its Hamiltonian is parametrized by a gate voltage V g and Hubbard interaction U , while the left and right leads are associated with fermion operators a Rp , respectively, and are described by tight-binding Hamiltonians. In the following calculations, the equilibrium state is generated by the Hamiltonian with parameters t leads = 1.0, t hyb = 0.4, and zero bias V = 0 for various specified temperatures and values of U . For all temperatures, the total system has the same number of particles as the number of sites with total S z = 0. A Hartree Fock calculation is performed at zero temperature, and these orbitals are used in the equilibrium CC calculations.
The dynamics is then generated by applying a small bias, V = −0.005, and the other parameters are kept fixed. The dynamics can be characterized by the timedependent current across the dot, which we compute as the average of the current between the dot and its closest left and right neighbors J(t) = (J L (t) + J R (t))/2, where and the bracket denotes the expectation value with respect to the Keldysh-OCCD density matrix. Figure 10 shows an example of the time-dependent current divided by the bias for several different temperatures at U = 1.0. As discussed, for example, in Ref. 73, the current will quickly reach its steady state in the infinitesize limit, but for finite leads, the finite system size produces an oscillatory behavior. In the case of 16 sites, we propagate for up to a time of t = 10.0 corresponding to half the oscillation period, which is sufficient to extract the physics of the system pertaining to the infinite-size steady state. Figure 11 shows the conductance G as a function of gate voltage V g for the 16 site model and U = 1.0 at various temperatures as computed from Keldysh-OCCD, as well as from reference density matrix renormalization group (DMRG) results. The conductance is computed as the current divided by bias averaged over the plateau region from t = 2.0 to t = 8.0. The zero-temperature Keldysh-OCCD result is computed using a zero temperature implementation similar to that described in Sato et al 17 . The reference DMRG results at both zero-and finite-temperature are computed using a time-step targetting time-dependent DMRG method 74-77 as implemented in the PyBlock3 software package 78,79 interfaced with the HPTT library 80 . The "Kondo peak" in the lowtemperature limit can be observed in the shape of a high conductance plateau, arising from many-body effects. As described elsewhere (see e.g. Ref. 81), the Kondo resonance marks the increase of the density of states of the impurity around the Fermi surface of the leads due to the spin interaction between a particle near the Fermisurface of the leads, and a particle on the impurity, and this resonance results in a high tunneling probability. The Kondo effect is only observed at temperatures below T K ∼ e −1/jρ where j is the (typically small) effective exchange coupling between a particle in a lead state and a particle on the impurity, and ρ is the lead density of states at the Fermi-energy. We find that the plateau at T = 0 does not reach the predicted G = 1/π unitary value at −U/2. The deviation from the unitary limit has been seen in previous calculations and can be attributed to finite-size errors 73 (for example, Figure 6(b) of Ref. 13 plots the conductance at V g = −U/2 as a function of system size, demonstrating the convergence to the unitary value G = 1/π (or more precisely, for the unit used in that work, 2e 2 /h) as the system size approaches infinity). is very sensitive to the window over which the current is averaged. This can be seen from computing the conductance in two different ways, i.e. as an average over the current divided by bias from t = 2.0 to t = 8.0 as shown in red, and as the current divided by the bias value at t = 2.0 as shown in green. As shown in the top panel, different averaging windows result in both an overall vertical shift of the conductance, and a small change in the shape of the curve. (ii) For the non-interacting case in the lower panel, where the Keldysh-OCCD method is exact, there is still a difference between the Keldysh-OCCD conductance and the DMRG conductance, coming from the bond dimension truncation. Thus given the sensitivity of G to the averaging window of the current, as well as the DMRG bond dimension truncation error, the Keldysh-OCCD and DMRG results in Figs. 11 and 12 are in very good agreement.
The temperature-dependence of the conductance at V g = −U/2 (U = 1.0) is shown in Figure 13, where the inset plots the Keldysh-OCCD and DMRG results with data points marked and the temperature axis on a log-scale. Analytic treatments indicate that the peak conductance has a logarithmic dependence on temperature 81 which is consistent with our data shown in the inset. Thus our results show that at moderate interaction strength, Keldysh-OCCD is able to correctly capture the physics of the single impurity Anderson model, including the non-equilibrium Kondo physics and its temperaturedependence.

V. CONCLUSIONS
In this work, we present a modification of the Keldysh coupled cluster theory that restores local and global conservation of one-particle quantities via an optimal orbital dynamics. On a variety of models and simple ab initio systems, we have demonstrated that such conservation laws are indeed obeyed within the Keldysh orbital-optimized coupled cluster doubles approximation (Keldysh-OCCD), with a concomitant improvement of the predicted dynamics, especially at longer times. In the single impurity Anderson model, we qualitatively reproduce the temperature dependent transport physics, including that associated with the Kondo plateau. We believe this will be useful in the ab initio modeling of Kondo transport in the future.
However, there remain important challenges in the practical application of the Keldysh coupled cluster formalism: (i) at the doubles level, the cost and memory scaling remains a barrier to many interesting applications. To decrease the computational cost, a potential future direction is to replace the optimal orbital dynam- ics by the orbital dynamics of time-dependent Hartree-Fock. This will not exactly preserve Ehrenfest's theorem for 1-particle dynamics as does the present Keldysh-OCCD approximation. However, for weakly interacting systems, the Hartree-Fock orbital update should still give better results than Keldysh coupled cluster approximations without orbital dynamics. (ii) At a formal level, the fact that the state generated by the finite-temperature coupled cluster theory is not in general a stationary state of the dynamical theory leads to an ambiguity in the definition of the equilibrium state. Here we show that Ehrenfest's theorem is restored for one-particle properties by including the optimal orbital dynamics into a zero-temperature time-dependent wavefunction ansatz. Let Ψ(y ν , p) be the time-dependent wavefunction ansatz, where y ν (t) are the variational parameters (e.g. such as the CI coefficients) and p(t) is the orbital basis. The equation of motion of y ν (t) and p(t) can be determined from a time-dependent variational principle 2,13,17,53,63 by making the action, stationary w.r.t small variations where the variation w.r.t. the orbitals, ∆|Ψ , is parameterized by an anti-Hermitian 1-body operator, ∆ 13 . Setting δS = 0 leads to where solving equation (A4) for each orbital pair is equivalent to enforcing Ehrenfest's theorem for the 1-particle density matrix elements, and hence for any 1-particle property. Furthermore, the energy is conserved: the timedependence of Ψ can be expressed as where X is anti-Hermitian and parameterizes the timedependence of the orbitals 13 . Then the energy derivative where we use equation (A3) and (A4) for the 3rd equality, and equation (A5) for the 4th equality. At finite temperature, the idea is analogous: orbital optimization based on an action principle can be used to satisfy Ehrenfest's theorem.
We first derive the conditions satisfied by the stationary orbital dynamics. The Lagrangian in equations 39 and 40 can be written as where and where the term −iX comes from the modification of the integrals as in equation 41. We use the notation ... 0 to denote an expectation value with respect to the mean-field 1-RDM p pq in equation D6 and 2-RDM p pq rs = p pr p qs − p ps p qr , and ... CC N to denote an expectation value with respect to the coupled-cluster normalordered 1-and 2-RDMs defined in equations D7-D15. As in the zero-temperature case, the orbital variation can be parameterized by an antihermitian matrix ∆ and then the variation of the modified Hamiltonian is which can be compactly denoted as [H − iX, ∆]. In the orbital variation of the action Q, the variation of terms A8, A11 comes from the variation of the modified Hamiltonian tensors and the variation of term A12 gives whereḋ is the time-derivative of the CC 1-RDM in Equation D4. Thus the orbital gradient of Q is where ... CC denotes the expectation value computed using the CC 1-and 2-RDMs given in equations D4, D5. Setting ∆Q = 0 for each orbital pair ∆ uv giveṡ where the F matrix is defined in equation 50. For Keldysh OCCD, only the "occupied-occupied" and "virtual-virtual" blocks of the 1-RDM are non-zero. This means that the only non-zero blocks of Equation A19 arė For an anti-hermitian matrix X, the final two equations are complex conjugates of each other and we only need to satisfy one of them, which we rewrite as where we have defined R ≡ −iX. This is the precise form of the orbital equation previously shown in Equation 47 that comes from stationarity of the Lagrangian. Furthermore, it can be shown that equations A20, A21 trivially hold for any occupied-occupied and virtual-virtual rotation, and hence those rotations can be set to zero. We can write Ehrenfest's theorem in a rotating basis at finite temperature in terms of the reduced density matrices and matrix elements of the operator A in the orbitals defined in Equation 17, Since this equation should hold for any operator, it must be separately satisfied for each uv pair. This immediately yields equation A19, where the reduced density matrices correspond to their definition in coupled cluster theory. Thus, we see that the stationarity of the coupled cluster Lagrangian with respect to orbital variations leads to the satisfaction of Ehrenfest's theorem.

Appendix B: Conserving approximations
There are several equivalent ways to define a conserving approximation in the theory of Green's functions. To draw a close correspondence with the approach in this work, we define a conserving approximation to be one where the Green's function on the Keldysh contour G pq (t 1 , t 2 ) = iT C a p (t 1 )a † q (t 2 ) (where T c indicates contour time-ordering) makes the Luttinger-Ward functional A stationary, defined as The density matrices which appear in the orbital equation (Equation 50) are also used to compute properties. They can be obtained from the derivative of the Lagrangian with respect to the potential: Here and throughout this work we have assumed that the matrix elements include factors of the square root of the occupation numbers as in Equation 42.