Time-dependent coupled cluster theory on the Keldysh contour for non-equilibrium systems
Alec F. White Garnet Kin-Lic Chan
cond-mat.stat-mech physics.chem-ph




whiteaf@berkeley.edu

gkc1000@gmail.com
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena,
California 91125, USA


We leverage the Keldysh formalism to extend our implementation of finite temperature coupled cluster theory [J. Chem. Theory Comput. 2018, 14, 5690-5700] to thermal systems that have been driven out of equilibrium. The resulting Keldysh coupled cluster theory is discussed in detail. We describe the implementation of the equations necessary to perform Keldysh coupled cluster singles and doubles calculations of finite temperature dynamics, and we apply the method to some simple systems including a Hubbard model with a Peierls phase and an ab initio model of warm-dense silicon subject to an ultrafact XUV pulse.

Time-dependent coupled cluster theory on the Keldysh contour for non-equilibrium systems
    Garnet Kin-Lic Chan
    July 2019
========================================================================================





§ INTRODUCTION

Understanding the behavior of condensed phase systems driven out of equilibrium is key to the development of technologies
that can make use of exotic non-equilibrium electronic properties  in real devices. Ultrafast spectroscopy provides the experimental means to probe non-equilibrium dynamics, and  has been used to study transient electronic phases of materials.<cit.>
However, the direct simulation of the dynamics is complicated by several factors including:

    
  * The quantum many-body problem for a realistic material
    
  * The statistical mixture of many-body states implied by a finite temperature
    
  * The generically non-thermal, time-dependent distribution of states generated by the perturbation
    
  * The coupling between electronic and nuclear degrees of freedom.


For simplicity in this study we will neglect the description of electron-nuclear coupling.
However, points 1.-3. already pose major challenges for an electronic treatment.
In ab initio real time dynamics, 
there have been applications of real-time density functional theory (DFT) to both
materials as well as molecules<cit.>,
while time-dependent ab initio wavefunction methods such as time-dependent configuration interaction <cit.> and coupled cluster methods have been introduced in molecules<cit.>.
These works however, have typically neglected temperature either due to the complications introduced, or because it is not of interest
to the phenomenon studied. Finite-temperature electron dynamics has been studied using model Hamiltonians with a variety of methods,
such as time-dependent exact diagonalization<cit.>, time-dependent density matrix renormalization group methods<cit.>, diagrammatic Monte Carlo<cit.>, and hierarchical
equations-of-motion<cit.>. However, the ab initio extensions of such
studies remain to be developed.
















Here, we present a method to model finite-temperature electronic dynamics that
is suited to ab initio simulations of condensed phase systems as well as molecules.
We work within the coupled cluster framework, which has enjoyed great popularity for solving the zero-temperature electronic structure problem,<cit.> and where extensions to periodic solids<cit.>, finite temperatures<cit.>, and open systems <cit.> have become areas of recent research.
From a theoretical perspective, one can either start from the equation-of-motion for the density matrix or from a thermofield formalism  which propagates a pure state in an extended Hilbert space<cit.>.


We work within the density matrix language, and use the Keldysh formalism<cit.>
to provide a straightforward generalization of the imaginary-time finite-temperature CC equations in Ref. White2018
to describe the dynamics of finite temperature systems.

We show some initial applications to model systems and a simple ab initio representation of high-temperature silicon under the influence of an ultrafast extreme ultraviolet (XUV) pulse.








§ THEORY

We will describe the theory in 4 steps:

    
  * In Section <ref> we describe the working equations of the FT-CC method presented in Ref. White2018
    
  * In Section <ref> we describe some important aspects of the Keldysh formalism
    
  * In Section <ref> we generalize the imaginary time FT-CC to out of equilibrium systems using the Keldysh formalism
    
  * In the final Sections (<ref>-<ref>) we discuss several important properties and extensions of the theory.




 §.§ Finite temperature coupled cluster

In Ref. White2018 we presented a finite temperature coupled cluster theory which provides an ansatz for the correlation contribution to the grand potential Ω such that

    Ω = Ω^(0) + Ω^(1) + Ω_CC

where Ω^(0) + Ω^(1) is the grand potential of the thermal Hartree-Fock theory and Ω_CC is the approximation
computed by coupled cluster.
Our strategy was to generalize the diagrammatic rules of finite temperature perturbation theory to the diagrammatic representation of the coupled cluster iteration. Interestingly, this statregy ultimately led to equations identical to those of Mukherjee's thermal cluster cumulant theory (TCC)<cit.> (the relationship between the thermal Wick expansion used in TCC and thermofield dynamics is also
discussed in Mukherjee's work <cit.>).
The coupled cluster correction to the grand potential is computed from time-dependent amplitudes, s_μ(τ), evaluated in imaginary time:

    Ω_CC   =   1/4β∑_ijab⟨ ij||ab⟩∫_0^βdτ
    	[s_ij^ab(τ) + 2s_i^a(τ)s_j^b(τ)] 
       +   1/β∑_iaf_ia∫_0^βdτ s_i^a(τ).

Here, β is the inverse temperature,  we use f and ⟨ pq || rs⟩ as the 1-particle and antisymmetrized 2-particle integrals of the interaction, and the sums run over all orbitals. The set of generalized excitation operators enumerated by μ defines the truncation of the theory, and the amplitudes are determined from a set of non-linear Volterra integral equations of the form

    s_μ(τ) = - ∫_0^τdτ' e^Δ_μ(τ' - τ)S_μ[s(τ')].

Δ_μ denotes a difference of orbital energies, and the kernel, S_μ, is a nonlinear function of the amplitudes that is local in imaginary time. It contains contractions that are identical in graphical representation to those of ground state coupled cluster theory, the only difference algebraically being the inclusion of Fermi-Dirac occupation numbers and sums that run over all orbitals. Explicit expressions for the FT-CCSD S_1 and S_2 kernels are given in the appendix of Ref. White2018.

Properties of the thermal system can then be computed by differentiation of Ω. In practice, this is accomplished by defining a Lagrangian, 


    ℒ≡1/β∫_0^βE(τ) -
    	1/β∑_μ∫_0^β dτλ_μ(τ)[
    	s_μ(τ) + ∫_0^τdτ' e^Δ_μ(τ' - τ)S_μ(τ')].


(Equation 37 of Ref. White2018) with Lagrange multipliers λ_μ(τ) such that ℒ coincides with the free energy at its stationary point but is a variational function of λ_μ(τ) and s_μ(τ). Variational optimization of ℒ with respect to λ_μ yields the amplitude equations (Equation <ref>) and variational optimization with respect to s_μ yields linear equations for λ_μ of the form 

    λ_μ(τ) = - L_μ[s(τ),λ̃(τ)].

The L_μ kernel is also local in imaginary time and contains the same contractions as the ground state λ equations. We use the notation given by

    λ̃_μ(τ) = ∫_τ^βdτ'e^Δ_μ(τ - τ')λ_μ(τ').

The discretized versions of these equations are given in Appendices A and C of Ref. White2018.



 §.§ The Keldysh formalism for non-equilibrium systems

Our approach to the non-equilibrium problem is based on the Keldysh formalism.<cit.> This approach is conceptually appealing because the 3 generic difficulties mentioned in Section <ref> are treated on the same footing. Many-body correlation, initial finite-temperature, and the dynamics of a driven system are all treated naturally in a single time-dependent framework. We introduce the basic ideas and notations in Appendix <ref>. 

The key to the Keldysh formalism is the compact representation of the out-of-equilibrium average of some property as a trace of a contour-ordered exponential (see Equation <ref>). In recent years, many-body methods using this formalism have gone by the name “Kadanoff-Baym dynamics"<cit.>. Such methods commonly compute the Green's functions along the Keldysh contour,<cit.> but we will instead use the non-equilibrium analog of the free energy to compute properties. We will refer to this quantity as the Keldysh cumulant generating function.

Consider the grand potential represented as the logarithm of the trace of a time-ordered exponential:

    Ω = -1/βlnTr[e^-β K_0𝒯exp(-i∫_0^-iβ dt V_I(t)) ].

K_0 is the zeroth order Hamiltonian including the chemical potential (Equation <ref>), and V_I(t) is the interaction in the interaction picture.
We can rewrite this expression as an integral along the Keldysh contour C because the contribution from the forward (0→ t_f) and backward (t_f → 0) branches will exactly cancel (see Fig. <ref>):

    Ω = -1/βlnTr[e^-β K_0𝒯_C
    	exp(-i∫_C dt V_I(t)) ].


We may now introduce a functional dependence on a parameter α such that

    Ω[α] = -1/βlnTr[e^-β K_0𝒯_C
    	exp(α O(t_f)-i∫_C dt V_I(t)) ].

Notice that for α = 0 this is just the equilibrium grand potential,

    Ω[0] = Ω,

but the derivative of Ω[α] evaluated α = 0 will reproduce Equation <ref> with an extra factor of β:

    .∂Ω/∂α|_α = 0   = -1/β ZTr{e^-β K_0
       ×𝒯_C[exp(-i∫_C dt V_I(t))O(t_f)]
    	}.

This means that the derivatives of Ω[α], the Keldysh cumulant generating function, are the observables of the non-equilibrium system

    O(t_f) = -β. ∂Ω[α]/∂α|_α = 0




 §.§ Non-equilibrium properties from Keldysh coupled cluster

In order to compute this quantity as a straightforward extension of finite temperature coupled cluster, we slightly rewrite Ω[α] as


    Ω[α] = -1/βlnTr[e^-β K_0𝒯_C
    	exp(-i∫_C dt [V_I(t) + iα O(t)δ_C(t - t_f)]) ].


The expression is now analogous to what we compute in the equilibrium version of finite temperature coupled cluster. The only difference is that we have a modified potential,

    V_I(t) + iα O(t)δ_C(t - t_f),

where δ_C indicates that the delta function acts only once along the contour. In practice there is great flexibility in how the contour is chosen as long as it passes through t_f at least once.

FT-CC provides an approximation to Ω by imaginary time-integration of non-linear amplitude equations. We now extend this same approximation to systems out of equilibrium by solving the same non-linear amplitude equations along the Keldysh contour. This leads to complex amplitudes s_μ(t) where the time argument, t, runs along the Keldysh contour. The coupled cluster contribution to the grand potential is given by

    Ω_CC   =   i/4β∑_ijab⟨ ij||ab⟩∫_Cdt
    	[s_ij^ab(t) + 2s_i^a(t)s_j^b(t)] 
       +   i/β∑_iaf_ia∫_C dt s_i^a(t).

For computational reasons, the integrals in the amplitude equations are broken up so that


    s_μ^f(t)    = -i∫_0^tdt' e^iΔ_μ(t' - t)S_μ[s^f(t')] 
    
        s_μ^b(t)    = e^iΔ_μ(t_f - t)s_μ^f(t_f)-i∫_t_f^t dt' e^iΔ_μ(t' - t)S_μ[s^b(t')] 
    
        s_μ^i(τ)    = s_μ^b(0)e^-Δ_μτ -
        ∫_0^τdτ'e^Δ_μ(τ' - τ)S_μ[s^i(τ')]


where we have used the superscripts f,b,i to refer to the forward, backward, and imaginary parts of the Keldysh contour (see Figure <ref>). Explicit expressions for the kernel, S_μ, of Keldysh-CCSD are given in Appendix <ref>.

We now define a Lagrangian in analogy with Equation <ref>:


    ℒ≡i/β∫_CdtE(t) -
    	i/β∫_C dt λ^μ(t)[
    	s_μ(t) + ∫_C(0)^C(t)dt' e^iΔ_μ(t' - t)S_μ(t')].


We have used C(t) to indicate the limits of the integral along the Keldysh contour. In Appendix <ref>, we discuss the equations determining the multipliers, λ^μ(t), along the Keldysh contour. From the amplitudes and Lagrange multipliers, we can construct (unrelaxed) coupled cluster density matrices at each time along the contour as described in Appendix <ref>. This means that once we solve the CC equations on the Keldysh contour, we can compute any 1-particle property for any real time t≤ t_f. Two particle properties can also be computed, though the additional computational cost is significant. Some details of the implementation are discussed in Section <ref>.



 §.§ Correlation functions from higher-order response

The theory can also be extended to the computation of correlation functions,

    O_1(t_1)O_2(t_2),…

in or out of time order. We must consider a more general Keldysh cumulant generating function


    Ω[α_1,α_2,…] = -1/βlnTr[e^-β K_0𝒯_C
    	exp(-i∫_C dt [V_I(t) + iα_1 O_1(t)δ_C(t - t_1) + iα_2 O_2(t)δ_C(t - t_2) + …]) ]

 
and use the flexibility of the Keldysh formalism to choose the contour and the locations of the contour delta functions (δ_C) such that contour ordering reproduces the ordering of the operators in the desired correlation function. One can then compute the desired n-point correlation function by the nth order derivative:

    O_1(t_1)O_2(t_2),…,O_n(t_n) = β^n .
        ∂^n Ω/∂α_1∂α_2…|_α_i = 0.

Therefore  higher order responses of a coupled cluster approximation to Ω[α_1, α_2, …], constructed along the appropriate contour,
can be used to obtain the coupled cluster approximation to arbitrarily general correlation functions.



 §.§ Differential form of the Keldysh CC equations

Thus far we have expressed everything as non-linear Volterra type integral equations, but there is an equally valid differential form (see for example the TCC representation of finite temperature coupled cluster given in references Sanyal1992,Sanyal1993,Mandal2001,Mandal2003). We start with the integral equations for s and λ̃ along the Keldysh contour:

    s_μ(t)    = -i∫_C(0)^C(t)dt'e^iΔ_μ(t' - t)S_μ[s(t')] 
    λ̃_μ(t)    = -i∫_C(t')^C(-iβ)dt'e^iΔ_μ(t - t')L_μ[s(t'),λ̃(t')].

Note that since the density (Equations <ref>-<ref>) depends only on λ̃, we never need to work with λ directly using this formulation. 

Taking the time derivative of Equations <ref> and <ref> yields equations of motion for s and λ̃:

    ds_μ/dt   = (-i){Δ_μs_μ(t) + S_μ[s(t)]}
    dλ̃_μ/dt   = i{Δ_μλ̃_μ(t) + 
        L_μ[s(t),λ̃(t)] }.

In both these equations, the term proportional to Δ_μ comes from the fact that we have formulated our equations in the Schrodinger picture.

These equations can be used to directly propagate the s-amplitudes and λ̃-amplitudes along the Keldysh contour. In this way it is possible to avoid storing the amplitudes as is necessary when working with the integral formulation of the equations. However, working with this differential form presents other difficulties. While the value of s_μ is initially known to be zero at t = 0, the value of λ̃_μ is known at t = -iβ. Considering the dependence of L_μ on s, this requires the calculation to be carried out in two steps. First s_μ is propagated from t = 0 to t = -iβ. This propagation can be carried out along the Keldysh contour or just along the imaginary axis. Second, both s and λ̃ are propagated back from t = -iβ to t = 0 along the Keldysh contour and the density matrix is computed along the way as per Appendix <ref>. Though this differential formulation is ultimately a more natural form for the dynamical problem, we have not yet fully investigated the numerical properties of such a method and instead use the integral formulation to investigate the basic properties and some initial applications.



 §.§ Violation of conservation laws in Keldysh-CC

When considering dynamics, and finite-temperature dynamics in particular, it is more difficult to construct approximations that yield dynamics that match even qualitatively those of the exact system. Kadanoff and Baym<cit.> stated conditions that approximate Green's functions must satisfy in order to yield dynamics consistent with all local conservation laws (such as the continuity equation relating density and current), and approximations that meet these criteria are called “conserving”. Coupled cluster theory is  not constructed as a conserving approximation. However, coupled cluster methods have been used to obtain accurate zero-temperature dynamics in small finite systems nonetheless<cit.>. Furthermore, even conserving approximations can yield qualitatively incorrect results<cit.>. We do not consider this property to be an absolute requirement for useful approximations to dynamics.

Unfortunately, the Keldysh CC presented in this work, truncated at some finite excitation level, can also violate global conservation laws. Perturbation theory approximations to the cumulant generating functional do conserve global symmetries in general as we will now show. 

Consider a cumulant generating functional, Ω, that admits a perturbation theory expansion in some interaction controlled by a coupling, λ, such that

    Ω[λ] = Ω[0] + λΩ'[0] + λ^2/2Ω”[0] + …

Note that nth order correction can be computed by finite difference in that

    Ω^[n] = .1/n!∂^n Ω[λ]/∂λ^n|_λ = 0 = ∑_xc_xΩ[λ_x]

where the λ_x are points slightly displaced from λ = 0, and the {c_x}, {λ_x} are defined by a particular finite difference approximation. If we now compute some time-dependent observable, O, by differentiating with respect to a time-dependent coupling, α(t), the nth order contribution to the observable is

    O^[n](t)    = ∂/∂α(t).Ω^[n][α(t)]|_α(t) = 0
       =∂/∂α(t).∑_x c_xΩ[λ_x,α(t)]|_α(t) = 0.

This implies that if O is globally conserved by the exact dynamics, it should also be conserved in perturbation theory. In Appendix <ref>, we present a detailed analysis of particle number conservation in an exactly solvable problem which suggests that in general, truncated Keldysh-CC does not have this property. This can limit the applicability of the method as we will discuss in Section <ref>.



§ IMPLEMENTATION

To enable a concrete implementation, we truncate the coupled cluster equations to the singles and doubles level (CCSD), where μ labels
at most four fermionic lines.
The implementation of Keldysh-CCSD closely mirrors the implementation of FT-CCSD, so we refer to Ref. White2018 for most of the details. There are two primary differences. First, everything must use complex arithmetic as a result of the real-time propagation. Second, the numerical time integration must be appropriately modified.



 §.§ Numerical integration

Conceptually, we integrate the amplitude equations (Eq. <ref>)
just as for FT-CCSD, except that we now consider each branch of the contour separately:

    s_μ^f(t_y)    = -i∑_x G_x^ye^iΔ_μ(t_x - t_y)S_μ[
        s^f(t_x)] 
    
        s_μ^b(t_y)    = s_μ^f(t_f)e^iΔ_μ(t_f - t_y)
       + i∑_xG_x^ye^iΔ_μ(t_x - t_y)S_μ[s^b(t_x)]
    
        s_μ^i(τ_y)    = s_μ^b(0)e^-Δ_μτ_y
        - ∑_xG_x^ye^Δ_μ(t_x - t_y)S_μ[s^i(τ_x)].

Again, we have used labels f, b, and i to refer to the different branches of the Keldysh contour (forward in real time, backward in real time, and imaginary time), and as in our FT-CCSD description, we have used G as a generic tensor of quadrature weights such that

    ∫_t_x^t_yf(t) dt ≈∑_xG_x^yf(t_x).


In the non-equilibrium case we compute by differentiation the response density matrix at a particular time. This means that, unlike in the equilibrium case, we must use a quadrature rule where the weights have a smooth continuum limit,

    G_x^y → G(t_x),    g_y → g(t_y)  as   n_g →∞,

where G(t) and g(t) are smooth functions.
If this condition is not satisfied, the response density matrix as a function of time will not have a well-defined limit as the grid becomes finer. This means that we cannot use quadratures like Simpson's rule or Boole's rule. In this work we use constant quadrature weights, i.e.

    g_y = 1/δ,      G_x^y =  1/δ,     x ≤ y
    
        0,    otherwise

where δ is the grid spacing. This has the advantage of simplicity, but it requires a large number of grid points to obtain high accuracy and this limits the current study to relatively short times. 



 §.§ Computational considerations

We will use n_r, n_i, and n_g = 2n_r + n_i to indicate the number of grid points in real time, imaginary time, and in total respectively. If n is the number of orbitals, the most expensive computational steps scale like O(n_gn^6) and O(n_g^2n^4). The O(n_gn^6) step arises from the most expensive CCSD contractions which must be performed for the amplitudes at each point in time. The (n_g^2n^4) step comes from the numerical integration in that the amplitudes at each point must be computed from summing quantities at all previous time points. For quadrature rules that strictly obey the integral property

    ∫_0^t_ydt'[…] = ∫_0^t_y - 1dt'[…] + ∫_t_y - 1^t_ydt'[…]

in discretized form,

    G_x^yI_x = G_x^y - 1I_x + G_y^yI_y,

the integrals to t_y can reuse the information from the previous point and the O(n_g^2n^4) step can be reduced to O(n_gn^4). 

The storage requirement scales like O(n_gn^4) if we store all the amplitudes as in FT-CCSD. Most of the amplitudes can be stored on disk at any given time, but the storage requirement becomes quickly prohibitive if n_g is large and this limits the lengths of real time that we can consider. Using the differential form of the equations described in Section <ref> can reduce the storage requirement to O(n^4) and eliminate the O(n_g^2) computational steps.

All computations must use complex arithmetic and, as a result, the grand potential and response properties may acquire an imaginary part. The grand potential and equilibrium properties should
be real for real-valued reference orbitals and exact integration because no operator appears on the forward and backward parts of the Keldysh contour, and any complex exponentials will cancel. 
However, this is not the case for dynamic properties,
and they can acquire an imaginary part due to errors in the coupled cluster approximation on the real axis. This is similar to the small imaginary parts of properties obtained from ground state CCSD with complex orbitals. Thus, if we use reference orbitals that are real, then we may attribute any imaginary part in the grand potential (Keldysh cumulant generating function) to  numerical integration, while for dynamic properties the imaginary part can be due to the numerical integration and to the coupled cluster approximation itself.



§ APPLICATIONS




 §.§ A simple benchmark

To demonstrate that the Keldysh CCSD method yields exact dynamics for a 2 orbital system, we present some benchmark calculations on a Hamiltonian of the form

    H(t)    = ∑_ij h_ija^†_ia_j + 1/2∑_ijkl
        V_ijklc^†_ia^†_ja_la_k 
       + ∑_ij
        d_ija^†_ia_j sin(ω t)

where i,j run over 2 spatial orbitals having the same spin degree of freedom.
We take h and V to be the 1-electron and 2-electron spatial-orbital integrals of a STO-3G hydrogen molecule at R = 0.6Å, and d to be the integrals of the dipole operator. We choose a laser frequency in the UV (ω = 0.2095588) and a temperature of k_BT = 1. In Figure <ref>, we show the error in the dipole moment as a function of time relative to the full CI value obtained by explicitly propagating the full CI, grand canonical density matrix.

As expected, the theory is exact in that the only error is due to the numerical integration. An analysis of the error on a log-log plot indicates that the error is approximately quadratic with time for a constant grid spacing. 



 §.§ Hubbard model with Peierls substitution

The Hubbard model is one of the simplest models of correlated electrons. Despite its simplicity, the Hubbard model is valuable as a testing ground for theory and as a minimal model of correlation-induced phenomena. The interaction of the underlying lattice with a field can be
  incorporated into such a model by adding a “Peierls phase” so that the Hamiltonian becomes time-dependent:

    H(t) = -t_H∑_iσ[ e^iA(t)a_iσ^†a_(i+1)σ 
        + h.c.] + U∑_in_i↑n_i↓.

The parameters t_H and U are the hopping and on-site repulsion respectively. We study the 1-dimensional Hubbard model with a Peierls phase of the form:

    A(t) = A_0 e^-(t - t_0)^2/2σ^2cos[ω(t - t_0)].

This functional form is taken almost exactly from Ref. Wang2017, where the dynamics of this model are studied at zero temperature.
We use U = 0.5, σ = 0.8, t_0 = 2, ω = 6.8, and we show full CI and Keldysh-CCSD results for A_0 = 0.5,1.0,2.0. The simulation starts from a temperature of T = 1.0 at half-filling. We use zero-temperature unrestricted Hartree-Fock (UHF) orbitals and orbital energies for the coupled cluster calculation. These calculations were done with 800 grid points which is enough such that the error due to numerical integration is less than 1%. The results for the difference in population between the two sites are shown in Figure <ref>.

Even for this simple problem, Keldysh-CCSD is no longer exact, but it still performs quite well relative to the exact the result. As expected, the error is larger for larger values of A_0, when the system is driven further out of equilibrium. For A_0 = 0.5 and A_0 = 1.0, the deviations from the exact dynamics are barely visible in Figure <ref>. Despite the fact that global symmetries are not generally conserved by Keldysh-CCSD, the particle number is conserved for this particular Hamiltonian by virtue of the particle-hole symmetry at half-filling.



 §.§ Ultrafast electron dynamics in Si at high temperature

We now apply the method to driven electron dynamics in high-temperature silicon in order to demonstrate clearly the advantages and drawbacks of the method when applied to a more realistic Hamiltonian. We use a single primitive cell of Si in a minimal basis (SZV<cit.> with GTH pseudopotentials<cit.>) 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<cit.>. We use the dipole approximation in the velocity gauge so that the coupling has the form:

    1/mcp·A(t).

A(t) is the time-dependent vector potential and it is related to the electric field by

    E(t) = -1/cȦ(t).


To mimic the application of an ultrafast XUV pulse, we use a pulse of the form

    1/cA(t) = A_0ẑ
        e^-(t - t_0)^2/2σ^2cos[ω(t - t_0)]

with the following pulse width and center: σ = 48.378 as (2 atomic units), t_0 = 362.83 as (15 atomic units). ω is chosen to correspond to a 46.9 nm XUV laser. The system is taken to initially have a temperature of T = 63155.01 K (k_BT = 0.2 E_h). These parameters are unrealistic in two ways. We have chosen an unrealistically short laser pulse so that we can avoid doing dynamics for very long times, and we have chosen an unrealistically high temperature so that finite-temperature effects are accentuated. Despite the unphysical conditions, the Hamiltonian itself should capture the qualitative features of ab initio Hamiltonians of crystalline solids.

We must first ask whether or not the particle number is conserved by Keldysh-CCSD dynamics, since we have already stated that this can be a problem. In Figure <ref>, we have plotted the number of electrons per unit cell as a function of real-time.

Clearly the number of electrons is not conserved by Keldysh-CCSD, and the drift in the particle number is qualitatively the same for different laser intensities. We expect this to lead to unphysical results at longer times.

If we are only interested in relatively short times, the laser-driven dynamics can be corrected by examining the difference induced by the laser field. The differences in the populations of the valence and conduction bands induced by the laser as a function of t are shown in Figure <ref>.

Here we see dynamics that is physically reasonable. At high intensity, the laser increases the population of the conduction band and
depopulates the valence band. Thus despite the violation of particle number conservation, Keldysh-CCSD can provide an approximation to the dynamics which yields qualitatively correct correlated many-body dynamics at finite temperature for short times. 



§ CONCLUSIONS AND FUTURE WORK

In this work, we have shown how the Keldysh formalism may be used to extend imaginary time finite temperature coupled cluster to treat the finite temperature dynamics of correlated fermionic systems. We have discussed the theory in detail and presented results for the singles and doubles variant of this method, Keldysh-CCSD. In addition to the integral equations that form the basis of our current implementation, we have shown how the theory may be extended to the computation of correlation functions with arbitrary time-orderings, and we have shown how a differential form, more natural for dynamics, may be used to avoid storing the amplitudes at all times. One difficulty in applying the method to realistic systems is that there is the potential for the violation of global symmetries. In Section <ref> we have discussed this issue, and in Appendix <ref>, we have presented a more detailed analysis for an exactly solvable system.

In order to demonstrate some of the potential applications and issues with the method, we have performed calculations on a 2-site Hubbard model with a Peierls phase and simple ab initio model of bulk silicon in the presence of an intense XUV laser pulse. For the small Hubbard model, Keldysh-CCSD performs well relative to the exact dynamics, while for the silicon model, we are limited to studying relatively short-time dynamics due to unphysical results at longer times.

The generalization of coupled cluster theory to finite-temperature dynamics presented here displays a number of practical and theoretical challenges which present opportunities for future work. Chief among them are

    
  * Practical improvements necessary to decrease the computational cost and memory requirements
    
  * Modifications of the theory to lessen the impact of particle-number symmetry breaking and other unphysical effects at longer times 



This work is supported by the US Department of Energy, Office of Science, via grant number SC0018140.




§ INTRODUCTION TO THE KELDYSH FORMALISM

We imagine a system in thermal equilibrium at time t = 0. At t = 0 some perturbation, v(t), is turned on. The average value of O at some time t_f is given by

    O(t_f) = Tr[e^-iH_v(t_f)ρ e^iH_v(t_f) O]

where ρ is the equilibrium density operator,

    ρ = e^-β(H - μ N),

for inverse temperature β and chemical potential μ, and H_v(t) is the full Hamiltonian including the perturbation,

    H_v(t) ≡ H + v(t).

The cyclic property of the trace allows us to write this expression as

    O(t_f) = Tr[ρ O_H(t_f)],

where the subscript H indicates the Heisenberg representation.

We can now rewrite this expression in the interaction picture by splitting the Hamiltonian as

    H_v(t) = K_0 + V(t)

where K_0 is the one-electron part,

    K_0 = H_0 - μ N,

and V(t) is the interaction and includes v(t) for t > 0. Note that if the Hamiltonian commutes with the number operator the chemical potential term can be included in the propagator from the beginning or trivially added at any stage. The expression for O becomes

    O(t_f) = 
    	1/ZTr[e^-β K_0
    	U_I(β)U_I^†(t_f)O_I(t_f)U_I(t_f)]

where U_I is the propagator in the interaction picture (as indicated by the “I" subscript), and Z is the grand canonical partition function. Inserting the time-ordered representations of the propagators, we obtain


    O(t_f)   = 1/ZTr[e^-β K_0𝒯exp(-∫_0^βdτ V_I(τ) )
    	𝒯^∗exp(-i∫_t_f^0dt V_I(t))O(t_f)
    	𝒯exp(-i∫_0^t_fdt V_I(t))] 
       = 1/ZTr{e^-β K_0𝒯_C[exp(-i∫_C dt V_I(t))O(t_f)]
    	}.


In the first line we have used the 𝒯 as the time-ordering operator in real or imaginary time and 𝒯^∗ as the anti-time-ordering operator. In the second line we have introduced the contour-ordering operator 𝒯_C which orders its arguments along the Keldysh (or Kadanoff-Baym) contour (see Figure <ref>) and allows us to write the out-of-equilibrium thermal average in a compact way. The contour C usually has three branches that go from 0→ t_f, t_f → 0, and 0→ -iβ respectively, though many different choicwes are possible.



§ KELDYSH CCSD AMPLITUDE EQUATIONS

The amplitude equations are determined by the kernels S_1 and S_2, which contain contractions identical to those in the CCSD T_1 and T_2 equations:


    S_i^a(t)    = f_ai(1 - n_a)n_i + ∑_b f_ab(1 - n_a)s_i^b(t) - ∑_j f_jin_is_j^a(t) + 
    	∑_jb⟨ja||bi⟩(1 - n_a)n_is_j^b(t) + ∑_jbf_jbs_ij^ab(t)
       + 1/2∑_jbc⟨aj||bc⟩(1 - n_a)s_ij^bc(t)
    	- 1/2∑_jkb⟨jk||ib⟩n_is_jk^ab(t)
    	-∑_jb f_jbs_i^b(t)s_j^a(t) 
       + ∑_jbc⟨ja||bc⟩(1 - n_a)
    	s_j^b(t)s_i^c(t) - ∑_jkb⟨jk||bi⟩n_is_j^b(t)s_k^a(t) 
    	-1/2∑_jkbc⟨jk||bc⟩s_i^b(t)s_jk^ac(t)
       - 1/2∑_jkbc⟨jk||bc⟩s_j^a(t)s_ik^bc(t)
    	+ ∑_jkbc⟨jk||bc⟩s_j^b(t)s_ki^ca(t)
    	+ ∑_jkcd⟨jk||bc⟩s_i^b(t)s_j^c(t)s_k^a(t)


    S_ij^ab(t)    = ⟨ab||ij⟩(1 - n_a)(1-n_b)n_jn_j + P_ij∑_c⟨ab||cj⟩s_i^c(t)(1 - n_a)(1 - n_b)n_j 
       - P_ab∑_k⟨kb||ij⟩(1 - n_b)n_in_js_k^a(t) + P_ab∑_cf_bc(1 - n_b)s_ij^ac(t) 
       - P_ij∑_kf_kjn_js_ik^ab(t) + 
    	1/2∑_cd⟨ab||cd⟩(1 - n_a)(1 - n_b)s_ij^cd(t) 
       +1/2∑_kl⟨kl||ij⟩n_in_js_kl^ab(t)
    	+ P_ijP_ab∑_kc⟨kb||cj⟩(1 - n_b)n_js_ik^ac(t)
       + 1/2P_ij∑_cd⟨ab||cd⟩(1 - n_a)(1 - n_b)s_i^c(t)s_j^d(t)
    	+ 1/2P_ab∑_kl⟨kl||ij⟩n_in_js_k^a(t)s_l^b(t) 
    	
       - P_ijP_ab∑_kc⟨ak||cj⟩(1 - n_a)n_js_i^c(t)s_k^b(t)
    	-P_ij∑_kcf_kcs_i^c(t)s_kj^ab(t)
    	- P_ab∑_kcf_kcs_k^as_ij^cb(t)
       + P_ab∑_kcd⟨ka||cd⟩(1 - n_a)s_k^c(t)s_ij^db(t)
    	- P_ij∑_klc⟨kl||ci⟩n_is_k^c(t)s_lj^ab(t) 
    	+ P_ijP_ab∑_kcd⟨ak||cd⟩(1 - n_a)s_i^c(t)s_kj^db(t) 
       - P_ijP_ab∑_klc⟨kl||ic⟩n_is_k^a(t)s_lj^cb(t) 
    	+ 1/2P_ij∑_klc⟨kl||cj⟩n_js_i^c(t)s_kl^ab(t)
    	- 1/2P_ab∑_kcd⟨kb||cd⟩(1 - n_b)s_k^a(t)s_ij^cd(t)
       + 1/4∑_klcd⟨kl||cd⟩s_ij^cd(t)s_kl^ab(t)
    	+ 1/2P_ijP_ab∑_klcd⟨kl||cd⟩s_ik^ac(t)s_lj^db(t)
    	- 1/2P_ab∑_klcd⟨kl||cd⟩s_kl^ca(t)s_ij^db(t)
       - 1/2P_ij∑_klcd⟨kl||cd⟩s_ki^cd(t)s_lj^ab(t)
    	-1/2P_ijP_ab∑_kcd⟨kb||cd⟩(1 - n_b)s_i^c(t)s_k^a(t)s_j^d(t)
    	+ 1/2P_ijP_ab∑_klc⟨kl||cj⟩n_js_i^c(t)s_k^a(t)s_l^b(t) 
       +1/4P_ij∑_klcd⟨kl||cd⟩s_i^c(t)s_j^d(t)s_kl^ab(t)
    	+ 1/4P_ab∑_klcd⟨kl||cd⟩s_k^a(t)s_l^b(t)s_ij^cd(t)
    	-P_ijP_ab∑_klcd⟨kl||cd⟩s_i^c(t)s_k^a(t)s_lj^db(t)	
       - P_ij∑_klcd⟨kl||cd⟩s_k^c(t)s_i^d(t) s_lj^ab(t)
    	- P_ab∑_klcd⟨kl||cd⟩s_k^c(t)s_l^a(t)s_ij^db(t)
    	+1/4P_ijP_ab∑_klcd⟨kl||cd⟩
    	s_i^c(t)s_k^a(t)s_l^b(t)s_j^d(t).


We have used n_p to indicate the Fermi-Dirac occupation number of the pth orbital:

    n_p ≡1/e^β(ε_p - μ) + 1.

Here and elsewhere, we use the notation:

    P_abF[ab] = F[ab] - F[ba].




§ KELDYSH CCSD LAMBDA EQUATIONS

In practice, we solve for the discretized Lagrange multipliers, λ(t_y), by taking derivatives of the discretized Lagrangian with respect to the amplitudes at each grid point. This leads to equations of the form

    λ^f_μ(t_x)    = -L_μ[s_f(t_x),λ̃_f(t_x)] 
    λ^b_μ(t_x)    = -L_μ[s_b(t_x),λ̃_b(t_x)] 
    λ^i_μ(τ_x)    = -L_μ[s_i(τ_x),λ̃_i(τ_x)]

where we have again used f,b,i to indicate forward, backward, and imaginary time branches. The L^μ kernel is local in time and contains precisely those contractions that appear in the ground state CC λ equations with the slight modification that those terms arising from the derivative of the energy are multiplied by a negative sign:


    L_a^i(t)    = -f_ia + ∑_b λ̃^i_b(t)(1 - n_b)f_ba - ∑_j λ̃^j_a(t)n_jf_ij
    	+ ∑_jbλ̃^j_b(t)(1 - n_b)n_j⟨bi||ja⟩ - ∑_jb⟨ij||ab⟩s_j^b(t)
    	
       - ∑_jbλ̃^j_a(t) f_ibs_j^b(t) - ∑_jkbλ̃^i_b(t) f_jas_j^b(t) + ∑_jbcλ̃^i_c(t)(1 - n_c)⟨cj||ab⟩s_j^b(t)
    	 -∑_jkbλ̃^k_a(t)n_k⟨ij||kb⟩s_j^b(t)
       + ∑_jbcλ̃^j_c(t)(1 - n_c)⟨ci||ba⟩s_j^b(t) - ∑_jkbλ̃^k_b(t)n_k⟨ji||ka⟩s_j^b(t)- 1/2∑_jkbcλ̃^j_a(t)n_k⟨ik||bc⟩s_jk^bc(t) 
       - 1/2∑_jkbcλ̃^i_b(t)⟨jk||ac⟩
    	s_jk^bc(t) + ∑_jkbcλ̃^j_b(t)⟨ki||ca⟩s_jk^bc(t)
    	- ∑_jkbcλ̃^j_a(t)⟨ik||bc⟩s_j^b(t)s_k^c(t) 
       - ∑_jkbcλ̃^i_b(t)⟨jk||ac⟩
    	s_j^b(t)s_k^c(t) - ∑_jkbcλ̃^k_b(t)⟨ji||ca⟩s_j^b(t)s_k^c(t)
    	+ 1/2∑_jbcλ̃_cb^ij(t)(1 - n_b)(1 - n_c)n_j⟨cb||aj⟩
       - 1/2∑_jkbλ̃_ab^kj(t)(1 - n_b)n_jn_k⟨ib||kj⟩
    	- ∑_jkbcλ̃^jk_ac(t)(1 - n_c)n_k⟨ic||bk⟩
    	s_j^b(t) - ∑_jkbcλ̃^ik_bc(t)(1 - n_c)n_k⟨jc||ak⟩s_j^b(t) 
       + 1/2∑_jbcdλ̃^ij_cd(t)(1 - n_b)(1 - n_c)⟨cd||ab⟩s_j^b(t)
    	+ 1/2∑_jklbλ̃^kl_ab(t)n_kn_l⟨ij||kl⟩s_j^b(t)
    	-1/2∑_jkbcλ̃^jk_ba(t)f_ics_jk^bc(t)
       - 1/2∑_jkbcλ̃^ji_bc(t)f_kas_jk^bc(t)
    	+1/2∑_jkbcdλ̃^jk_bd(t)(1 - n_d)⟨di||ca⟩s_jk^bc(t)
    	-1/2∑_jklbcλ̃^jl_bc(t)n_l⟨ki||la⟩s_jk^bc(t)
    	
       + ∑_jkbcdλ̃^ji_bd(t)(1 - n_d)⟨kd||ca⟩s_jk^bc(t)
    	- ∑_jklbcλ̃^jl_ba(t)n_l⟨ki||cl⟩s_jk^bc(t)
    	- 1/4∑_jkbcdλ̃^jk_ad(t)(1 - n_d)⟨id||bc⟩t_jk^bc(t)
    	
       +1/4∑_jklbcλ̃^il_bc(t)n_l⟨jk||al⟩s_jk^bc(t)
    	-∑_jkbcdλ̃^ik_db(t)(1 - n_d)⟨dj||ac⟩s_j^b(t)s_k^c(t)
    	+ ∑_jklbcλ̃^lk_ab(t)n_l⟨ij||lc⟩s_j^b(t)s_k^c(t)
       - 1/2∑_jkbcdλ̃^jk_ad(t)(1 - n_d)⟨id||bc⟩s_j^bs_k^c(t)
    	+ 1/2∑_jklbcλ̃^il_bc(t)⟨jk||ad⟩s_j^b(t)s_k^c(t)
    	-1/2∑_jklbcdλ̃^kl_ca(t)⟨ij||db⟩s_j^b(t)s_kl^cd(t)
       - 1/2∑_jklbcdλ̃^ki_cd(t)⟨lj||ab⟩s_j^b(t)s_kl^cd(t)
    	- ∑_jklbcdλ̃^jl_ad(t)⟨ik||bc⟩s_j^b(t)s_kl^cd(t)
    	-∑_jklbcdλ̃^il_bd(t)⟨jk||ac⟩s_j^b(t)s_kl^cd(t)
       +1/4∑_jklbcdλ̃^kl_ab(t)n_k⟨ij||cd⟩s_j^b(t)s_kl^cd(t)
    	+1/4∑_jklbcdλ̃^ij_cd(t)⟨kl||ab⟩s_j^b(t)s_kl^bc(t)
    	-1/2∑_jklbcdλ̃^kl_cb(t)⟨ji||da⟩s_j^b(t)s_kl^cd(t)
    	
       -1/2∑_jklbcdλ̃^kj_cd(t)⟨li||ba⟩s_j^b(t)s_kl^cd(t)
    	+ 1/2∑_jklbcdλ̃^jl_ac(t)
    	⟨ik||bd⟩s_j^b(t)s_k^c(t)s_l^d(t) + 1/2∑_jklbcdλ̃^il_bc(t)
    	⟨jk||ad⟩s_j^b(t)s_k^c(t)s_l^d(t)


    L_ab^ij(t)    = -⟨ij||ab⟩ + P_ijP_abf_iaλ̃^j_b(t) + 
    	P_ij∑_c λ̃^i_c(t)(1 - n_c)⟨cj||ab⟩ - P_ab∑_kλ̃^k_a(t)n_k⟨ij||kb⟩
       + P_ijP_ab∑_kcλ̃^j_b(t)⟨ik||ac⟩
    	s_k^c(t) - P_ij∑_kcλ̃^i_c(t)⟨kl||ab⟩s_k^c(t) - P_ab∑_kcλ̃^k_a⟨ij||cb⟩s_k^c(t) 
       + P_ab∑_c λ̃^ij_ac(t)(1 - n_c)f_cb 
    	- P_ij∑_kλ̃^ik_ab(t)n_kf_jk 
    	+ 1/2∑_cdλ̃^ij_cd(t)(1 - n_c)(1 - n_d) ⟨cd||ab⟩ + 
    	1/2∑_klλ̃^kl_ab(t)n_kn_l⟨ij||kl⟩
       + P_ijP_ab∑_kcλ̃^ik_ac(t)(1 - n_c)n_k⟨cj||kb⟩
    	-P_ij∑_kcλ̃^ik_ab(t)f_jcs_k^c(t)
    	- P_ab∑_kcλ̃^ij_ac(t)f_kbs_k^c(t)
       + P_ab∑_kcdλ̃^ij_ad(t)(1 - n_d)⟨dk||bc⟩s_k^c(t) - 
    	P_ij∑_klcλ̃^il_ab(t)n_l⟨jk||lc⟩s_k^c(t)
    	+ P_ijP_ab∑_kcdλ̃^ik_ad(t)(1 - n_d)⟨dj||cb⟩s_k^c(t) 
       - P_ijP_ab∑_klcλ̃^il_ac(t)n_l⟨kj||lb⟩s_k^c(t) 
    	- ∑_kcdλ̃^ij_cd(t)(1 - n_d)⟨kd||ab⟩s_k^c(t)
    	+ ∑_klcλ̃^kl_ab(t)⟨ij||cd⟩s_k^c(t) 
       - P_ij1/2∑_klcdλ̃^ik_ab(t)
    	⟨jl||cd⟩s_kl^cd (t)- P_ab1/2∑_klcdλ̃^ij_ac(t)
    	⟨kl||bd⟩s_kl^cd(t) + P_ijP_ab∑_klcdλ̃^ik_ac(t)
    	⟨lj||db⟩s_kl^cd(t)
       - P_ab1/2∑_klcdλ̃^kl_ca(t)⟨ij||db⟩s_kl^cd(t)
    	- P_ij1/2∑_klcdλ̃^ki_cd(t)⟨lj||ab⟩s_kl^cd(t)
    	+ 1/4∑_klcdλ̃^kl_ab(t)⟨ij||cd⟩s_kl^cd(t)
    	
       + 1/4∑_klcdλ̃^ij_cd(t)⟨kl||ab⟩s_kl^cd(t)
    	-P_ij∑_klcdλ̃^ik_ab(t)
    	⟨jl||cd⟩s_k^c(t)s_l^d(t) - P_ab∑_klcdλ̃^ij_ac(t)⟨kl||bd⟩s_k^c(t)s_l^d(t) 
       - ∑_klcdλ̃^ik_ad(t)⟨lj||cb⟩s_k^c(t)s_l^d(t)
    	+ 1/2∑_klcdλ̃^kl_ab(t)⟨ij||cd⟩s_k^c(t)s_l^d(t)
    	+ 1/2∑_klcdλ̃^ij_cd(t)⟨kl||ab⟩s_k^c(t)s_l^d(t)


The λ̃ are integrals of the λ amplitudes with an exponential factor. We compute them separately along each branch of the contour:

    λ̃^f_μ(t_x)    =  ∑_yg_yλ^i_μ(τ_y)G_x^t_f/g_x
        e^Δ_μ(it_x - τ_y)
        - i∑_yg_yλ^b_μ(t_y)G_x^t_0/g_x
        e^iΔ_μ(t_x - t_y)
       + i∑_y g_y λ^f_μ(t_y)G_x^y/g_xe^iΔ_μ(t_x - t_y)
    λ̃^b_μ(t_x)    = ∑_yg_yλ^i_μ(τ_y)G_x^t_f/g_x
        e^Δ_μ(it_x - τ_y)
        -i∑_y g_y λ^b_μ(t_y)G_x^y/g_xe^iΔ_μ(t_x - t_y)
    λ̃^i_μ(τ_x)    = ∑_y g_y λ^i_μ(τ_y)G_x^y/g_xe^Δ_μ(τ_x - τ_y)

The t_f and t_0 appearing as indices of grid point tensors are those grid points occuring at t = t_f and t = 0 respectively.



§ KELDYSH CCSD 1-RDM

Given the s_μ and λ_μ amplitudes on the Keldysh contour, we can efficiently form the one-electron density matrix as a function of real time so that any one-electron observable can be computed at any real time. For Keldysh-CCSD, the coupled cluster contribution to the density has 4 parts:


    γ_ia(t_y)    = -λ̃^i_a(t_y) 
    γ_ba(t_y)    = -λ̃^i_a(t_y)s_i^b(t_y) - 1/2λ̃^ki_cb(t_y)s^ca_ki(t_y) 
    γ_ji(t_y)    = λ̃_a^j(t_y)s_i^a(t_y) + 1/2λ̃^kj_ca(t_y)s_ki^ca(t_y)
    γ_ai(t_y)    = s_i^a(t_y) - λ̃_b^j(t_y)s_ji^ba(t_y) +λ̃^j_b(t_y)s_i^b(t_y)s_j^a(t_y)
       + 1/2λ̃^jk_bc(t_y)
    	s_i^b(t_y)s_jk^ac(t_y)
    	+ 1/2λ̃^jk_bc(t_y)
    	s_j^a(t_y)s_ik^bc(t_y)

The value of some one-electron observable at time t_y can be computed as

    O(t_y) = O_0 
        + ∑_iaγ_ia(t_y)O_ain_i(1 - n_a) 
        + ∑_baγ_ba(t_y)O_ab(1 - n_a)
        + ∑_ijγ_ji(t_y)O_ijn_j
        + ∑_aiγ_ai(t_y)O_ia


We have used …_0 as the average in the thermal reference. Note that the density matrix as a function of real time can be computed from either the forward or backward parts of the contour. In the limit of exact integration the result will be the same, but in practice there could be differences. Averaging the values obtained from the forward and backward parts of the contour may provide slightly higher accuracy at the cost of computing additional density matrices. We have not explored this possibility and compute the densities from the forward branch of the contour.



§ AN ANALYSIS OF PARTICLE NUMBER CONSERVATION IN A 1-PARTICLE PROBLEM


In this appendix we will show some illuminating analytic results for a 1-particle Hamiltonian. In particular, we will show explicitly that perturbation theories truncated at 2nd order and at 3rd order preserve the particle number, and, though Keldysh coupled cluster singles (CCS) is exact, linearized Keldysh coupled cluster singles (LCCS) is not exact and does not preserve particle number at finite temperature. These results are supported by numerical calculations, and we use the features of this simple problem to infer some properties of more general Keldysh CC models.



 §.§ Preliminaries

We will consider a 1-particle Hamiltonian of the form

    H = h_0 + f(t).

The corresponding grand potential is denoted by Ω, and 
and we will can do perturbation theory in powers of f:

    Ω = Ω^[0] + Ω^[1] + Ω^[2] + …

The value of some observable at time t is given by

    O(t_1) = .∂/∂αΩ[α]|_α = 0

where

    H[α] = h_0 + f(t) + αδ_C(t - t_1)O

as discussed in Section <ref>. For the purpose of this discussion, we will assume the perturbation, f, is independent of time. This will simplify the algebra and make violations of particle-number conservation more obvious.

We will consider several theories based on perturbation theory and coupled cluster expansions. We will not derive these theories in detail, as they can be found elsewhere (see for example Ref. White2018), but they are most easily considered from a diagrammatic perspective. The 0th and 1st order perturbation theory contributions are trivial, and are not considered here. The diagrams representing 2nd, 3rd, and 4th order perturbation theory contributions are shown in Figure <ref>.  


The CCS and LCCS equations are shown diagrammatically in Figure <ref>. 

A cursory analysis of these diagrams indicates that


  * CCS is complete to all orders in perturbation theory and is therefore exact

  * LCCS is complete only to 3rd order in perturbation theory

  * The 4th order terms missing in LCCS are those arising from the 3rd skeleton in Figure <ref>(c).




 §.§ Second order perturbation theory

The free energy at second order can be expressed as a contour integral

    Ω^[2] = i/β∫_C dt ∑_ia f_ia s_i^a(t)^[1]

where s_i^a(t)^[1] are the first order amplitudes, and C is a contour integral along the Keldysh contour. If we allow the perturbation to become a function of a coupling to the number operator acting at a particular point on the Keldysh contour, we get that

    f_pq(α)    = f_pq + iαδ_C(t - t_1)δ_pq
       ⇒  Ω^[2][α] 
    	= i/β∫_C dt ∑_ia f_ia(α) s_i^a(t)^[1]

where the amplitudes are now implicitly a function of α as well. The amplitudes are equal to 

    s_i^a(t)^[1] = -i∫_C(0)^C(t)dt'e^iΔ_i^a(t' - t)f_ai(α)n_i(1 - n_a)

which implies that the free energy as a function of α is

    Ω^[2](α) 
    	   = 1/β∫_C dt ∑_ia f_ia(α)
       ×∫_C(0)^C(t)dt'
    	e^iΔ_i^a(t' - t)f_ai(α)n_i(1 - n_a).

Expanding Ω in α results in

    Ω^[2](α) = Ω^[2]_0 + Ω^[2]_1α + 1/2Ω^[2]_2 α^2 + …

where Ω^[2]_n is the nth derivative wrt α. We find two contributions to the first derivative:


    Ω_1^[2] = 1/β∑_ian_i(1 - n_a)[
    	∫_C dt f_ia∫_C(0)^C(t)dt'
    	e^iΔ_i^a(t' - t)iδ_C(t' - t_1)δ_pq + 
    	∫_C dt iδ_C(t - t_1)δ_pq∫_C(0)^C(t)dt'
    	e^iΔ_i^a(t' - t)f_ai].

The sum over i,a reduces to a single sum because of the delta functions and we obtain

    Ω_1^[2]   = i/β∑_pn_p(1 - n_p)[
    	∫_C dt f_pp∫_C(0)^C(t)dt'
    	δ_C(t' - t_1) + 
    	∫_C dt δ_C(t - t_1)∫_C(0)^C(t)dt'
    	f_pp] 
       = i/β∑_pn_p(1 - n_p)f_pp[
    	∫_C(t_1)^C(-iβ) dt  + ∫_C(0)^C(t_1)dt'] 
       = i/β∑_pn_p(1 - n_p)f_pp[
    	∫_Cdt] = ∑_pn_p(1 - n_p)f_pp.



This derivative is the 2nd order correction to the particle number in that

    N_0 - βΩ_1^[2] = 
    	N_0 - β∑_pn_p(1 - n_p)f_pp

Notice that this correction is non-zero at finite temperature, equal to the FT-MP2 correction to the particle number (cf. Equation 49 of Ref. Santra2017), but independent of time (t_1).



 §.§ Third order perturbation theory

We can perform the same analysis for the 3rd order terms. We must use the second order amplitudes


    s_i^a(t)^[2]   = -i∫_C(0)^C(t) dt' e^iΔ_i^a(t' - t)[∑_b f_ab(α)s_i^b(t')^[1] - ∑_j f_ji(α)s_j^a(t')^[1]]
       = -∫_C(0)^C(t) dt'e^iΔ_i^a(t' - t)[∑_b (1 - n_a)f_ab(α)∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')f_bi(α)n_i(1 - n_b) 
        - ∑_j n_if_ji(α)∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')f_aj(α)n_j(1 - n_a)].

The corresponding contribution to the grand potential is

    Ω^[3] = -i/β∫_C dt ∑_ia f_ia(α)
    	∫_C(0)^C(t) dt'e^iΔ_i^a(t' - t)∑_b (1 - n_a)f_ab(α)
    	∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')f_bi(α)n_i(1 - n_b)


    + i/β∫_C dt ∑_ia f_ia(α)
    	∫_C(0)^C(t) dt'e^iΔ_i^a(t' - t)∑_j n_if_ji(α)∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')f_aj(α)n_j(1 - n_a).

The two terms in the expression for the grand potential are the two 3rd order diagrams (particle and hole-type diagrams) arising from the 3rd order skeleton depicted in Figure <ref>(b).
Each term has 3 factors of f which means that there will be 6 total terms that are 1st order in α:

    (1)    = -i/β∫_C dt ∑_ia iαδ_C(t - t_1)δ_ia∫_C(0)^C(t) dt'e^iΔ_i^a(t' - t)∑_b (1 - n_a)f_ab∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')f_bin_i(1 - n_b)
    
    	(2)    = -i/β∫_C dt ∑_ia f_ia∫_C(0)^C(t) dt'e^iΔ_i^a(t' - t)∑_b (1 - n_a)iαδ_C(t' - t_1)δ_ab∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')f_bin_i(1 - n_b)
    
    	(3)    = -i/β∫_C dt ∑_ia f_ia∫_C(0)^C(t) dt'e^iΔ_i^a(t' - t)∑_b (1 - n_a)f_ab∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')iαδ_C(t” - t_1)δ_bin_i(1 - n_b)
    
    	(4)    = i/β∫_C dt ∑_ia iαδ_C(t - t_1)δ_ia∫_C(0)^C(t) dt'e^iΔ_i^a(t' - t)∑_j n_if_ji∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')f_ajn_j(1 - n_a)
    
    	(5)    = i/β∫_C dt ∑_ia f_ia∫_C(0)^C(t) dt'e^iΔ_i^a(t' - t)∑_j n_iiαδ_C(t' - t_1)δ_ji∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')f_ajn_j(1 - n_a)
    
    	(6)    = i/β∫_C dt ∑_ia f_ia∫_C(0)^C(t) dt'e^iΔ_i^a(t' - t)∑_j n_if_ji∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')iαδ_C(t” - t_1)δ_ajn_j(1 - n_a)

We will now perform algebra similar to that of previous section to simplify terms (1)-(3):

    (1)    = α/iβΔ_p^b∑_pbf_pb f_bp
    	(1 - n_p)n_p(1 - n_b)[t_1 + 1/iΔ_p^b(e^-iΔ_p^bt_1 - 1)]
    
    	(2)    = α/β∑_ia f_iaf_ai(1 - n_a)n_i(1 - n_a)
    	1/(iΔ_i^a)^2[-e^iΔ_i^at_1e^-βΔ_i^a + e^-βΔ_i^a + 1 -e^-iΔ_i^at_1]
    
    	(3)    = α/β∑_pa f_pa(1 - n_a)f_apn_p(1 - n_p)
    	1/iΔ_p^a[-t_1 - 1/iΔ_p^a 
    	-iβ + e^iΔ_i^at_1/iΔ_i^ae^-βΔ_p^a]

Terms (4)-(6) are analogous.

Examining these expressions, we see that each contains terms dependent on t_1 and that the prefactors of (1) and (3) are the same, but the prefactor of (2) is different. We will now use the notation (1)[t_1] to indicate those terms in (1) that are dependent on t_1, and we will show that

    (1)[t_1] + (2)[t_1] + (3)[t_1] = 0.

There is a slight subtlety here in that while we know that (1), (2), and (3) all have well-defined, finite, limits as Δ→ 0, the same is not true of (1)[t_1] for example. So, for the remainder of this section we will assume that Δ is finite and show that all the t_1 dependence cancels in this case. Since we know that (1) + (2) + (3) has a finite limit, this should also be sufficient to show that the t_1-dependence vanishes in the limit as Δ→ 0, but the expressions that we will be be working with will themselves not be valid in this limit. If we wanted to be perfectly rigorous, we would consider this limit separately.

The t_1-dependent parts are

    (1)[t_1]    = α/β∑_paf_pa f_ap
    	(1 - n_p)n_p(1 - n_a)1/iΔ_p^a[t_1 + 
    	1/iΔ_p^ae^-iΔ_p^at_1] 
    
    	(2)[t_1]    = α/β∑_ia f_iaf_ai(1 - n_a)n_i(1 - n_a)
    	1/(iΔ_i^a)^2[-e^iΔ_i^at_1e^-βΔ_i^a -
    	e^-iΔ_i^at_1]
    
    	(3)[t_1]    = α/β∑_pa f_pa(1 - n_a)f_apn_p(1 - n_p)
    	1/iΔ_p^a[-t_1 + 
    	e^iΔ_p^at_1/iΔ_p^ae^-βΔ_p^a].

Each expression has two terms and the 1st term of (1)[t_1] cancels the first term of (3)[t_1] so that the sum of all of them is given by

    (1)[t_1] + (2)[t_1] + (3)[t_1]   = α/β∑_paf_pa f_ap
    	(1 - n_p)n_p(1 - n_a)1/(iΔ_p^a)^2e^-iΔ_p^at_1
       - α/β∑_ia f_iaf_ai(1 - n_a)n_i(1 - n_a)
    	1/(iΔ_i^a)^2e^iΔ_i^at_1e^-βΔ_i^a
       - α/β∑_ia f_iaf_ai(1 - n_a)n_i(1 - n_a)
    	1/(iΔ_i^a)^2e^-iΔ_i^at_1
       + α/β∑_pa f_pa(1 - n_a)f_apn_p(1 - n_p)
    	1/(iΔ_p^a)^2e^iΔ_p^at_1e^-βΔ_p^a

A critical simplification is acheived by multiplying the real exponential factors into the occupation numbers so that, for example, 

    (1 - n_a)n_i(1 - n_a)e^-βΔ_i^a = (1 - n_a)(1 - n_i)n_a.

The results are that

    (1)[t_1] + (2)[t_1] + (3)[t_1]   = α/β∑_paf_pa f_ap
    	(1 - n_p)n_p(1 - n_a)1/(iΔ_p^a)^2e^-iΔ_p^at_1
       - α/β∑_ia f_iaf_ai(1 - n_a)(1 - n_i)n_a
    	1/(iΔ_i^a)^2e^iΔ_i^at_1
       - α/β∑_ia f_iaf_ai(1 - n_a)n_i(1 - n_a)
    	1/(iΔ_i^a)^2e^-iΔ_i^at_1
       + α/β∑_pa f_paf_apn_a(1 - n_p)(1 - n_p)
    	1/(iΔ_p^a)^2e^iΔ_p^at_1

We can now relabel indices so that in the second line a,i → p,a and in the 4th line a,p → i,a:

    (1)[t_1] + (2)[t_1] + (3)[t_1]   = α/β∑_paf_pa f_ap
    	(1 - n_p)n_p(1 - n_a)1/(iΔ_p^a)^2e^-iΔ_p^at_1
       - α/β∑_ap f_apf_pa(1 - n_p)(1 - n_a)n_p
    	1/(iΔ_a^p)^2e^iΔ_a^pt_1
       - α/β∑_ia f_iaf_ai(1 - n_a)n_i(1 - n_a)
    	1/(iΔ_i^a)^2e^-iΔ_i^at_1
       + α/β∑_ai f_aif_ian_i(1 - n_a)(1 - n_a)
    	1/(iΔ_a^i)^2e^iΔ_a^it_1

Everything is now close to canceling completely. We note that

    Δ_p^q = -Δ_q^p

which allows us to write

    (1)[t_1] + (2)[t_1] + (3)[t_1]   = α/β∑_paf_pa f_ap
    	(1 - n_p)n_p(1 - n_a)1/(iΔ_p^a)^2e^-iΔ_p^at_1
       - α/β∑_ap f_apf_pa(1 - n_p)n_p(1 - n_a)
    	1/(iΔ_p^a)^2e^-iΔ_p^at_1
       - α/β∑_ia f_iaf_ai(1 - n_a)n_i(1 - n_a)
    	1/(iΔ_i^a)^2e^-iΔ_i^at_1
       + α/β∑_ai f_aif_ia(1 - n_a)n_i(1 - n_a)
    	1/(iΔ_i^a)^2e^-iΔ_i^at_1.

Finally, it is clear that the first line cancels the second and the 3rd line cancels the 4th. Thus, we can say that (1) + (2) + (3) makes a contribution to the particle number that is independent of t_1. An exactly analogous argument can presumably be made for (4) + (5) + (6). This means that for 3rd order perturbation theory each of the two 3rd order diagrams individually conserves particle number.




 §.§ Some features of 4th order perturbation theory

4th order perturbation theory becomes quite tedious, but it is sufficient for our purposes to consider the pair of 4th order diagrams that arise from the quadratic term in the CCS equations. These diagrams correspond to the skeleton listed 3rd in Figure <ref>(c):


    Ω^[4]∗(α)    = i/β∫_Cdt ∑_ia
    	f_ia(α) s_i^a(t)^[3∗]
       = i/β∫_Cdt ∑_ia
    	f_ia(α) (-i)∫_C(0)^C(t)dt' e^iΔ_i^a(t' - t)∑_jb f_jb(α)s_j^a(t')^[1]s_i^b(t')^[1]

We have indicated this particular pair of 4th order contributions by an asterisk. Plugging in the definition of the 1st order s amplitudes given in Equation <ref> yields:

    Ω^[4]∗(α)    =1/β∫_Cdt ∑_ia
    	f_ia(α) ∫_C(0)^C(t)dt' e^iΔ_i^a(t' - t)∑_jbf_jb(α)
       × (-i)∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')f_aj(α)n_j(1 - n_a) (-i)∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')f_bi(α)n_i(1 - n_b)

Just as before, we expand this expression as a function of α to obtain

    Ω^[4]∗(α) = Ω^[4]∗_0 + Ω^[4]∗_1α + 1/2Ω^[4]∗_2 α^2 + …


We are interested in the coefficient of α in this expansion since it gives the 4th order (∗) correction to the particle number. Clearly there are 4 terms:

    (1)    =1/β∫_Cdt ∑_iaδ_iaδ_C(t - t_1) ∫_C(0)^C(t)dt' e^iΔ_i^a(t' - t)∑_jbf_jb
       × (-i)∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')f_ajn_j(1 - n_a) (-i)∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')f_bin_i(1 - n_b) 
    
    	(2)    =1/β∫_Cdt ∑_ia
    	f_ia∫_C(0)^C(t)dt' e^iΔ_i^a(t' - t)∑_jbδ_jbδ_C(t' - t_1)
       × (-i)∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')f_ajn_j(1 - n_a) (-i)∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')f_bin_i(1 - n_b)
    	
    
    	(3)    =1/β∫_Cdt ∑_ia
    	f_ia∫_C(0)^C(t)dt' e^iΔ_i^a(t' - t)∑_jbf_jb
       × (-i)∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')δ_ajδ_C(t” - t_1)n_j(1 - n_a) (-i)∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')f_bin_i(1 - n_b) 
    
    	(4)    =1/β∫_Cdt ∑_ia
    	f_ia∫_C(0)^C(t)dt' e^iΔ_i^a(t' - t)∑_jbf_jb
       × (-i)∫_C(0)^C(t')dt”e^iΔ_j^a(t” - t')f_ajn_j(1 - n_a) (-i)∫_C(0)^C(t')dt”e^iΔ_i^b(t” - t')δ_biδ_C(t” - t_1)n_i(1 - n_b)


After a great deal of algebra, we obtain

    (1)    =-1/β∑_pjbn_pn_j(1 - n_p)(1 - n_b)f_jbf_bpf_pj1/iΔ_j^biΔ_p^b
       ×[t_1 + 1/iΔ_j^p(e^-iΔ_j^pt_1 - 1)  + 1/iΔ_p^b(e^-iΔ_p^bt_1 - 1) - 1/iΔ_j^b( e^-iΔ_j^bt_1 - 1)]
    
    	(2)    =-1/β∑_iap(1 - n_i)n_pn_a(1 - n_p)f_iaf_apf_pi1/iΔ_p^a1/iΔ_i^p1/iΔ_i^a[1 - e^-iΔ_p^it_1 - e^-iΔ_a^pt_1 + e^iΔ_i^at_1] 
       -1/β∑_iapn_in_p(1 - n_a)(1 - n_p)f_iaf_apf_pi1/iΔ_p^a1/iΔ_i^p1/iΔ_i^a[e^-Δ_i^at_1 - e^-iΔ_p^at_1 - e^-iΔ_i^pt_1 + 1]
    
    	(3)    =-1/β∑_ipbn_in_p(1 - n_p)(1 - n_b)f_bif_ipf_pb1/iΔ_i^b
       ×{1/iΔ_i^p[-iβ - t_1 - 1/iΔ_i^p(e^iΔ_i^pt_1e^-Δ_i^pβ) - 1)] - 1/iΔ_b^p[1/iΔ_b^i(e^βΔ_b^i - e^iΔ_b^it_1) - 1/iΔ_i^p(e^iΔ_b^pt_1e^-βΔ_i^p - 1)]}
    
    	(4)    =-1/β∑_japn_pn_j(1 - n_a)(1 - n_p)f_ajf_jpf_pa1/iΔ_j^a
       ×{1/iΔ_p^a[-iβ - t_1 + 1/iΔ_p^a(e^-βΔ_p^ae^iΔ_p^at_1 - 1)]  
    	- 1/iΔ_j^p[1/iΔ_j^a(e^-βΔ_j^a - e^-iΔ_j^at_1) - 1/iΔ_j^p(e^-βΔ_p^ae^-iΔ_j^pt_1 - 1)]}

These expressions are complicated and difficult to verify, but it is sufficient to examine the coefficients of the terms linear in t_1 to see that for (1) + (2) + (3) + (4), the t_1-dependence does not cancel. Since the sum of all 4th order terms should have no dependence on t_1 as per the argument in Section <ref>, this also implies that Ω^[4] - Ω^[4]∗ does not generate dynamics that conserve particle number. Since LCCS includes the Ω^[4] - Ω^[4]∗ contribution at 4th order, this suggests that Keldysh-LCCS will not conserve particle number.




 §.§ Numerical results

From this discussion and the discussion given in Section <ref>, we conclude several things about particle number conservation in the 1-particle problem:

	
  * nth order perturbation theory should conserve particle number
	
  * For 2nd and 3rd order perturbation theory, each diagram individually conserves the particle number
	
  * But this is not true of 4th order perturbation theory where the (∗) terms (and therefore the remaining terms also) do not conserve particle number
	
  * LCCS does not conserve particle number
	
  * CCS does conserve particle number by virtue of being exact: it is complete to all orders in perturbation theory.

We will now verify numerically these statements. Because there is error in the particle number due to the numerical integration, we will say that a given theory does not conserve particle number if the quantity changes in time and this change does not decrease as the grid is made finer. Our example Hamiltonian is a two-level system,

    H = [ E_1   0;   0 E_2 ] + 
    	[ V_11 V_12; V_21 V_22;      ].

For this example, we use

    E_1 = 0.1      E_2 = 0.4


    V_11 = V_22 = 0.1      V_12 = V_21^∗ = 1 + i/2.

The relatively large perturbation is used to accentuate the defects in the theory. We use μ = 0.0 and T = 0.5 in these calculations, and compute N at points between t = 0 and t = 1. For the exact dynamics of the system, N is conserved. 

In Table <ref>, we show the number of particles for different approximate methods for the equilibrium grand potential. As expected, different many-body methods will lead to different numbers of electrons at finite temperature, and CCS is exact for this 1-particle problem.


In the following time-dependent calculations, separate calculations are performed at each time point and the property of interest, N, is computed by finite difference differentiation. Since separate calculations have to be performed for each time point anyway, we use the same number of grid points at each time point. This leads to an integration error which is easy to recognize and separate from errors in the many-body approximations themselves. In Figures <ref>-<ref> we confirm numerically point (1).



These figures indicate that for perturbation theory at any order, N will be conserved as a function of time in the limit of exact numerical integration.

Point (2) is trivially true at 2nd order where there is only a single diagrammatic contribution. At third order however, there are two diagrams which we refer to as  the “particle” and “hole” diagrams for the diagrams containing more particle or hole propagators respectively. In Figure <ref> we show that each of these contributions individually conserves N.

This confirms numerically the result of Section <ref>.

Point (3) is confirmed numerically in Figure <ref> where the ∗ contributions at 4th order are included and the particle number is no longer conserved.


Point(4) is confirmed in Figure <ref> which qualitatively resembles Figure <ref>. 


Finally, point (5) is corroborated in Figure <ref>.


In this particular example, the fact that particle number is not conserved does not create a significant problem because the differences on this time-scale are small compared to the error in the various methods at equilibrium (see Table <ref>), but this will not be true in general.

In this Appendix, we have examined particle-number conservation of different perturbation theory and coupled cluster models for finite-temperature dynamics without two-particle interactions. The problem is trivially solvable, but it nonetheless allows us to show explicitly conclusions (1) - (5). In Section <ref> we verified these results numerically. This discussion suggests that in more general Keldysh-CC models, such as the Keldysh-CCSD presented in this article, global symmetries will not be preserved by the dynamics. The seriousness of this flaw will, in general, depend on the system.


