Time-dependent coupled cluster theory on the Keldysh contour for non-equilibrium systems

TL;DR

This paper extends finite temperature coupled cluster theory using the Keldysh formalism to study non-equilibrium thermal systems, enabling detailed simulations of driven quantum systems with applications to models like Hubbard and silicon under ultrafast pulses.

Contribution

The authors develop and implement a Keldysh coupled cluster method for non-equilibrium, finite-temperature quantum systems, expanding the applicability of coupled cluster techniques.

Findings

Successfully applied to Hubbard model with Peierls phase

Demonstrated capability on warm-dense silicon under XUV pulse

Extended coupled cluster theory to non-equilibrium thermal conditions

Abstract

We leverage the Keldysh formalism to extend our implementation of finite temperature coupled cluster theory [\textit{J. Chem. Theory Comput.} 2018, \textit{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 {\it ab initio} model of warm-dense silicon subject to an ultrafact XUV pulse.

Tables

Table 1. Table 1: The average number of particles, ⟨ N ⟩ delimited-⟨⟩ 𝑁 \left<N\right> , for each method in equilibrium. Note that the small error in CCS is due to the numerical integration and to error in the finite difference determination of ⟨ N ⟩ delimited-⟨⟩ 𝑁 \left<N\right> .

method	$N$	error	% relative error
PT2	0.6679063	-0.20734	23.69
PT3	0.9500828	0.07484	8.55
PT4	1.0446668	0.16942	19.36
PT4$\ast$	1.0580504	0.18281	20.89
LCCS	0.9436346	0.06839	7.81
CCS	0.8752419	-3.6189E-07	4.13E-05
exact	0.8752423	-	-

Equations

$$\Omega ={\Omega }^{(0)}+{\Omega }^{(1)}+{\Omega }_{CC}$$

$${\Omega }_{CC}$$

$$s_{\mu }(\tau )=-{\int }_0^{\tau }d{\tau }^{\prime }{e}^{{\Delta }_{\mu }({\tau }^{\prime }-\tau )}{\text{S}}_{\mu }[\mathbf{s}({\tau }^{\prime })].$$

$$\mathcal{L}\equiv \frac{1}{\beta }{\int }_0^{\beta }\mathrm{E}(\tau )-\frac{1}{\beta }\sum _{\mu }{\int }_0^{\beta }d\tau {\lambda }_{\mu }(\tau )\left[s_{\mu }(\tau )+{\int }_0^{\tau }d{\tau }^{\prime }{e}^{{\Delta }_{\mu }({\tau }^{\prime }-\tau )}{\mathrm{S}}_{\mu }({\tau }^{\prime })\right].$$

$${\lambda }_{\mu }(\tau )=-{\text{L}}_{\mu }[\mathbf{s}(\tau ),\tilde{\lambda }(\tau )].$$

$${\tilde{\lambda }}_{\mu }(\tau )={\int }_{\tau }^{\beta }d{\tau }^{\prime }{e}^{{\Delta }_{\mu }(\tau -{\tau }^{\prime })}{\lambda }_{\mu }({\tau }^{\prime }).$$

$$\Omega =-\frac{1}{\beta }\ln \text{Tr}\left[e^{-\beta K_0}\mathcal{T}\exp \left(-i{\int }_0^{-i\beta }dt{V}_I(t)\right)\right].$$

$$\Omega =-\frac{1}{\beta }\ln \text{Tr}\left[e^{-\beta K_0}{\mathcal{T}}_C\exp \left(-i{\int }_{C}dt{V}_I(t)\right)\right].$$

$$\Omega [\alpha ]=-\frac{1}{\beta }\ln \text{Tr}\left[e^{-\beta K_0}{\mathcal{T}}_C\exp \left(\alpha O(t_f)-i{\int }_{C}dt{V}_I(t)\right)\right].$$

$$\Omega [0]=\Omega ,$$

$${\frac{\partial \Omega }{\partial \alpha }|}_{\alpha =0}$$

$$\displaystyle\times\mathcal{T}_{C}\left[\exp\left(-i\int_{C}dtV_{I}(t)\right)O(t_{f})\right]\bigg{\}}.$$

$$⟨O(t_f)⟩=-\beta {\frac{\partial \Omega [\alpha ]}{\partial \alpha }|}_{\alpha =0}$$

$$\Omega [\alpha ]=-\frac{1}{\beta }\ln \text{Tr}\left[e^{-\beta K_0}{\mathcal{T}}_C\exp \left(-i{\int }_{C}dt[V_I(t)+i\alpha O(t){\delta }_C(t-t_f)]\right)\right].$$

$$V_I(t)+i\alpha O(t){\delta }_C(t-t_f),$$

$${\Omega }_{CC}$$

$$s_{\mu }^f(t)$$

$$s_{\mu }^b(t)$$

$$s_{\mu }^i(\tau )$$

$$\mathcal{L}\equiv \frac{i}{\beta }{\int }_{C}dt\mathrm{E}(t)-\frac{i}{\beta }{\int }_{C}dt{\lambda }^{\mu }(t)\left[s_{\mu }(t)+{\int }_{C(0)}^{C(t)}d{t}^{\prime }{e}^{i{\Delta }_{\mu }(t^{\prime }-t)}{\mathrm{S}}_{\mu }(t^{\prime })\right].$$

$$⟨O_1(t_1)O_2(t_2),\dots ⟩$$

$$\Omega [{\alpha }_1,{\alpha }_2,\dots ]=-\frac{1}{\beta }\ln \text{Tr}\left[e^{-\beta K_0}{\mathcal{T}}_C\exp \left(-i{\int }_{C}dt[V_I(t)+i{\alpha }_1{O}_1(t){\delta }_C(t-t_1)+i{\alpha }_2{O}_2(t){\delta }_C(t-t_2)+\dots ]\right)\right]$$

$$⟨O_1(t_1)O_2(t_2),\dots ,O_n(t_n)⟩={\beta }^n{\frac{{\partial }^n\Omega }{\partial {\alpha }_1\partial {\alpha }_2\dots }|}_{{\alpha }_i=0}.$$

$$s_{\mu }(t)$$

$${\tilde{\lambda }}_{\mu }(t)$$

$$\frac{d{s}_{\mu }}{dt}$$

$$\frac{d{\tilde{\lambda }}_{\mu }}{dt}$$

$$\Omega [\lambda ]=\Omega [0]+\lambda {\Omega }^{\prime }[0]+\frac{{\lambda }^2}{2}{\Omega }^{\prime \prime }[0]+\dots$$

$${\Omega }^{[n]}={\frac{1}{n!}\frac{{\partial }^n\Omega [\lambda ]}{\partial {\lambda }^n}|}_{\lambda =0}=\sum _x{c}_x\Omega [{\lambda }_x]$$

$${⟨O⟩}^{[n]}(t)$$

$$=\frac{\partial }{\partial \alpha (t)}{\sum _x{c}_x\Omega [{\lambda }_x,\alpha (t)]|}_{\alpha (t)=0}.$$

$$s_{\mu }^f(t_y)$$

$$s_{\mu }^b(t_y)$$

$$+i\sum _x{G}_x^y{e}^{i{\Delta }_{\mu }(t_x-t_y)}{\text{S}}_{\mu }[{\mathbf{s}}^b(t_x)]$$

Full text

I 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.Kübler et al. (2007); Lu et al. (2008); Fausti et al. (2011); Kim et al. (2012); Morrison et al. (2014) However, the direct simulation of the dynamics is complicated by several factors including:

The quantum many-body problem for a realistic material 2. 2.

The statistical mixture of many-body states implied by a finite temperature 3. 3.

The generically non-thermal, time-dependent distribution of states generated by the perturbation 4. 4.

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 moleculesMarques and Gross (2004); Botti et al. (2007); Otobe, Yabana, and Iwata (2009); Burgdörfer et al. (2014); Sato et al. (2015); Ullrich, Reinhard, and Suraud (1998); Tong and Chu (2001); Baer et al. (2003); Marques et al. (2003); Cheng, Evans, and Van Voorhis (2006), while time-dependent ab initio wavefunction methods such as time-dependent configuration interaction Krause, Klamroth, and Saalfrank (2005); Schlegel, Smith, and Li (2007); Greenman et al. (2010) and coupled cluster methods have been introduced in moleculesSchönhammer and Gunnarsson (1978); Hoodbhoy and Negele (1978, 1979); Huber and Klamroth (2011); Kvaal (2012); Nascimento and Deprince (2017). 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 diagonalizationChen et al. (2009), time-dependent density matrix renormalization group methodsCazalilla and Marston (2002); Verstraete, Garcia-Ripoll, and Cirac (2004); Kokalj and Prelovšek (2009); Wolf, McCulloch, and Schollwöck (2014); Ren, Shuai, and Kin-Lic Chan (2018), diagrammatic Monte CarloWerner, Oka, and Millis (2009); Schiró and Fabrizio (2009); Segal, Millis, and Reichman (2010); Werner et al. (2010); Antipov et al. (2017), and hierarchical equations-of-motionJin, Zheng, and Yan (2008); Huo and Coker (2011); Schinabeck, Härtle, and Thoss (2018). 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,Cizek and Paldus (1980); Bartlett (1981); Crawford and Schaefer III (2000); Shavitt and Bartlett (2009) and where extensions to periodic solidsHirata et al. (2001, 2004); Booth et al. (2013); McClain et al. (2017), finite temperaturesSanyal, Mandal, and Mukherjee (1992); Mandal et al. (2003); White and Chan (2018); Hummel (2018), and open systems Dzhioev and Kosov (2015) 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 spaceSeminoff and Umezama (1983); Matsumoto et al. (1983); Matsumoto and Umezawa (1985). We work within the density matrix language, and use the Keldysh formalismSchwinger (1961); Kadanoff and Baym (1962); Keldysh (1965) to provide a straightforward generalization of the imaginary-time finite-temperature CC equations in Ref. White and Chan, 2018 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.

II Theory

We will describe the theory in 4 steps:

In Section II.1 we describe the working equations of the FT-CC method presented in Ref. White and Chan, 2018 2. 2.

In Section II.2 we describe some important aspects of the Keldysh formalism 3. 3.

In Section II.3 we generalize the imaginary time FT-CC to out of equilibrium systems using the Keldysh formalism 4. 4.

In the final Sections (II.4-II.6) we discuss several important properties and extensions of the theory.

II.1 Finite temperature coupled cluster

In Ref. White and Chan, 2018 we presented a finite temperature coupled cluster theory which provides an ansatz for the correlation contribution to the grand potential $\Omega$ such that

[TABLE]

where $\Omega^{(0)}+\Omega^{(1)}$ is the grand potential of the thermal Hartree-Fock theory and $\Omega_{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)Sanyal, Mandal, and Mukherjee (1992); Sanyal et al. (1993); Mandal, Ghosh, and Mukherjee (2001); Mandal et al. (2003) (the relationship between the thermal Wick expansion used in TCC and thermofield dynamics is also discussed in Mukherjee’s work Sanyal et al. (1993)). The coupled cluster correction to the grand potential is computed from time-dependent amplitudes, $s_{\mu}(\tau)$, evaluated in imaginary time:

[TABLE]

Here, $\beta$ is the inverse temperature, we use $f$ and $\langle pq||rs\rangle$ 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 $\mu$ defines the truncation of the theory, and the amplitudes are determined from a set of non-linear Volterra integral equations of the form

[TABLE]

$\Delta_{\mu}$ 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 S1 and S2 kernels are given in the appendix of Ref. White and Chan, 2018.

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

[TABLE]

(Equation 37 of Ref. White and Chan, 2018) with Lagrange multipliers $\lambda_{\mu}(\tau)$ such that $\mathcal{L}$ coincides with the free energy at its stationary point but is a variational function of $\lambda_{\mu}(\tau)$ and $s_{\mu}(\tau)$. Variational optimization of $\mathcal{L}$ with respect to $\lambda_{\mu}$ yields the amplitude equations (Equation 3) and variational optimization with respect to $s_{\mu}$ yields linear equations for $\lambda_{\mu}$ of the form

[TABLE]

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

[TABLE]

The discretized versions of these equations are given in Appendices A and C of Ref. White and Chan, 2018.

II.2 The Keldysh formalism for non-equilibrium systems

Our approach to the non-equilibrium problem is based on the Keldysh formalism.Schwinger (1961); Kadanoff and Baym (1962); Keldysh (1965); van Leeuwen et al. (2006) This approach is conceptually appealing because the 3 generic difficulties mentioned in Section I 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 A.

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 52). In recent years, many-body methods using this formalism have gone by the name “Kadanoff-Baym dynamics”Kwong and Bonitz (2000); Von Friesen, Verdozzi, and Almbladh (2009); Puig Von Friesen, Verdozzi, and Almbladh (2010). Such methods commonly compute the Green’s functions along the Keldysh contour,Kwong and Bonitz (2000); Balze and Bonitz (2013); Aoki et al. (2014); Špička, Velický, and Kalvová (2014); Do (2014) 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:

[TABLE]

$K_{0}$ is the zeroth order Hamiltonian including the chemical potential (Equation 49), 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\rightarrow t_{f}$) and backward ($t_{f}\rightarrow 0$) branches will exactly cancel (see Fig. 1):

[TABLE]

We may now introduce a functional dependence on a parameter $\alpha$ such that

[TABLE]

Notice that for $\alpha=0$ this is just the equilibrium grand potential,

[TABLE]

but the derivative of $\Omega[\alpha]$ evaluated $\alpha=0$ will reproduce Equation 52 with an extra factor of $\beta$:

[TABLE]

This means that the derivatives of $\Omega[\alpha]$, the Keldysh cumulant generating function, are the observables of the non-equilibrium system

[TABLE]

II.3 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 $\Omega[\alpha]$ as

[TABLE]

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,

[TABLE]

where $\delta_{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 $\Omega$ 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_{\mu}(t)$ where the time argument, $t$, runs along the Keldysh contour. The coupled cluster contribution to the grand potential is given by

[TABLE]

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

[TABLE]

where we have used the superscripts $f,b,i$ to refer to the forward, backward, and imaginary parts of the Keldysh contour (see Figure 1). Explicit expressions for the kernel, $\mathrm{S}_{\mu}$, of Keldysh-CCSD are given in Appendix B.

We now define a Lagrangian in analogy with Equation 4:

[TABLE]

We have used $C(t)$ to indicate the limits of the integral along the Keldysh contour. In Appendix C, we discuss the equations determining the multipliers, $\lambda^{\mu}(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 D. 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\leq 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 III.

II.4 Correlation functions from higher-order response

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

[TABLE]

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

[TABLE]

and use the flexibility of the Keldysh formalism to choose the contour and the locations of the contour delta functions ($\delta_{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 $n$th order derivative:

[TABLE]

Therefore higher order responses of a coupled cluster approximation to $\Omega[\alpha_{1},\alpha_{2},\ldots]$, constructed along the appropriate contour, can be used to obtain the coupled cluster approximation to arbitrarily general correlation functions.

II.5 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 Sanyal, Mandal, and Mukherjee, 1992; Sanyal et al., 1993; Mandal, Ghosh, and Mukherjee, 2001; Mandal et al., 2003). We start with the integral equations for $s$ and $\tilde{\lambda}$ along the Keldysh contour:

[TABLE]

Note that since the density (Equations 65-68) depends only on $\tilde{\lambda}$, we never need to work with $\lambda$ directly using this formulation.

Taking the time derivative of Equations 23 and 24 yields equations of motion for $s$ and $\tilde{\lambda}$:

[TABLE]

In both these equations, the term proportional to $\Delta_{\mu}$ 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 $\tilde{\lambda}$-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_{\mu}$ is initially known to be zero at $t=0$, the value of $\tilde{\lambda}_{\mu}$ is known at $t=-i\beta$. Considering the dependence of Lμ on $s$, this requires the calculation to be carried out in two steps. First $s_{\mu}$ is propagated from $t=0$ to $t=-i\beta$. This propagation can be carried out along the Keldysh contour or just along the imaginary axis. Second, both $s$ and $\tilde{\lambda}$ are propagated back from $t=-i\beta$ to $t=0$ along the Keldysh contour and the density matrix is computed along the way as per Appendix D. 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.

II.6 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 BaymBaym and Kadanoff (1961); Baym (1962) 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 nonethelessSchönhammer and Gunnarsson (1978); Hoodbhoy and Negele (1978, 1979); Huber and Klamroth (2011); Kvaal (2012); Nascimento and Deprince (2017). Furthermore, even conserving approximations can yield qualitatively incorrect resultsPuig Von Friesen, Verdozzi, and Almbladh (2010). 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, $\Omega$, that admits a perturbation theory expansion in some interaction controlled by a coupling, $\lambda$, such that

[TABLE]

Note that $n$th order correction can be computed by finite difference in that

[TABLE]

where the $\lambda_{x}$ are points slightly displaced from $\lambda=0$, and the $\{c_{x}\}$, $\{\lambda_{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, $\alpha(t)$, the $n$th order contribution to the observable is

[TABLE]

This implies that if $O$ is globally conserved by the exact dynamics, it should also be conserved in perturbation theory. In Appendix E, 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 IV.

III Implementation

To enable a concrete implementation, we truncate the coupled cluster equations to the singles and doubles level (CCSD), where $\mu$ labels at most four fermionic lines. The implementation of Keldysh-CCSD closely mirrors the implementation of FT-CCSD, so we refer to Ref. White and Chan, 2018 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.

III.1 Numerical integration

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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.

[TABLE]

where $\delta$ 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.

III.2 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_{g}n^{6})$ and $O(n_{g}^{2}n^{4})$. The $O(n_{g}n^{6})$ step arises from the most expensive CCSD contractions which must be performed for the amplitudes at each point in time. The $(n_{g}^{2}n^{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

[TABLE]

in discretized form,

[TABLE]

the integrals to $t_{y}$ can reuse the information from the previous point and the $O(n_{g}^{2}n^{4})$ step can be reduced to $O(n_{g}n^{4})$.

The storage requirement scales like $O(n_{g}n^{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 II.5 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.

IV Applications

IV.1 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

[TABLE]

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 ($\omega=0.2095588$) and a temperature of $k_{B}T=1$. In Figure 2, 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.

IV.2 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:

[TABLE]

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:

[TABLE]

This functional form is taken almost exactly from Ref. Wang et al., 2017, where the dynamics of this model are studied at zero temperature. We use $U=0.5$, $\sigma=0.8$, $t_{0}=2$, $\omega=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 3.

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 3. 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.

IV.3 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 (SZVVandeVondele and Hutter (2007) with GTH pseudopotentialsGoedecker and Teter (1996); Hartwigsen, Goedecker, and Hutter (1998)) 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 fittingSun et al. (2018). We use the dipole approximation in the velocity gauge so that the coupling has the form:

[TABLE]

$\mathbf{A}(t)$ is the time-dependent vector potential and it is related to the electric field by

[TABLE]

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

[TABLE]

with the following pulse width and center: $\sigma=48.378$ as (2 atomic units), $t_{0}=362.83$ as (15 atomic units). $\omega$ 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_{B}T=0.2E_{\text{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 4, 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 5.

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.

V 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 II.6 we have discussed this issue, and in Appendix E, 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 2. 2.

Modifications of the theory to lessen the impact of particle-number symmetry breaking and other unphysical effects at longer times

Acknowledgements.

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

Appendix A 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

[TABLE]

where $\rho$ is the equilibrium density operator,

[TABLE]

for inverse temperature $\beta$ and chemical potential $\mu$, and $H_{v}(t)$ is the full Hamiltonian including the perturbation,

[TABLE]

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

[TABLE]

where the subscript $H$ indicates the Heisenberg representation.

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

[TABLE]

where $K_{0}$ is the one-electron part,

[TABLE]

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

[TABLE]

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

[TABLE]

In the first line we have used the $\mathcal{T}$ as the time-ordering operator in real or imaginary time and $\mathcal{T}^{\ast}$ as the anti-time-ordering operator. In the second line we have introduced the contour-ordering operator $\mathcal{T}_{C}$ which orders its arguments along the Keldysh (or Kadanoff-Baym) contour (see Figure 1) 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\rightarrow t_{f}$, $t_{f}\rightarrow 0$, and $0\rightarrow-i\beta$ respectively, though many different choicwes are possible.

Appendix B Keldysh CCSD amplitude equations

The amplitude equations are determined by the kernels S1 and S2, which contain contractions identical to those in the CCSD T1 and T2 equations:

[TABLE]

[TABLE]

We have used $n_{p}$ to indicate the Fermi-Dirac occupation number of the $p$th orbital:

[TABLE]

Here and elsewhere, we use the notation:

[TABLE]

Appendix C Keldysh CCSD lambda equations

In practice, we solve for the discretized Lagrange multipliers, $\lambda(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

[TABLE]

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

[TABLE]

[TABLE]

The $\tilde{\lambda}$ are integrals of the $\lambda$ amplitudes with an exponential factor. We compute them separately along each branch of the contour:

[TABLE]

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.

Appendix D Keldysh CCSD 1-RDM

Given the $s_{\mu}$ and $\lambda_{\mu}$ 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:

[TABLE]

The value of some one-electron observable at time $t_{y}$ can be computed as

[TABLE]

We have used $\left<\ldots\right>_{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.

Appendix E 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.

E.1 Preliminaries

We will consider a 1-particle Hamiltonian of the form

[TABLE]

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

[TABLE]

The value of some observable at time $t$ is given by

[TABLE]

where

[TABLE]

as discussed in Section II.3. 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. White and Chan, 2018), 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 6.

The CCS and LCCS equations are shown diagrammatically in Figure 7.

A cursory analysis of these diagrams indicates that

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

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

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

E.2 Second order perturbation theory

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

[TABLE]

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

[TABLE]

where the amplitudes are now implicitly a function of $\alpha$ as well. The amplitudes are equal to

[TABLE]

which implies that the free energy as a function of $\alpha$ is

[TABLE]

Expanding $\Omega$ in $\alpha$ results in

[TABLE]

where $\Omega^{[2]}_{n}$ is the $n$th derivative wrt $\alpha$. We find two contributions to the first derivative:

[TABLE]

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

[TABLE]

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

[TABLE]

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. Santra and Schirmer, 2017), but independent of time ($t_{1}$).

E.3 Third order perturbation theory

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

[TABLE]

The corresponding contribution to the grand potential is

[TABLE]

[TABLE]

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 6(b). Each term has 3 factors of $f$ which means that there will be 6 total terms that are 1st order in $\alpha$:

[TABLE]

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

[TABLE]

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

[TABLE]

There is a slight subtlety here in that while we know that (1), (2), and (3) all have well-defined, finite, limits as $\Delta\rightarrow 0$, the same is not true of $(1)[t_{1}]$ for example. So, for the remainder of this section we will assume that $\Delta$ 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 $\Delta\rightarrow 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

[TABLE]

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

[TABLE]

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

[TABLE]

The results are that

[TABLE]

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

[TABLE]

Everything is now close to canceling completely. We note that

[TABLE]

which allows us to write

[TABLE]

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.

E.4 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 6(c):

[TABLE]

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 76 yields:

[TABLE]

Just as before, we expand this expression as a function of $\alpha$ to obtain

[TABLE]

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

[TABLE]

After a great deal of algebra, we obtain

[TABLE]

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 II.6, this also implies that $\Omega^{[4]}-\Omega^{[4]\ast}$ does not generate dynamics that conserve particle number. Since LCCS includes the $\Omega^{[4]}-\Omega^{[4]\ast}$ contribution at 4th order, this suggests that Keldysh-LCCS will not conserve particle number.

E.5 Numerical results

From this discussion and the discussion given in Section II.6, we conclude several things about particle number conservation in the 1-particle problem:

$n$th order perturbation theory should conserve particle number 2. 2.

For 2nd and 3rd order perturbation theory, each diagram individually conserves the particle number 3. 3.

But this is not true of 4th order perturbation theory where the ($\ast$) terms (and therefore the remaining terms also) do not conserve particle number 4. 4.

LCCS does not conserve particle number 5. 5.

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,

[TABLE]

For this example, we use

[TABLE]

[TABLE]

The relatively large perturbation is used to accentuate the defects in the theory. We use $\mu=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, $\left<N\right>$ is conserved.

In Table 1, 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 8-10 we confirm numerically point (1).

These figures indicate that for perturbation theory at any order, $\left<N\right>$ 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 11 we show that each of these contributions individually conserves $N$.

This confirms numerically the result of Section E.3.

Point (3) is confirmed numerically in Figure 12 where the $\ast$ contributions at 4th order are included and the particle number is no longer conserved.

Point(4) is confirmed in Figure 13 which qualitatively resembles Figure 12.

Finally, point (5) is corroborated in Figure 14.

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 1), 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 E.5 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.