Real-time density-matrix coupled-cluster approach for closed and open systems at finite temperature

TL;DR

This paper introduces a real-time density-matrix coupled-cluster method using thermo-field formalism for simulating quantum dynamics of closed and open systems at finite temperature, ensuring trace-preservation and computational simplicity.

Contribution

It extends coupled-cluster theory to finite-temperature quantum dynamics with an exponential ansatz and Keldysh contour, maintaining trace-preservation and simplifying calculations.

Findings

Successfully captures correlated dynamics of the Anderson model

Maintains exact trace-preservation during evolution

Demonstrates effectiveness for both closed and open systems

Abstract

We extend the coupled-cluster method to correlated quantum dynamics of both closed and open systems at finite temperatures using the thermo-field formalism. The approach expresses the time-dependent density matrix in an exponential ansatz and describes time-evolution along the Keldysh path contour. A distinct advantage of the approach is exact trace-preservation as a function of time, ensuring conservation of probability and particle number. Furthermore, the method avoids the computation of correlated bra-states, simplifying the computational implementation. We develop the method in a thermal quasi-particle representation, which allows seamless connection to the projection method and diagrammatic techniques of the traditional coupled-cluster formalism. For comparison, we also apply the thermo-field framework to renormalization-group methods to obtain reference results for closed and open systems at finite temperatures. We test the singles and doubles approximation to the density-matrix coupled-cluster method on the correlated electronic dynamics of the single-impurity Anderson model, demonstrating that the new method successfully captures the correlated dynamics of both closed systems at finite temperature and driven-dissipative open systems. This encouraging performance motivates future applications to non-equilibrium quantum many-body dynamics in realistic systems.

Equations

$${\left[a_i,a_j^{†}\right]}_{\sigma }={\delta }_{ij},{\left[a_i^{†},a_j^{†}\right]}_{\sigma }={\left[a_i,a_j\right]}_{\sigma }=0,$$

$$\rho =\sum _{n,m}{\rho }_{nm}|\mathbf{n}⟩⟨\mathbf{m}|$$

$$\begin{array}{rl}{\tilde{a}}_i|\mathbf{n}⟩⟨\mathbf{m}|& ={\sigma }^{\mu +1}|\mathbf{n}⟩⟨\mathbf{m}|a_i^{†},\\ {\tilde{a}}_i^{†}|\mathbf{n}⟩⟨\mathbf{m}|& ={\sigma }^{\mu }|\mathbf{n}⟩⟨\mathbf{m}|a_i\end{array}$$

$$\begin{array}{rl}{\left[{\tilde{a}}_i,{\tilde{a}}_j^{†}\right]}_{\sigma }={\delta }_{ij},{\left[{\tilde{a}}_i^{†},{\tilde{a}}_j^{†}\right]}_{\sigma }& ={\left[{\tilde{a}}_i,{\tilde{a}}_j\right]}_{\sigma }=0,\\ {\left[{\tilde{a}}_i^{†},a_j^{†}\right]}_{\sigma }& ={\left[{\tilde{a}}_i,a_j\right]}_{\sigma }=0,\end{array}$$

$$\begin{array}{rl}{\left(AB\right)}^{\stackrel{~}{}}& =\tilde{A}\tilde{B},\\ {\left(c_{a}A+c_{b}B\right)}^{\stackrel{~}{}}& =c_a^{\ast }\tilde{A}+c_b^{\ast }\tilde{B},\end{array}$$

$$\rho =\sum _{n,m}{\rho }_{nm}|\mathbf{n},\tilde{\mathbf{m}}⟩.$$

$$|{\rho }_0⟩=\frac{1}{Q_0}\sum _n{e}^{-\beta {\sum }_i{n}_i{ϵ}_i}|\mathbf{n},\tilde{\mathbf{n}}⟩=\frac{1}{Q_0}\text{exp}\left(\sum _i{e}^{-\beta {ϵ}_i}{a}_i^{†}{\tilde{a}}_i^{†}\right)|0⟩,$$

$$|1⟩=\sum _n|\mathbf{n},\tilde{\mathbf{n}}⟩=\text{exp}\left(\sum _i{a}_i^{†}{\tilde{a}}_i^{†}\right)|0⟩.$$

$$\text{tr}\left[\rho O\right]=⟨1|O|\rho ⟩,$$

$$i|{\dot{\rho }}_t⟩=\hat{H}|{\rho }_t⟩,$$

$$\hat{H}=H-\tilde{H}.$$

$$|{\rho }_t⟩=e^{-i\hat{H}t}|{\rho }_{t_0}⟩$$

$$\begin{array}{rl}⟨1|a_i^{†}& =⟨1|{\tilde{a}}_i,\\ ⟨1|{\tilde{a}}_i^{†}& =\sigma ⟨1|a_i.\end{array}$$

$$⟨1|\hat{H}=0.$$

$$⟨1|O|{\rho }_t⟩=⟨1|O_t|{\rho }_{t_0}⟩$$

$$\begin{array}{rl}{a}_i|{\rho }_0⟩& =e^{-\beta {ϵ}_i}{\tilde{a}}_i^{†}|{\rho }_0⟩,\\ {\tilde{a}}_i|{\rho }_0⟩& =\sigma e^{-\beta {ϵ}_i}{a}_i^{†}|{\rho }_0⟩.\end{array}$$

$$\begin{array}{rl}⟨1|b_i^{†}& =⟨1|{\tilde{b}}_i^{†}=0,\\ b_i|{\rho }_0⟩& ={\tilde{b}}_i|{\rho }_0⟩=0.\end{array}$$

$$\begin{array}{rl}{b}_i^{†}=a_i^{†}-{\tilde{a}}_i,& {\tilde{b}}_i^{†}={\tilde{a}}_i^{†}-\sigma a_i,\\ b_i=u_i{a}_i-v_i{\tilde{a}}_i^{†},& {\tilde{b}}_i=u_i{\tilde{a}}_i-\sigma v_i{a}_i^{†},\end{array}$$

$$|{\rho }_t⟩=e^{T_t}|{\rho }_{t_0}⟩,$$

$$T_t=\sum _{n=1}^N{T}_n,$$

$$T_n=\frac{1}{n!}\sum _{ij...,kl...}{t}_{jl...}^{ik...}{b}_i^{†}{b}_k^{†}...{\tilde{b}}_l^{†}{\tilde{b}}_j^{†},$$

$$T_t=\sum _{ij}{t}_j^i{b}_i^{†}{\tilde{b}}_j^{†}+\frac{1}{4}\sum _{ijkl}{t}_{jl}^{ik}{b}_i^{†}{b}_k^{†}{\tilde{b}}_l^{†}{\tilde{b}}_j^{†}.$$

$$G=\frac{\text{tr}\left[e^{-\beta H}O(t_N)..O(t_1)\right]}{\text{tr}\left[e^{-\beta H}\right]}=\frac{⟨T_c{U}_c{O}_N..O_1⟩}{⟨T_c{U}_c⟩}.$$

$$G=\frac{⟨1|T{e}^{-i\hat{H}t}{O}_N..O_1|{\rho }_0⟩}{⟨1|e^{-i\hat{H}t}|{\rho }_0⟩}=⟨1|O(t_N)..O(t_1)|{\rho }_0⟩,$$

$$i{e}^{-T_t}{\partial }_t{e}^{T_t}|{\rho }_0⟩=e^{-T_t}\hat{H}{e}^{T_t}|{\rho }_0⟩,$$

$$e^{-T_t}{\partial }_t{e}^{T_t}={\partial }_t+{\dot{T}}_t,$$

$$i{\dot{T}}_t|{\rho }_0⟩={\left(\hat{H}\sum _{n=0}^{n=N}\frac{1}{n!}{T}_t^n\right)}_c|{\rho }_0⟩,$$

$$i⟨0|{\tilde{b}}_j..b_i{\dot{T}}_t|{\rho }_0⟩=⟨0|{\tilde{b}}_j..b_i{\left(\hat{H}\sum _{n=0}^{n=N}\frac{1}{n!}{T}_t^n\right)}_c|{\rho }_0⟩.$$

$$H_0=\sum _{i=0s}^{N_b}{ϵ}_i{a}_{is}^{†}{a}_{is}+\sum _{i=1,s}^{N_b}\left(V_i{a}_{0s}^{†}{a}_{is}+a_{is}^{†}{a}_{0s}\right),$$

$$W=U{a}_{0\alpha }^{†}{a}_{0\alpha }{a}_{0\beta }^{†}{a}_{0\beta }$$

Full text

I Introduction

The simulation of non-equilibrium quantum many-body dynamics is an outstanding challenge in chemical physics. The theoretical description of non-equilibrium phenomena is essential for the understanding of many processes, including chemical dynamics of molecules,Krausz and Ivanov (2009); Domcke and Yarkony (2012) defects in extended systems,Doherty et al. (2013); Falk et al. (2015); Freysoldt et al. (2014) quantum transport in molecular junctions,Thoss and Evers (2018); Cui et al. (2017) and chemical reactions at surfaces and interfaces,Zhang et al. (2017); Jiang and Guo (2019); Maurer et al. (2016) as well as the description of novel driven-dissipative phenomena in open systems, such as discrete time crystals.Choi et al. (2017); Else et al. (2019) These open questions motivate the extension of well-established wavefunction approaches of quantum chemistry to correlated quantum dynamics of closed and open systems at finite temperature.

Among the methods of quantum chemistry, the coupled-cluster approach has become the golden standard for incorporation of correlation effects in molecular systems.Bartlett and Musiał (2007); Crawford and Schaefer (2000) The approach relies on an exponential ansatz for the wavefunction of the system, which is parameterized by connected cluster amplitudes. The singles and doubles approximation to the cluster operator together with a perturbative triples correction has set the accuracy standard for ground state properties of molecules, and excited state properties have been targeted via linear response theory.Bartlett (2012); Nooijen and Bartlett (1997); Krylov (2008) Furthermore, time-dependent coupled-cluster wavefunction methods have been recently applied to the non-equilibrium electronic dynamics of molecules,Huber and Klamroth (2011); Sato et al. (2018); Kvaal (2012) and beyond the single-reference framework, the coupled-cluster approach has been extended to more challenging, multi-reference systems.Lyakh et al. (2011)

Thermal effects have been recently included in the coupled-cluster formalismWhite and Chan (2018); Hummel (2018) using the method of thermal cumulants.Sanyal, Mandal, and Mukherjee (1992); Sanyal et al. (1993); Mandal et al. (2003) Alternatively, thermo-field theory has been combined with the method of configuration interaction to simulate finite-temperature molecular properties.Harsha, Henderson, and Scuseria (2019) The thermo-field methodTakahashi and Umezawa (1996); Chu and Umezawa (1994) was developed to incorporate thermal mixed states in quantum field theory and the approach was later extended to open systems.Arimitsu and Umezawa (1985, 1987a); Arimitsu, Guida, and Umezawa (1988); Arimitsu and Umezawa (1987b) Besides configuration interaction treatment, the method has allowed direct application of coupled-cluster theory to non-equilibrium stationary states,Dzhioev and Kosov (2011a, 2012, b, 2014a, 2014b) and it has been also combined with renormalization-group methodsde Vega and Bañuls (2015); Schwarz et al. (2018); Borrelli (2019); Borrelli and Gelin (2016, 2017) to simulate the real-time finite-temperature dynamics of closed systems, including impurity-bath models and coupled electron-nuclear dynamics.

Our goal is to extend the coupled-cluster approach to the real-time correlated quantum dynamics of closed and open systems at finite temperatures, bridging the gap between recent finite-temperature and real-time coupled-cluster extensions. We provide a unified framework to describe both equilibrium and non-equilibrium stationary states and non-equilibrium quantum dynamics. We use the thermo-field theory to represent the mixed quantum states of closed and open systems in a wavefunction form, and we express this wavefunction in an exponential ansatz. We test the singles and doubles approximation to the new density-matrix coupled-cluster method on the correlated electronic dynamics of the single-impurity Anderson modelNewns (1969); Anderson (1961); Wilson (1975) and compare the results with the renormalization-group method.

The paper is organized as follows. We summarize thermo-field theory and present the real-time finite-temperature density-matrix coupled-cluster method for closed and open systems in Sec. II. We present the results for the quantum dynamics of the single-impurity Anderson model at a range of interaction strengths and temperatures in Sec. III, including also dynamics with driving and dissipation. We conclude in Sec. IV with remarks on the future directions for method development.

II Theory

We start this section by reviewing the thermo-field formalism and the thermal Bogoliubov transformation, which we use to define temperature-dependent quasi-particles. We then present the main result of this paper, the extension of the real-time coupled-cluster method to closed and open systems at finite temperatures using the thermal quasi-particle representation. We focus next on the test system for our quantum dynamics approximations, the single-impurity Anderson model, and the extension of renormalization group approaches to closed and open systems, which we compare with the coupled-cluster approximation.

II.1 Thermo-field theory

Thermo-field theoryTakahashi and Umezawa (1996); Arimitsu and Umezawa (1985, 1987b); Schmutz (1978) expresses the mixed states of quantum systems in the form of pure state wavefunctions. Let us define a basis in the Hilbert space $\mathcal{H}$ of a system by the ket-states $\ket{\bm{n}}$, generated via the action of the creation operators $a^{\dagger}_{i}$ on the vacuum ket $\ket{0}$, $\ket{\bm{n}}=\prod_{i}a^{\dagger}_{i}\ket{0}$, where $\bm{n}$ is the set of occupation numbers and the index $i$ enumerates the states of the single-particle basis. The creation operators $a^{\dagger}_{i}$ and their Hermitian conjugate annihilation operators $a_{i}$ obey the commutation relations

[TABLE]

where the commutator is $\left[A,B\right]_{\sigma}=AB-\sigma BA$ with $\sigma=1$ for bosons and $\sigma=-1$ for fermions. Using the basis $\ket{\bm{n}}$ of $\mathcal{H}$, the density matrix $\rho$ is given by

[TABLE]

with the ket-bra couples $\ket{\bm{n}}\bra{\bm{m}}$ forming an orthonormal basis. To write Eq. (2) in the form of a wavefunction, thermo-field theory uses a second replica $\tilde{\mathcal{H}}$ of the original Hilbert space $\mathcal{H}$ to express the ket-bra couples $\ket{\bm{n}}\bra{\bm{m}}$ as ket-basis states $\ket{\bm{n},\tilde{\bm{m}}}$, effectively doubling the number of degrees of freedom of the system. We define tilde creation $\tilde{a}^{\dagger}_{i}$ and annihilation $\tilde{a}_{i}$ operators in the tilde Hilbert space $\tilde{\mathcal{H}}$ via the relations

[TABLE]

with $\mu=\sum_{i}(n_{i}-m_{i})$, the difference of ket- and bra-state occupation numbers. From Eqs. (1) and  (3), the non-tilde and tilde operators satisfy the commutation relations

[TABLE]

which show that the non-tilde and tilde operators represent independent fermionic or bosonic degrees of freedom. Given that the tilde Hilbert space is Hermitian conjugate to the original Hilbert space, thermo-field theory relates operators acting on the two subspaces via the anti-linear tilde conjugation operation with the properties

[TABLE]

with $c_{a}$ and $c_{b}$ complex numbers and $*$ denoting complex conjugation. The double tilde conjugation rules $\tilde{\tilde{a}}_{i}=\sigma a_{i}$ and $\tilde{\tilde{a}}^{\dagger}_{i}=\sigma a^{\dagger}_{i}$ supplement Eq. (5).

In the doubled Hilbert space $\mathcal{H}\crossproduct\tilde{\mathcal{H}}$ of non-tilde (physical) and tilde (hidden) degrees of freedom, the density matrix, Eq. (2), takes the form of a pure state wavefunction

[TABLE]

The thermal density matrix $\ket{\rho_{0}}$ of a non-interacting system, for instance, can be written in the exponential form

[TABLE]

where $\beta$ is the reciprocal temperature, $\beta=(k_{B}T)^{-1}$ with $T$ the temperature, $Q_{0}$ is the partition function, $\epsilon_{i}$ are the single-particle energy levels, and $\ket{0}=\ket{0,\tilde{0}}$ is the vacuum ket-state in the doubled Hilbert space. Similarly, the resolution of identity, the unit ket-state $\ket{1}$, also has an exponential form

[TABLE]

The average value of an operator $O$, given by the trace with the density matrix $\rho$, becomes an expectation value with the ket-state $\ket{\rho}$, Eq. (6), in thermo-field theory. Specifically,

[TABLE]

where the density matrix is normalized $\bra{1}\ket{\rho}=1$. The bra-state and ket-state in Eq. (9) are different due to the choice of the path contour in the complex time plane as shown in Sec. II.3. This choice is most suitable for extending the coupled-cluster approach to finite temperatures because it circumvents the need to compute a correlated bra-state. We provide further details in Sec. II.3 devoted to the coupled-cluster formalism.

The Liouville-von Neumann equation for the time-evolving density matrix, $i\dot{\rho}=\left[H,\rho\right]$, becomes Schrödinger equation for the pure ket-state $\ket{\rho}$

[TABLE]

where the super-Hamiltonian operator $\hat{H}$ is the difference of the original Hamiltonian $H$ and its tilde conjugate

[TABLE]

The super-Hamiltonian in Eq. (11) does not mix the non-tilde and tilde degrees of freedom for a closed system, resulting in a unitary evolution. For open, dissipative systems, the super-Hamiltonian entangles the non-tilde and tilde degrees of freedom by terms that resemble pairing of opposite spin electrons in the theory of superconductivity. We will give an example of such a super-Hamiltonian in Sec. II.4.

Formal integration of Eq. (10) gives the time-dependent density matrix $\ket{\rho_{t}}$ as

[TABLE]

with $\ket{\rho_{t_{0}}}$ the initial density matrix. In quantum-quench simulations, we initiate the quantum dynamics from an uncorrelated state and take into account the effects of interactions by the subsequent time evolution. In Sec. III, we present results of such simulations where most of the system, the bath, is initially at thermal equilibrium with a density matrix given by Eq. (7). Non-equilibrium states relax to stationary states, such as equilibrium states, with time-independent density matrices $\ket{\rho_{s}}$, which are annihilated by the super-Hamiltonian $\hat{H}\ket{\rho_{s}}=0$.

II.2 Thermal quasi-particles

The creation $a^{\dagger}_{i}$ and $\tilde{a}^{\dagger}_{i}$ and annihilation $a_{i}$ and $\tilde{a}_{i}$ operators have a complementary action on the unit bra-state

[TABLE]

Because the super-Hamiltonian, Eq. (11), conserves the number of particles and contains equal numbers of non-tilde and tilde operators, Eq. (13) gives

[TABLE]

This bra-state condition ensures the conservation of probability and the normalization of the time-dependent density matrix, Eq. (12), at all times, and it allows us to express the expectation value, Eq. (9), in the Heisenberg picture

[TABLE]

with $O_{t}=e^{i\hat{H}t}Oe^{-i\hat{H}t}$ – the Heisenberg-evolved operator $O$.

The action of the $a_{i}$ and $\tilde{a}_{i}$ operators on the non-interacting thermal density matrix is similarly complementary

[TABLE]

Because $a_{i}$ annihilates an electron, the thermal ket-state condition, Eq. (16), shows that $\tilde{a}^{\dagger}_{i}$ creates a hole, and vice versa because $a^{\dagger}_{i}$ creates an electron, $\tilde{a}_{i}$ annihilates a hole.

Eqs. (13) and  (16) imply that we can transform the creation $a^{\dagger}_{i}$ and $\tilde{a}^{\dagger}_{i}$ and annihilation $a_{i}$ and $\tilde{a}_{i}$ operators to a new set of creation $b^{\dagger}_{i}$ and $\tilde{b}^{\dagger}_{i}$ and annihilation $b_{i}$ and $\tilde{b}_{i}$ operators, which respectively annihilate the unit bra-state and non-interacting thermal ket-state

[TABLE]

The transformation that satisfies Eq. (17) and preserves the commutation relations, Eqs. (1) and  (4), is given by

[TABLE]

where $v_{i}=\left(e^{\beta\epsilon_{i}}-\sigma\right)^{-1}$ is the Bose-Einstein or the Fermi-Dirac distribution and $u_{i}=1+\sigma v_{i}$. Eq. (18) preserves the commutation relations among the new set of creation and annihilation operators but does not preserve the Hermitian conjugation relation between creation and annihilation operators. Such canonical, non-unitary transformations are known in the theory of superconductivity as Bogoliubov transformations.Blaizot and Ripka (1986) The thermal Bogoliubov transformation, Eq. (18), defines new thermal quasi-particles $b^{\dagger}_{i}$ and $b_{i}$ and thermal quasi-holes $\tilde{b}^{\dagger}_{i}$ and $\tilde{b}_{i}$ with respect to the unit bra-state and non-interacting thermal ket-state, Eq. (17), which we refer to as the new left $\bra{0}$ and right $\ket{0}$ vacuum states. We use the thermal quasi-particle representation to introduce correlated dynamics via coupled-cluster theory at finite temperatures in closed and open systems in the following section.

II.3 Density-matrix coupled-cluster approach

In our coupled-cluster approach, we express the time-dependent density matrix of the system in the exponential form

[TABLE]

where $\ket{\rho_{t_{0}}}=\ket{\rho_{0}}$ is the non-interacting thermal density matrix in Eq. (7), and $T_{t}$ is the time-dependent cluster operator, which we expand in the number of excited thermal quasi-particle-hole pairs $n$ as

[TABLE]

where $N$ is the highest excitation level. We express the cluster operators $T_{n}$ via the quasi-particle $b^{\dagger}_{i}$ and quasi-hole $\tilde{b}^{\dagger}_{i}$ creation operators from Sec. II.2 as

[TABLE]

where $t^{ik...}_{jl...}$ are the time-dependent cluster amplitudes, and the indeces $i$, $j$, $k$, and $l$ run over the entire single-particle basis. We test the performance of the density-matrix coupled-cluster (DMCC) method keeping terms with up to two quasi-particle-hole pairs in Eq. (20). The explicit form of the cluster operator for this singles and doubles approximation is

[TABLE]

Eq. (19) together with Eqs. (20) and  (21) parameterize a class of trace-preserving transformations that include the subset of dynamical completely positive transformations of open systems and the narrower subset of unitary transformations of closed systems. The trace-preservation directly follows from the definition of the $b^{\dagger}$ and $\tilde{b}^{\dagger}_{i}$ creation operators, Eq. (17), $\bra{1}\ket{\rho_{t}}=\bra{1}e^{T_{t}}\ket{\rho_{0}}=\bra{1}\ket{\rho_{0}}$. The time evolution governed by the Schrödinger equation, Eq. (10), imposes further constraints on the density matrix in addition to trace-preservation, such as the Hermiticity and positivity of the density matrix. These constraints translate to relations among the cluster amplitudes. The Hermiticity of the density matrix, for instance, gives invariance upon tilde conjugation in thermo-field theory, $\ket{\tilde{\rho_{t}}}=\ket{\rho_{t}}$, resulting in a tilde invariant cluster operator $\tilde{T_{t}}=T_{t}$. For a singles-only approximation, this leads to the Hermiticity of the singles amplitude matrix $t^{i}_{j}=t^{j*}_{i}$ and similar relations hold for the higher rank cluster amplitudes. The dynamical constraints determined by the Schrödinger equation need not be incorporated in the parameterization of the cluster amplitudes; rather, their satisfaction may be used as diagnostics for the solutions to the Schrödinger equation.

The trace-preservation property of the coupled-cluster density matrix allows us to express quantum-statistical averages of time-dependent operators in a simple way.Matsumoto, Nakano, and Umezawa (1984) Thermal averages of Heisenberg-evolved operators $O(t)$ are given in general by the ratio

[TABLE]

In diagrammatic terms, the denominator in Eq. (23) cancels the disconnected vacuum diagrams in the numerator, such that only connected diagrams contribute to the thermal average. In Eq. (23), we have also introduced a more succinct expression for the average that involves a time-ordered propagation with $T_{c}$ the chronological operator and $U_{c}$ the full quantum propagator along a contour in the complex time plane, which initiates at time $t_{i}$ before time $t_{1}$, evolves along the real time axis to $t_{f}$ past the time $t_{N}$, comes back along the real time axis to $t_{i}$, and propagates downwards along the imaginary time axis to time $t_{i}-i\beta$. We can take respectively the initial $t_{i}$ and final $t_{f}$ times to the infinite past, $t_{i}\rightarrow-\infty$, and infinite future, $t_{f}\rightarrow\infty$, where the interaction is turned off, resulting in cancellation of the correlated contribution to the imaginary axis propagation. Such a contour in the complex time plane underlies non-equilibrium Keldysh diagrammatic theoryKeldysh et al. (1965); Altland and Simons (2010) and the density-matrix coupled-cluster theory presented here.

Using the definition of expectation value Eq. (9) and the super-Hamiltonian Eq. (11), we cast Eq. (23) in thermo-field theory as

[TABLE]

where the imaginary axis propagation is represented by the non-interacting thermal state $\ket{\rho_{0}}$ and the contour runs only along the real time axis from $t_{i}$ to $t_{f}$. The second equality in Eq. (24) follows from the trace-preservation of the super-Hamiltonian, allowing us to express the thermal average in an explicitly connected form. In our density-matrix coupled-cluster method, we rigorously impose the trace-preservation of the time-evolved density matrix, resulting in explicit cancellation of the disconnected contributions to the thermal average, which avoids the construction of a correlated bra-state.

Other choices for the path contour in the complex time plane are also possible because of the analytic property of thermal averages. These choices lead to equivalent formulations of thermal averages in thermo-field theory but may not lead to equivalent coupled-cluster approximations. One other choice, for instance, is to place the backward portion of the real-time propagation at imaginary time $-i\beta/2$ rather than immediately below the real time axis.Matsumoto, Nakano, and Umezawa (1984) This contour gives unitary formulation of thermal averages that has bra-states Hermitian conjugate to the ket-states. If we express the time-dependent density matrix in this unitary formulation in an exponential form, then we also need to compute the correlated bra-state to evaluate the thermal average, unless we use a unitary form for the cluster operator.

We derive the equation of motion for the cluster amplitudes by using Eq. (19) in the Schrödinger equation Eq. (10) for the density matrix ket-state.Crawford and Schaefer (2000) Following multiplication from the left with $e^{-T_{t}}$, we obtain

[TABLE]

where $\hat{H}$ is the super-Hamiltonian of the system expressed in the thermal quasi-particle representation. Because $\partial_{t}$ is a single-particle operator and the creation quasi-particle operators are time-independent, the Baker-Campbell-Hausdorff identity gives

[TABLE]

where the dot stands for the time derivative. Because the super-Hamiltonian has a finite number of creation and annihilation quasi-particle operators and the cluster operators, Eqs. (20) and  (21), commute, application of the Baker-Campbell-Hausdorff identity shows that the expansion of the similarity-transformed super-Hamiltonian $e^{-T_{t}}\hat{H}e^{T_{t}}$ contains a finite number of nested commutator terms. For the standard electronic structure super-Hamiltonian, which has terms with at most four annihilation quasi-particle operators, the commutator series exactly truncates at fourth order. Thus, the equation of motion for the cluster operators, Eq. (25), simplifies to

[TABLE]

where $N$ stands for the largest order non-vanishing term in the series, and we used the time-independence of the ket-state $\partial_{t}\ket{\rho_{0}}=0$. The subscript "c" in Eq. (27) denotes that only the connected part of the expression in the brackets has to be considered as a result of the cancellation of the disconnected part by the nested commutators. To derive differential equations for the cluster amplitudes, we project Eq. (27) onto states with finite number of quasi-particle excitations, $\bra{0}\tilde{b}_{j}..b_{i}$,

[TABLE]

In our singles and doubles approximation to the DMCC method, Eq. (22), we project onto bra-states that have only single and double quasi-particle excitations. These equations contain vacuum expectation values of strings of creation and annihilation operators, which when evaluated give a set of coupled non-linear ordinary differential equations for the cluster amplitudes that we solve using the Runge-Kutta method.Press et al. (2007)

A distinct advantage of the thermal quasi-particle representation is the existence of traditional normal ordering of operators and the applicability of Wick’s theorem for evaluation of expectation values of strings of creation and annihilation operators.Blaizot and Ripka (1986) We define normal-ordering as the string of operators, which has all creation operators to the left of all annihilation operators. Contractions of quasi-particle operators, the difference of any operator ordering from the normal order, are straightforward to evaluate because the quasi-particle operators satisfy the usual commutation relations. Given that the creation operators $b^{\dagger}_{i}$ and $\tilde{b}^{\dagger}_{i}$ annihilate the left vacuum $\bra{1}$, vacuum expectation values of normal-ordered strings of operators are identically zero. Using Wick’s theorem, we can thus show that the vacuum expectation value of any string of creation and annihilation quasi-particle operators is equal to the sum of all possible pairwise contractions. If sub-strings are already in normal-ordered form, then contractions need to be considered only among the normal-ordered sub-strings. The number of contractions, however, grows quickly with the complexity of the operator products and diagrammatic representations, which express the contributions to the coupled-cluster equations in a graphical form, become useful. The same set of rules for drawing zero-temperature time-independent coupled-cluster diagramsBartlett and Musiał (2007); Crawford and Schaefer (2000) applies to the density-matrix coupled-cluster method with upward going lines representing non-tilde particles and downward going lines – tilde particles according to the interpretation of non-tilde particles as thermal quasi-particles and tilde particles as thermal quasi-holes. Summation over "internal" lines runs over the entire single-particle basis, and the sign of diagrams is computed by $(-1)^{l}$, where $l$ is the number of loops determined in the usual way. The direct applicability of the diagrammatic techniques provides powerful methods for evaluating real-time finite-temperature coupled-cluster approximations.

II.4 Single-impurity Anderson model

We apply the DMCC method with single and double quasi-particle excitations to the single-impurity Anderson model (SIAM), which is a well-established system to understand the role of electron correlation in dynamical processes at surfaces. SIAM consists of a localized impurity with a repulsive interaction between opposite spin electrons that is coupled to a non-interacting electronic bath.Newns (1969); Anderson (1961); Wilson (1975) The SIAM Hamiltonian $H=H_{0}+W$ has a non-interacting part $H_{0}$ given by

[TABLE]

where $\epsilon_{i}$ are the single-particle energy levels, $V_{i}$ – the impurity-bath coupling elements, $i=0,...,N_{b}$ enumerates the single-particle basis where $i=0$ corresponds to the impurity and $N_{b}$ is the number of bath levels, and $s$ sums over the spins of the electrons. The interaction $W$ is given by

[TABLE]

with $U$ the Hubbard parameter, a measure of the interaction strength. Eqs. (29) and  (30) represent a closed system with a super-Hamiltonian that is the difference of $H$ and its tilde conjugate, $\hat{H}=H-\tilde{H}=H_{0}-\tilde{H}_{0}+W-\tilde{W}$. Dissipative coupling of the bath levels to a continuum can be incorporated in the SIAM super-Hamiltonian via terms that mix non-tilde and tilde operators,Arimitsu and Umezawa (1985) such as

[TABLE]

with $\gamma^{(1)}_{i}=\gamma(1-v_{i})$, $\gamma^{(2)}_{i}=\gamma v_{i}$, and $\gamma$ the dissipation strength. The relation among $\gamma^{(1)}_{i}$, $\gamma^{(2)}_{i}$ and $\gamma$ follows from the detailed balance condition, and $\hat{D}$ conforms with the dissipative part of a Lindblad operatorLindblad (1976) that represents open systems with Markovian dynamics. Non-Markovian effects of the impurity dynamics are captured by explicitly incorporating a finite electronic bath. Together with dissipation, we include driving in SIAM via a time-dependent modulation of the electronic level of the impurity, $H_{t}=\epsilon_{0}+\delta\epsilon\text{sin}(\Omega t)$. We use logarithmic discretization of the electronic bath in our calculations as suggested by Wilson for the Kondo modelWilson (1975) with a discretization parameter $\Lambda=1.1$ and $N_{b}=100$. We consider also the more challenging asymmetric case of SIAM with $\epsilon_{0}=-0.08$ eV, impurity-bath coupling strength $V=0.04$ eV, and a range of interaction strengths $U$ and temperatures $T$.

Using the thermal Bogoliubov transformation, Eq. (18), the SIAM super-Hamiltonian $\hat{H}=\hat{H}_{0}+\hat{W}+\hat{D}$ becomes $\hat{H^{\prime}}=\hat{H^{\prime}_{0}}+\hat{W^{\prime}}$, where the non-interacting super-Hamiltonian $\hat{H^{\prime}_{0}}$ is given by

[TABLE]

with the single-particle Hamiltonian matrix $h_{ijs}$,

[TABLE]

and the pairing matrix $\Delta_{ijs}$,

[TABLE]

and the interaction $\hat{W^{\prime}}$ is written as

[TABLE]

with t.c. standing for tilde conjugation. In Eq. (33), $\epsilon^{\prime}_{i}=\epsilon_{i}-i\gamma$ are the broadened levels of the electronic bath due to dissipation, and the single-particle energy of the impurity is modified by a thermal mean-field term $v_{0}U$. The Hubbard interaction in Eq. (35) also incorporates thermal effects as the interaction strength depends on the thermal population of the impurity. Because $\hat{W^{\prime}}$ includes contributions with up to three annihilation operators, the series in the equation of motion for the cluster operator, Eq. (27), truncates at third order, simplifying the differential equations for the cluster amplitudes. In our quantum quench simulations, the impurity is initially decoupled from the bath, $V=0$. Thus, the thermal quasi-particles are defined with respect to the non-interacting thermal state of the decoupled system. At time $t=0$, the coupling to the electronic bath is turned on, and we compute the non-equilibrium electronic dynamics of the impurity.

II.5 Finite-temperature renormalization-group approach for open systems

The quasi-particle representation of thermo-field theory is a suitable framework to extend time-dependent renormalization-group methodsVidal (2003, 2004); White and Feiguin (2004); Schollwöck (2011) to closed and open systems at finite temperature. We apply the time-evolving block-decimation approachVidal (2003, 2004) to SIAM and compare the results to DMCC in Sec. III. We represent the time-evolving density-matrix ket-state as a matrix product state, for which we arrange the electronic sites in a linear array, starting with the $\alpha$ impurity site, followed by $N_{b}$ $\alpha$ bath sites, the $\beta$ impurity site, and $N_{b}$ $\beta$ bath sites. Each electronic site carries both non-tilde and tilde degrees of freedom. The terms of the SIAM super-Hamiltonian $\hat{H^{\prime}}$ act on two sites at a time with the interaction $\hat{W^{\prime}}$, Eq. (35), coupling the impurity sites and the non-interacting super-Hamiltonian $\hat{H^{\prime}_{0}}$, Eq. (32), acting between the impurity site and each of the bath sites. The super-Hamiltonian terms, however, are non-local as they couple distant sites of the linear array, and we use fermionic swap operators $\hat{S}$ to rearrange sites given by

[TABLE]

where

[TABLE]

We apply the second-order Suzuki-Trotter decompositionSuzuki (1985) to the super-Hamiltonian terms in the time-evolution operator, which together with the swap operators results in the following sequence of two-site gates for a half time-step $\delta t$

[TABLE]

where $i$ runs over the sequence of $\alpha$ sites and $j$ over the sequence of $\beta$ sites, respectively. The full evolution operator is the product of Eq. (38) and its transpose. The application of each of the two-site gates in Eq. (38) requires the singular value decomposition of the resultant tensor to obtain the time-evolved matrix product state. Only singular values above a certain threshold are kept in the matrix product state, which represents the decimation step of the renormalization-group algorithm. At zero temperature, the doubled set of tilde and non-tilde degrees of freedom is not necessary because we can evolve only the wavefunction of the system. We apply the same Suzuki-Trotter decomposition using the SIAM Hamiltonian, instead of the super-Hamiltonian, in this case.

III Results

We present the DMCC results for the electronic dynamics of SIAM and compare them with the renormalization-group approach, which we regard as numerically exact. We carry out quantum quench simulations where the impurity is initially decoupled from the bath and is occupied by a single electron, a state with unit total population and unit spin polarization. We turn the coupling to the bath on at $t=0$ and monitor the electronic population dynamics of the impurity.

Figure  1 shows the time evolution of the total population (A) and spin polarization (B) of the impurity for three values of the interaction strength $U$ at zero temperature $T=0$. The DMCC singles (DMCC-S) approximation is exact for the non-interacting SIAM $U=0$ as demonstrated by the correspondence to the reference data. With increasing interaction strength, the DMCC-S approximation overestimates the total population and spin polarization of the impurity with larger deviations at higher Hubbard parameter values, revealing the importance of correlation effects. Inclusion of double quasi-particle excitations via the DMCC singles and doubles (DMCC-SD) approximation improves the electronic population dynamics of the impurity. Both the total population and spin polarization of the impurity decrease as a result of correlated dynamics, which brings the coupled-cluster performance in close correspondence to the exact renormalization-group results.

Figure 2 depicts the time evolution of the total population (A) and spin polarization (B) of the interacting impurity ($U=0.1$ eV) at three values of temperature $T$. Whereas the zero-temperature SIAM can also be simulated via the time-dependent coupled-cluster wavefunction method, the finite temperature cases require a density-matrix treatment. Both the total population and spin polarization of the impurity decrease with increasing temperature. Similarly, the correlation contribution to the electronic dynamics of the impurity decrease with temperature as measured by the DMCC-S deviation from the reference results. Thus, the DMCC-S approximation, which is equivalent to the time-dependent finite-temperature mean-field method, accurately describes the population dynamics of the impurity at high temperatures. Correlation effects, however, are significant in the low and intermediate temperature regimes and the close correspondence of DMCC-SD with the exact results demonstrates that the method successfully captures the correlation contribution at finite temperatures, as well.

The models considered so far represent closed quantum systems with dynamics characterized by relaxation to the equilibrium state of the interacting system. The DMCC method also allows us to investigate the electronic dynamics of driven-dissipative open quantum systems. In particular, we carry out quantum quench simulations of the driven-dissipative SIAM with a super-Hamiltonian that includes a dissipative contribution, Eq. (31), and a harmonic driving contribution, $H_{t}=\epsilon_{0}+\delta\epsilon\text{sin}(\Omega t)$ with $\delta\epsilon=0.08$ eV. Figure 3 shows the total population dynamics of the interacting impurity ($U=0.1$ eV) at temperature $T=0.04$ eV and two values of the modulation frequency, $\Omega=\pi V$ and $\Omega=4\pi V$. Figure 3A also includes dissipation with strength $\gamma=0.04$ eV, and Figure 3B – dissipation with strength $\gamma=0.2$ eV. Figure 3A reveals that the amplitude of impurity population oscillations decreases with increasing modulation frequency. This non-adiabatic effect results from the inefficient equilibration of the impurity with the electronic bath, and is reproduced by both DMCC approximations. The DMCC-S results, however, show asymmetric deviations from the reference data, especially in the low-frequency case where the DMCC-S deviations are larger near the peaks of the population oscillations. This time-evolving correlation effect is captured by DMCC-SD and results from the decrease of the correlation contribution as the impurity depopulates. Comparison of the results in Figures 3A and 3B further shows that both the amplitude of oscillations and average value of impurity population decrease with increasing dissipation strength. The coincidence of DMCC-S, DMCC-SD and exact results demonstrates that dissipation depresses correlation effects, similarly to increasing temperature. Overall, the DMCC-SD approximation reproduces well the rich non-equilibrium behavior of the model, underscoring the promise of the approach to study systems far from equilibrium.

Finally, Figure 4 demonstrates the particle number conservation during the coupled-cluster time evolution. The number of electrons is conserved within numerical convergence thresholds for the closed system at zero and finite temperature, and shows characteristic oscillations for the driven-dissipative, open quantum system. Taken together, the DMCC-SD method successfully captures the non-equilibrium many-body dynamics of SIAM in close correspondence with the renormalization-group method.

IV Conclusions

The current work presents an extension of coupled-cluster theory to real-time finite-temperature correlated quantum dynamics. The method uses the Keldysh path contour in the complex time plane, which incorporates thermal correlation effects via thermalization from a non-interacting thermal state. The density matrix is expressed in an exponential ansatz and the equation of motion for the cluster amplitudes is obtained via projection from the Schrödinger equation for the cluster operator. The thermal quasi-particle presentation of thermo-field theory is applied to the cluster operator, ensuring the trace-preservation of the density matrix. In addition, the diagrammatic techniques of the traditional coupled-cluster approach are directly transferable, which facilitates the computational implementation of the method.

Investigation of the single-impurity Anderson model demonstrates the accuracy of the approach in reproducing the electronic dynamics of the impurity compared to renormalization-group methods. The coupled-cluster time evolution captures the balance of correlation effects with increasing interaction strength and temperature in closed systems, and also efficiently includes driving and dissipation effects in open systems. Thus, the method provides a unified framework to simulate the non-equilibrium many-body dynamics of both closed and open systems at finite temperature.

This work presents a study of a coupled-cluster formalism with many avenues for further developments. One possibility is the more general definition of the thermal Bogoliubov transformation. Furthermore, the time-dependent approach can serve as a basis for more approximate treatments that explicitly include thermal and dissipative effects, and possibly more general single-particle density matrices can also be accommodated.

Acknowledgements.

P.S. would like to gratefully acknowledge financial support by the German Science Foundation.