Real-time density-matrix coupled-cluster approach for closed and open systems at finite temperature
Philip Shushkov Thomas F. Miller III
quant-ph cond-mat.str-el





[]shushkov@caltech.edu

Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, USA




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.




Real-time density-matrix coupled-cluster approach for closed and open systems at finite temperature
    Thomas F. Miller III
    December 30, 2023
===================================================================================================





§ 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,<cit.> defects in extended systems,<cit.> quantum transport in molecular junctions,<cit.> and chemical reactions at surfaces and interfaces,<cit.> as well as the description of novel driven-dissipative phenomena in open systems, such as discrete time crystals.<cit.> 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.<cit.> 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.<cit.> Furthermore, time-dependent coupled-cluster wavefunction methods have been recently applied to the non-equilibrium electronic dynamics of molecules,<cit.> and beyond the single-reference framework, the coupled-cluster approach has been extended to more challenging, multi-reference systems.<cit.>

Thermal effects have been recently included in the coupled-cluster formalism<cit.> using the method of thermal cumulants.<cit.> Alternatively, thermo-field theory has been combined with the method of configuration interaction to simulate finite-temperature molecular properties.<cit.> The thermo-field method<cit.> was developed to incorporate thermal mixed states in quantum field theory and the approach was later extended to open systems.<cit.> Besides configuration interaction treatment, the method has allowed direct application of coupled-cluster theory to non-equilibrium stationary states,<cit.> and it has been also combined with renormalization-group methods<cit.> 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 model<cit.> 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. <ref>. We present the results for the quantum dynamics of the single-impurity Anderson model at a range of interaction strengths and temperatures in Sec. <ref>, including also dynamics with driving and dissipation. We conclude in Sec. <ref> with remarks on the future directions for method development. 



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



 §.§ Thermo-field theory

Thermo-field theory<cit.> expresses the mixed states of quantum systems in the form of pure state wavefunctions. Let us define a basis in the Hilbert space ℋ of a system by the ket-states |n⟩, generated via the action of the creation operators a^†_i on the vacuum ket |0⟩, |n⟩ = ∏_i a^†_i |0⟩, where n is the set of occupation numbers and the index i enumerates the states of the single-particle basis. The creation operators a^†_i and their Hermitian conjugate annihilation operators a_i obey the commutation relations

    [a_i,a^†_j]_σ = δ_ij, [a^†_i,a^†_j]_σ = [a_i,a_j]_σ = 0,

where the commutator is [A,B]_σ = AB - σ BA with σ=1 for bosons and σ=-1 for fermions. Using the basis |n⟩ of ℋ, the density matrix ρ is given by

    ρ = ∑_n,mρ_nm|n⟩⟨m|

with the ket-bra couples |n⟩⟨m| forming an orthonormal basis. To write Eq. (<ref>) in the form of a wavefunction, thermo-field theory uses a second replica ℋ̃ of the original Hilbert space ℋ to express the ket-bra couples |n⟩⟨m| as ket-basis states |n,m̃⟩, effectively doubling the number of degrees of freedom of the system. We define tilde creation ã^†_i and annihilation ã_i operators in the tilde Hilbert space ℋ̃ via the relations

    ã_i |n⟩⟨m|   = σ^μ+1|n⟩⟨m| a^†_i, 
    ã^†_i |n⟩⟨m|   = σ^μ|n⟩⟨m| a_i

with μ=∑_i (n_i-m_i), the difference of ket- and bra-state occupation numbers. From Eqs. (<ref>) and  (<ref>), the non-tilde and tilde operators satisfy the commutation relations

    [ã_i,ã^†_j]_σ = δ_ij, [ã^†_i,ã^†_j]_σ    = [ã_i,ã_j]_σ = 0, 
    [ã^†_i,a^†_j]_σ    = [ã_i,a_j]_σ = 0,

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

    (AB)^   = ÃB̃, 
    (c_a A + c_b B)^   = c^*_a Ã + c^*_b B̃,

with c_a and c_b complex numbers and * denoting complex conjugation. The double tilde conjugation rules ã̃_i = σ a_i and ã̃^†_i = σ a^†_i supplement Eq. (<ref>).

In the doubled Hilbert space ℋℋ̃ of non-tilde (physical) and tilde (hidden) degrees of freedom, the density matrix, Eq. (<ref>), takes the form of a pure state wavefunction

    ρ = ∑_n,mρ_nm|n,m̃⟩.

The thermal density matrix |ρ_0⟩ of a non-interacting system, for instance, can be written in the exponential form

    |ρ_0⟩ = 1/Q_0∑_n e^-β∑_i n_i ϵ_i|n,ñ⟩ = 1/Q_0exp(∑_i e^-βϵ_i a^†_i ã^†_i ) |0⟩,

where β is the reciprocal temperature, β=(k_BT)^-1 with T the temperature, Q_0 is the partition function, ϵ_i are the single-particle energy levels, and |0⟩ = |0,0̃⟩ is the vacuum ket-state in the doubled Hilbert space. Similarly, the resolution of identity, the unit ket-state |1⟩, also has an exponential form

    |1⟩ = ∑_n |n,ñ⟩ = exp( ∑_i a^†_i ã^†_i ) |0⟩.


The average value of an operator O, given by the trace with the density matrix ρ, becomes an expectation value with the ket-state |ρ⟩, Eq. (<ref>), in thermo-field theory. Specifically,

    tr[ρ O ] = ⟨1| O |ρ⟩,

where the density matrix is normalized ⟨1||ρ⟩=1. The bra-state and ket-state in Eq. (<ref>) are different due to the choice of the path contour in the complex time plane as shown in Sec. <ref>. 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. <ref> devoted to the coupled-cluster formalism.

The Liouville-von Neumann equation for the time-evolving density matrix, iρ = [H,ρ], becomes Schrödinger equation for the pure ket-state |ρ⟩

    i|ρ_t⟩ = Ĥ|ρ_t⟩,

where the super-Hamiltonian operator Ĥ is the difference of the original Hamiltonian H and its tilde conjugate

    Ĥ = H - H̃.

The super-Hamiltonian in Eq. (<ref>) 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. <ref>.

Formal integration of Eq. (<ref>) gives the time-dependent density matrix |ρ_t⟩ as

    |ρ_t⟩ = e^-i Ĥ t|ρ_t_0⟩

with |ρ_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. <ref>, 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. (<ref>). Non-equilibrium states relax to stationary states, such as equilibrium states, with time-independent density matrices |ρ_s⟩, which are annihilated by the super-Hamiltonian Ĥ|ρ_s⟩=0.



 §.§ Thermal quasi-particles

The creation a^†_i and ã^†_i and annihilation a_i and ã_i operators have a complementary action on the unit bra-state

    ⟨1| a^†_i    = ⟨1|ã_i, 
    ⟨1|ã^†_i    = σ⟨1| a_i.

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

    ⟨1|Ĥ = 0.

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

    ⟨1| O |ρ_t⟩ = ⟨1| O_t |ρ_t_0⟩

with O_t = e^i Ĥ t O e^-i Ĥ t – the Heisenberg-evolved operator O.

The action of the a_i and ã_i operators on the non-interacting thermal density matrix is similarly complementary

    a_i |ρ_0⟩   = e^-βϵ_iã^†_i |ρ_0⟩, 
    ã_i |ρ_0⟩   = σ e^-βϵ_i a^†_i |ρ_0⟩.

Because a_i annihilates an electron, the thermal ket-state condition, Eq. (<ref>), shows that ã^†_i creates a hole, and vice versa because a^†_i creates an electron, ã_i annihilates a hole.

Eqs. (<ref>) and  (<ref>) imply that we can transform the creation a^†_i and ã^†_i and annihilation a_i and ã_i operators to a new set of creation b^†_i and b̃^†_i and annihilation b_i and b̃_i operators, which respectively annihilate the unit bra-state and non-interacting thermal ket-state

    ⟨1| b^†_i    = ⟨1|b̃^†_i = 0, 
    
            b_i |ρ_0⟩   = b̃_i |ρ_0⟩ = 0.

The transformation that satisfies Eq. (<ref>) and preserves the commutation relations, Eqs. (<ref>) and  (<ref>), is given by 

    b^†_i = a^†_i - ã_i,      b̃^†_i = ã^†_i - σ a_i, 
    
        b_i = u_i a_i - v_i ã^†_i,      b̃_i = u_i ã_i - σ v_i a^†_i,

where v_i = ( e^βϵ_i - σ)^-1 is the Bose-Einstein or the Fermi-Dirac distribution and u_i = 1 + σ v_i. Eq. (<ref>) 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.<cit.> The thermal Bogoliubov transformation, Eq. (<ref>), defines new thermal quasi-particles b^†_i and b_i and thermal quasi-holes b̃^†_i and b̃_i with respect to the unit bra-state and non-interacting thermal ket-state, Eq. (<ref>), which we refer to as the new left ⟨0| and right |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.



 §.§ Density-matrix coupled-cluster approach

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

    |ρ_t⟩ = e^T_t|ρ_t_0⟩,

where |ρ_t_0⟩=|ρ_0⟩ is the non-interacting thermal density matrix in Eq. (<ref>), and T_t is the time-dependent cluster operator, which we expand in the number of excited thermal quasi-particle-hole pairs n as

    T_t = ∑^N_n=1 T_n,

where N is the highest excitation level. We express the cluster operators T_n via the quasi-particle b^†_i and quasi-hole b̃^†_i creation operators from Sec. <ref> as

    T_n = 1/n!∑_ij ... , kl ...  t^ik..._jl... b^†_i b^†_k ... b̃^†_l b̃^†_j,

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. (<ref>). The explicit form of the cluster operator for this singles and doubles approximation is

    T_t = ∑_ij t^i_j b^†_i b̃^†_j + 1/4∑_ijkl t^ik_jl b^†_i b^†_k b̃^†_l b̃^†_j.


Eq. (<ref>) together with Eqs. (<ref>) and  (<ref>) 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^† and b̃^†_i creation operators, Eq. (<ref>), ⟨1||ρ_t⟩ = ⟨1| e^T_t|ρ_0⟩ = ⟨1||ρ_0⟩. The time evolution governed by the Schrödinger equation, Eq. (<ref>), 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, |ρ̃_̃t̃⟩ = |ρ_t⟩, resulting in a tilde invariant cluster operator 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.<cit.> Thermal averages of Heisenberg-evolved operators O(t) are given in general by the ratio

    G = tr[e^-β H O(t_N) .. O(t_1)]/tr[e^-β H] = ⟨ T_c U_c O_N .. O_1 ⟩/⟨ T_c U_c ⟩.

In diagrammatic terms, the denominator in Eq. (<ref>) cancels the disconnected vacuum diagrams in the numerator, such that only connected diagrams contribute to the thermal average. In Eq. (<ref>), 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β. We can take respectively the initial t_i and final t_f times to the infinite past, t_i → -∞, and infinite future, t_f →∞, 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 theory<cit.> and the density-matrix coupled-cluster theory presented here. 

Using the definition of expectation value Eq. (<ref>) and the super-Hamiltonian Eq. (<ref>), we cast Eq. (<ref>) in thermo-field theory as

    G = ⟨1| Te^-iĤt O_N .. O_1 |ρ_0⟩/⟨1|e^-iĤt|ρ_0⟩ = ⟨1| O(t_N) .. O(t_1) |ρ_0⟩,

where the imaginary axis propagation is represented by the non-interacting thermal state |ρ_0⟩ and the contour runs only along the real time axis from t_i to t_f. The second equality in Eq. (<ref>) 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β/2 rather than immediately below the real time axis.<cit.> 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. (<ref>) in the Schrödinger equation Eq. (<ref>) for the density matrix ket-state.<cit.> Following multiplication from the left with e^-T_t, we obtain

    ie^-T_t∂_te^T_t|ρ_0⟩ = e^-T_tĤe^T_t|ρ_0⟩,

where Ĥ is the super-Hamiltonian of the system expressed in the thermal quasi-particle representation. Because ∂_t is a single-particle operator and the creation quasi-particle operators are time-independent, the Baker-Campbell-Hausdorff identity gives

    e^-T_t∂_te^T_t = ∂_t + T_t,

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. (<ref>) and  (<ref>), commute, application of the Baker-Campbell-Hausdorff identity shows that the expansion of the similarity-transformed super-Hamiltonian e^-T_tĤ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. (<ref>), simplifies to

    i T_t |ρ_0⟩ = ( Ĥ∑^n=N_n=01/n! T^n_t )_c |ρ_0⟩,

where N stands for the largest order non-vanishing term in the series, and we used the time-independence of the ket-state ∂_t |ρ_0⟩ = 0. The subscript "c" in Eq. (<ref>) 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. (<ref>) onto states with finite number of quasi-particle excitations, ⟨0|b̃_j .. b_i,

    i ⟨0|b̃_j .. b_i T_t |ρ_0⟩ = ⟨0|b̃_j .. b_i ( Ĥ∑^n=N_n=01/n! T^n_t )_c |ρ_0⟩.

In our singles and doubles approximation to the DMCC method, Eq. (<ref>), 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.<cit.>

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.<cit.> 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^†_i and b̃^†_i annihilate the left vacuum ⟨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 diagrams<cit.> 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.



 §.§ 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.<cit.> The SIAM Hamiltonian H=H_0+W has a non-interacting part H_0 given by

    H_0 = ∑^N_b_i=0sϵ_i a^†_is a_is + ∑^N_b_i=1,s( V_i a^†_0s a_is + a^†_is a_0s),

where ϵ_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

    W = U a^†_0αa_0αa^†_0βa_0β

with U the Hubbard parameter, a measure of the interaction strength. Eqs. (<ref>) and  (<ref>) represent a closed system with a super-Hamiltonian that is the difference of H and its tilde conjugate, Ĥ = H - H̃ = H_0 - H̃_0 + W - 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,<cit.> such as

    D̂ = -i ∑^N_b_i=1,s    [ (γ^(1)_i - γ^(2)_i) ( a^†_is a_is + ã^†_isã_is) 
        - 2 γ^(1)_iã_is a_is + 2 γ^(2)_iã^†_is a^†_is + 2 γ^(2)_i ],

with γ^(1)_i = γ (1-v_i), γ^(2)_i = γ v_i, and γ the dissipation strength. The relation among γ^(1)_i, γ^(2)_i and γ follows from the detailed balance condition, and D̂ conforms with the dissipative part of a Lindblad operator<cit.> 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 = ϵ_0 + δϵsin(Ω t). We use logarithmic discretization of the electronic bath in our calculations as suggested by Wilson for the Kondo model<cit.> with a discretization parameter Λ=1.1 and N_b = 100. We consider also the more challenging asymmetric case of SIAM with ϵ_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. (<ref>), the SIAM super-Hamiltonian Ĥ = Ĥ_0 + Ŵ + D̂ becomes Ĥ^̂'̂ = Ĥ^̂'̂_̂0̂ + Ŵ^̂'̂, where the non-interacting super-Hamiltonian Ĥ^̂'̂_̂0̂ is given by

    Ĥ^̂'̂_̂0̂ = ∑_i,j,s h_ij b^†_is b_js - h^*_ijsb̃^†_isb̃_js + Δ_ij( b^†_isb̃^†_js - b^†_jsb̃^†_is),

with the single-particle Hamiltonian matrix h_ijs,

    h_ij = ( ϵ + v_0 U )    δ_ij=0 + ϵ^'_i δ_ij>0
        + V_i δ_i>0, j=0 + V_j δ_i=0,j>0,

and the pairing matrix Δ_ijs,

    Δ_ij = V_i (u_i v_0 - u_0 v_i) δ_i>0, j=0 + V_j (u_j v_0 - u_0 v_j) δ_i=0, j>0,

and the interaction Ŵ^̂'̂ is written as

    Ŵ^̂'̂ =     U (u_0-v_0) b^†_α b_α b^†_β b_β + U b^†_α b_αb̃_β b_β + 
        U b̃_α b_α b^†_β b_β + U u_0 v_0 b^†_α b_α b^†_βb̃^†_β + 
        U u_0 v_0 b^†_αb̃^†_α b^†_β b_β - t.c.

with t.c. standing for tilde conjugation. In Eq. (<ref>), ϵ^'_i = ϵ_i - i γ 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_0U. The Hubbard interaction in Eq. (<ref>) also incorporates thermal effects as the interaction strength depends on the thermal population of the impurity. Because Ŵ^̂'̂ includes contributions with up to three annihilation operators, the series in the equation of motion for the cluster operator, Eq. (<ref>), 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.  



 §.§ 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 methods<cit.> to closed and open systems at finite temperature. We apply the time-evolving block-decimation approach<cit.> to SIAM and compare the results to DMCC in Sec. <ref>. 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 α impurity site, followed by N_b α bath sites, the β impurity site, and N_b β bath sites. Each electronic site carries both non-tilde and tilde degrees of freedom. The terms of the SIAM super-Hamiltonian Ĥ^̂'̂ act on two sites at a time with the interaction Ŵ^̂'̂, Eq. (<ref>), coupling the impurity sites and the non-interacting super-Hamiltonian Ĥ^̂'̂_̂0̂, Eq. (<ref>), 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 Ŝ to rearrange sites given by

    Ŝ_i+1,i = s_i+1,i,α s_i+1,i,β - s̃_i+1,i,αs̃_i+1,i,β,

where

    s_pqσ = 1 + a^†_p,σ a_q,σ + a^†_q,σ a_p,σ - a^†_p,σ a_p,σ - a^†_q,σ a_q,σ.

We apply the second-order Suzuki-Trotter decomposition<cit.> 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 δ t

    e^-iĤ^̂'̂δ t/2 = ∏_i ∈α[e^-i(Ĥ^̂'̂_̂0̂)_0iδ t/2Ŝ_i,i+1]     e^-iŴ^̂'̂δ t/2
       ∏_j ∈β[e^-i(Ĥ^̂'̂_̂0̂)_0jδ t/2Ŝ_j,j+1],

where i runs over the sequence of α sites and j over the sequence of β sites, respectively.
The full evolution operator is the product of Eq. (<ref>) and its transpose. The application of each of the two-site gates in Eq. (<ref>) 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.



§ 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  <ref> 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 <ref> 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. (<ref>), and a harmonic driving contribution, H_t = ϵ_0 + δϵsin (Ω t) with δϵ = 0.08 eV. Figure <ref> 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, Ω=π V and Ω=4 π V. Figure <ref>A also includes dissipation with strength γ=0.04 eV, and Figure <ref>B – dissipation with strength γ=0.2 eV.
Figure <ref>A 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 <ref>A and <ref>B 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 <ref> 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.



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


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


*


