Thermofield theory for finite-temperature coupled cluster
Gaurav Harsha Thomas M. Henderson Gustavo E. Scuseria
physics.chem-ph cond-mat.str-el



Department of Physics and Astronomy, Rice University, Houston TX 77005

Department of Physics and Astronomy, Rice University, Houston TX 77005
Department of Chemistry, Rice University, Houston TX 77005

Department of Physics and Astronomy, Rice University, Houston TX 77005
Department of Chemistry, Rice University, Houston TX 77005


    We present a coupled cluster and linear response theory to compute properties of many-electron systems at non-zero temperatures. For this purpose, we make use of the thermofield dynamics, which allows for a compact wavefunction representation of the thermal density matrix, and extend our recently developed framework [J. Chem. Phys. 150, 154109 (2019)] to parameterize the so-called thermal state using an exponential ansatz with cluster operators that create thermal quasiparticle excitations on a mean-field reference. As benchmark examples, we apply this method to both model (one-dimensional Hubbard and Pairing) as well as ab-initio (atomic Beryllium and molecular Hydrogen) systems, while comparing with exact results.


Thermofield theory for finite-temperature coupled cluster
    Gustavo E. Scuseria
    July 24, 2019
=========================================================





§ INTRODUCTION

Temperature is an important parameter of physical systems. For many problems, the temperature scale of interest is far below the optical gap (or the excitation energies). For example, typical electronic excitation energies in molecules are of the order of a few electron volts (or eV) , which is much larger than room temperature (∼ 25 meV). In such problems where we are interested only in electronic degrees of freedom, it suffices to know the ground electronic state and perhaps a few low-lying excited states. To access these states, we solve the time-independent Schrödinger equation. As is well known,<cit.> this is a very complicated problem and one generally relies on a series of approximate methods such as Hartree-Fock, perturbation theory, configuration interaction (CI), coupled cluster theory<cit.> (CC), Monte Carlo methods,<cit.> among many more.<cit.>

There are, however, many interesting problems and applications where one may be interested in temperature scales that are comparable or even larger than the excitation gap. Examples include metallic compounds with small gap that host an unconventional superconductivity persisting at relatively high temperatures,<cit.> ultra-cold chemistry,<cit.> geochemical processes which generally involve very high temperatures and pressure,<cit.> etc. In these problems, we can no longer make do with a few electronic states and must evaluate properties as thermal averages weighted over an appropriate ensemble of states. For a system in thermal equilibrium at inverse temperature β, the ensemble of choice is generally canonical or grand-canonical and the expectation value of an observable A is defined by

    ⟨ A ⟩_β = 1/𝒵Tr ( A ρ ),

where ρ is the thermal density matrix, and 𝒵 is the partition function, given by

    ρ = e^-β H',    𝒵 = Tr ( ρ ),

with H' = H for the canonical ensemble and H' = H-μ N for the grand canonical ensemble, where μ is the chemical potential. In this paper, we shall work explicitly with the grand canonical ensemble.

Exactly computing ρ (or equivalently 𝒵) requires information about the entire spectrum of the Hamiltonian, which is far from feasible. Accordingly, just as for zero temperature, a series of approximate methods is needed. Several methods have been proposed over the years to evaluate thermal averages of operators.<cit.> Most of these methods can be broadly categorized into deterministic methods such as diagrammatic perturbation theory based on the Matsubara formalism,<cit.> and stochastic quantum Monte Carlo methods.
Wavefunction methods are particularly convenient in the study of zero-temperature ground-state properties of finite sized systems and clearly, their thermal equivalents are highly desirable. However, the development of such thermal wavefunctions methods has been rather challenging, primarily because the thermal density matrix cannot be expressed in terms of a single wavefunction in the original Hilbert space. Nevertheless, several such methods have been introduced over the years such as the Ancilla density matrix renormalization group<cit.> and finite-temperature perturbation theories,<cit.> among many more.

Given its features, especially size-extensivity and success with weakly correlated systems, the coupled cluster ansatz is an ideal candidate to study finite-temperature properties.
A thermal analogue of the CC method<cit.> was proposed by Mukherjee et. al. and has been further elaborated recently in independent works by White et al.<cit.> and Hummel.<cit.> This formulation uses a thermal Wick's theorem to compactly represent the imaginary time evolution operator as a thermal normal ordered exponential of some cluster operator and a number.

In this paper, we present an alternative approach to thermal coupled cluster based on the thermofield dynamics (TFD),<cit.> using a framework we recently explored in Ref. . Thermofield dynamics provides a convenient way to represent the thermal density matrix via a wavefunction which evolves in temperature according to the imaginary-time evolution Schrödinger equation. Undoubtedly, TFD has the potential to study many-electron systems in quantum chemistry and condensed-matter physics.<cit.> A connection between thermofield formalism and thermal cluster-cumulant theory was also explored by Mukherjee et. al. in Ref. . Here, we parametrize this so-called thermal wavefunction as an exponential ansatz to integrate the evolution equation.



§ COUPLED CLUSTER THEORY


The coupled cluster method is one of the most widely used methods in quantum chemistry. Introduced first in nuclear physics by Coester and Kümmel<cit.> and later reformulated for electronic structure theory by Čižek and Paldus,<cit.> it uses an exponential wavefunction ansatz

    |Ψ⟩ = e^T |Φ⟩,

to solve the time-independent Schrödinger equation. Here |Φ⟩ is some Hartree-Fock (HF) Slater determinant reference and T contains particle-hole excitations defined on this reference state,

    
    T    = T_1 + T_2 + …, 
    
            T_1    = ∑_i,a t_i^a c^†_a c_i,
    
            T_2    = 1/4∑_i,j,a,b t_ij^ab c^†_a c^†_b c_j c_i, …


where we have followed the standard notation for labeling orbital indices, i.e., occupied orbitals are denoted by indices i,j,k,…, while unoccupied orbitals are denoted by a, b, c, …. Since the wavefunction |Ψ⟩ is constructed as excitations on a single Slater determinant |Φ⟩, the traditional CC is also known as single-reference CC. In what follows, unless explicitly mentioned, we consider only the single-reference CC. The equations for the ground state energy and the unknown t-amplitudes can be obtained by left-projecting the Schrödinger equation

    H̅| Φ⟩ = E  | Φ⟩,   H̅ = e^-T H e^T

with various Slater determinants. For example, if the cluster operator T is truncated to single and double excitations only (CCSD), the energy and amplitude equations are obtained by solving the following equations

    E    = ⟨Φ|H̅|Φ⟩, 
    
        0    = ⟨Φ_i^a|H̅|Φ⟩, 
    
        0    = ⟨Φ_ij^ab|H̅|Φ⟩,

where ⟨Φ_i^a| and ⟨Φ_ij^ab| are singly- and doubly-excited slater determinants.

As introduced so far, CC describes a correlated ansatz for the ket wavefunction. In order to compute expectation values other than that of the Hamiltonian, one also needs a correlated bra state. A linear response wavefunction is generally employed for this purpose,<cit.> i.e. one makes the energy functional

    E = ⟨Φ| (1 + Z) e^-T H e^T|Φ⟩

stationary with respect to T and Z, where

    
    
    Z    = Z_1 + Z_2 + …, 
    
            Z_1    = ∑_i,a z_i^a c_i^† c_a, 
    
            Z_2    = 1/4∑_i,a z_ij^ab c_i^† c_j^† c_b c_a.


By realizing that the cluster operator T and CI operator Z are composed of particle-hole excitation and de-excitation operators respectively, the bra state can be re-written as an explicit CI wavefunction

    ⟨Ψ'| = ⟨Φ| (1 + Z) e^-T = ⟨Φ| (1 + W) e^w_0,

where w_0 is a constant and W has the same operator-form as Z.

A similar formulation known as equation of motion CC<cit.> (EOM-CC) can be used for excited states. Coupled cluster gives highly accurate results for weakly correlated systems and its success can be attributed to the polynomial computational scaling (𝒪(N^6) for CCSD) as well as to the fact that the computed properties are size-extensive (scale linearly with the particle number) in the thermodynamic limit.



§ THERMOFIELD DYNAMICS

Thermofield dynamics is a real-time thermal field theory that treats both time and temperature on equal footing, and was proposed as an alternative to the Keldysh approach in the Matsubara imaginary time formalism. It provides a prescription for purification of the thermal density matrix, allowing us to construct a wavefunction |Ψ(α, β) ⟩, with α=βμ, generally known as the thermal vacuum, thermofield double state or simply as the thermal state, so that the trace over an ensemble of states in Eq. <ref> can be replaced by an expectation value over this wavefunction, i.e.

    ⟨ A ⟩ = ⟨Ψ (α, β) | A |Ψ (α, β) ⟩/⟨Ψ (α, β) |Ψ (α, β) ⟩.

That the thermal state |Ψ (α, β) ⟩ cannot be a pure state in the physical Hilbert space ℋ is easily established since the density matrix ρ represents mixed states to begin with. In TFD, one therefore introduces a fictitious, conjugate copy of the original Hilbert space, known as the tilde-conjugate space or ℋ̃, and the thermal state is then defined in the doubled space ℋ⊗ℋ̃ as

    |Ψ (α, β) ⟩ = e^(α N-β H)/2|𝕀⟩,   |𝕀⟩ = ∑_m | m, m̃⟩,

where H is the Hamiltonian, {| m ⟩} is some orthonormal basis in the Hilbert space ℋ and |m̃⟩ is the tilde-state corresponding to | m ⟩, while we have used the shorthand notation | m, m̃⟩ = | m ⟩⊗|m̃⟩. The state |𝕀⟩ is the state with maximal entanglement between ℋ and ℋ̃ and is the exact thermal state at infinite temperature (i.e. β=0) and α=0. Moreover, it is invariant under any transformation of the basis. Further details about the TFD formalism and the structure of the tilde-conjugate space ℋ̃ can be found in Ref. harsha_thermofield_2019 and the references therein.

By construction, the thermal state satisfies the following imaginary time evolution equations,

    
    
    ∂/∂β|Ψ (α, β) ⟩   = -1/2 H |Ψ (α, β) ⟩.
            
    ∂/∂α|Ψ (α, β) ⟩   = 1/2 N |Ψ (α, β) ⟩.


One can solve for |Ψ (α, β)⟩ by integrating Eq. <ref>, generally starting from (α,β)=(0,0), where the initial thermal state is known exactly.
Exactly evolving the thermal state |Ψ (α, β)⟩ is equivalent to the exact diagonalization of the density matrix ρ. Clearly, approximations need to be introduced in the process.

The simplest approximation that we can invoke is to use the mean-field Hamiltonian H_0 instead of H. For a many-electron system, the state |𝕀⟩ can be expressed in terms of single-particle Fock states,

    |𝕀⟩ = ∏_p∈levels ( | 0,0̃⟩_p + | 1,1̃⟩_p  ),

where | 0 ⟩_p and | 1 ⟩_p mean that the orbital p  is empty or occupied, respectively.
Moreover, if we chose to work with the eigen-basis of H_0, i.e.

    H_0 = ∑_p ϵ_p c^†_p c_p

the normalized mean-field thermal state can be written as

    | 0 (α, β) ⟩ = ∏_p ∈levels ( x_p | 0, 0̃⟩_p + y_p | 1, 1̃⟩ ),

where x_p and y_p are related to the Fermi-Dirac statistics,

    
    x_p    = 1/√(1 + e^(α - βϵ_p)), 
    
            y_p    =  e^(α - βϵ_p)/2/√(1 + e^(α - βϵ_p)),


with x_p^2 + y_p^2=1. The  mean-field thermal state in Eq. <ref> allows us to introduce a thermal Bogoliubov transformation,

    [   a_p; ã_p^† ] =
        [  x_p -y_p;  y_p  x_p ][    c_p; c̃_p^† ],

such that

    a_p | 0 (α, β) ⟩ = 0 = ã_p | 0(α, β) ⟩.

Correlated methods can be built with either |𝕀⟩ or | 0(α, β) ⟩ as the reference, while integrating Eq. <ref>. As we have discussed in Ref , the former choice, which we call the fixed-reference formalism, performs well only in the vicinity of β=0, while the latter, called the covariant formalism, yields accurate results for the entire range of β. In the next section, we present details to integrate Eq. <ref> through the CC ansatz in the covariant formalism. We have also included a short discussion on the fixed-reference approach in Appendix <ref>. For brevity of notation, we will henceforth refer the  mean-field thermal state | 0(α, β)⟩ by |Ψ_0 ⟩.



§ THERMAL COUPLED CLUSTER


As explained in Eq. <ref>, the CC expectation value of any operator A can be evaluated as an asymmetric expectation value,

    ⟨ A ⟩_CC = ⟨Ψ' | A |Ψ⟩/⟨Ψ' |Ψ⟩,

where both the ket |Ψ⟩ and the bra ⟨Ψ' | states are approximations to the same thermal state, and consequently evolve according to Eq. <ref> and its adjoint respectively.

Given that the Bogoliubov transformation in Eq. <ref> is BCS-like, i.e. orbitals p (in ℋ) and p̃ (in ℋ̃) in the mean-field thermal state (Eq. <ref>) are coupled in just the same way as two opposite momentum single-particle levels in a BCS-wavefunction, we parametrize the ket state as an exponential of quasiparticle creation operators<cit.> acting on an α- and β-dependent mean-field thermal reference, |Ψ_0 ⟩, defined in Eq. <ref>,

    
    |Ψ⟩   = e^S(α, β)|Ψ_0 ⟩,
            
    
            S    = s_0
            + ∑_p,q s_pq a^†_p ã^†_p
            + 1/(2!)^2∑_p,q,r,s s_pqrs a^†_p a^†_q ã^†_s ã^†_r
            + …


On the other hand, the bra state, as in the traditional CC formalism, is approximated as a linear CI-like wavefunction, i.e.

    ⟨Ψ'| = ⟨Ψ_0 |  (1 + Z)e^z_0e^-S,

and as explained in Eq. <ref>, it can be expressed as an effective CI wavefunction

    
    ⟨Ψ'|   = ⟨Ψ_0 | (1 + W) e^w_0, 
    
            W    = ∑_p,q w_pqã_q a_p + 1/4∑_p,q,r,s w_pqrsã_r ã_s a_q a_p + …


For both the bra and the ket states, the reference |Ψ_0 ⟩ evolves continuously as we evolve the Schrödinger Eq. <ref> (hence the name `covariant'). Accordingly, both the amplitudes (s_pq, w_pq, etc.) and the quasiparticle operators a^†, ã^† carry α- and β-dependence.
The α- and β-evolution of ⟨Ψ' |, a thermal CI wavefunction, is governed by

    
    
    ⟨Ψ_0 | ( ∂ W/∂α + (1 + W) ∂ w_0/∂α )
               = 1/2⟨Ψ_0 |  N_CI,
            
    ⟨Ψ_0 | ( ∂ W/∂β + (1 + W) ∂ w_0/∂β )
               = -1/2⟨Ψ_0 | H_CI,


where H_CI and N_CI are effective CI Hamiltonian and Number operators respectively, and are given by

    H_CI   = (1 + W) H - H_0 (1 + W) ,
        
    
        N_CI   = W  N - N  W.

A detailed discussion on thermal CI and the derivation of these equations can be found in Ref. harsha_thermofield_2019.

For the evolution of |Ψ⟩, substituting the CC ansatz from Eq. <ref> into the Schrödinger Eq. <ref> gives

    
    
    e^-S ( ∂/∂α e^S ) |Ψ_0 ⟩   =
            1/2 ( e^-S N e^S - N  )  |Ψ_0 ⟩,
            
    
            e^-S ( ∂/∂β e^S ) |Ψ_0 ⟩   =
            - 1/2 ( e^-S H e^S - H_0  )  |Ψ_0 ⟩.


The evolution equations for the amplitudes can be obtained by left projecting Eq. <ref> with the respective determinants.



 §.§ Wilcox identity

The process of reducing Eq. <ref> to evolution equations for the amplitudes is complicated by the fact that the derivative of the cluster operator does not commute with the operator itself, i.e.

    [ ∂ S/∂ x , S  ] ≠ 0,

where x = α, β. The derivative of the exponential cluster operator is appropriately performed by making use of the Wilcox identity,<cit.> which states that the derivative of the exponential of an operator M with respect to some parameter λ can be evaluated as

    ∂/∂λe^M(λ) = ∫_0^1 dy   e^(1-y)M∂ M/∂λ e^y M.

With this, the left-hand side of Eq. <ref> becomes

    
    e^-S ( ∂_x   e^S )    = ∫_0^1 dy   e^-y S  (∂_x S)    e^y S, 
       = (∂_x S) + 1/2! [ (∂_x S), S ] 
          + 1/3![ [(∂_x S), S], S ] + …,


where we have used the shorthand ∂_x for ∂/∂ x, and made use of the Baker-Campbell-Hausdorff expansion in going from Eq. <ref> to <ref>. Finally, breaking the derivative ∂_x S into the amplitude (∂_ampS) and operator (∂_op S) derivatives,

    ∂_x S = ∂_ampS + ∂_op S,

and realizing that the former commutes with S, we can compactly write the left-hand side of Eq. <ref> as

    e^-S ( ∂_x   e^S ) = ∂_amp S/∂ x+ S_x,   
        S_x = ∫_0^1 dy   e^-yS∂_op S/∂ x   e^yS,

where the integration over y in the second equation can be carried out analytically, as explained in Eq. <ref>. With these details, Eq. <ref> can be further simplified as

    
    
    ∂_amp S/∂α |Ψ_0 ⟩   =  [ 1/2 ( e^-S N e^S - N  ) - S_α ]  |Ψ_0 ⟩
    ∂_amp S/∂β |Ψ_0 ⟩   = - [ 1/2 ( e^-S H e^S - H_0  ) + S_β ]  |Ψ_0 ⟩,


which can then be left-projected with various thermal quasiparticle states to yield a set of differential equations governing the evolution of the s-amplitudes in the chemical potential - temperature or α-β space. Complete expressions for the CCSD evolution equations are included in Appendix <ref>.









§ IMPLEMENTATION



 §.§ Integration

An explicit Runge-Kutta method of order (4)5 with step-size control<cit.> has been employed for integrating the resulting set of differential equations, where step-size is adaptively modified so as to keep the relative error in evolution within 10^-8. Furthermore, for applications that require a fixed average filling fraction, a tolerance of 10^-5 in the number of electrons is used to find the target chemical potential.



 §.§ Choice of H_0

As mentioned earlier, for the covariant thermal CC, a mean-field Hamiltonian H_0 = ∑_p _p a_p^† a_p is used to construct the mean-field thermal reference state |Ψ_0 ⟩. There are several alternatives to pick this mean-field Hamiltonian – in thermal HF theory,<cit.> one uses an H_0 that optimizes the applicable free energy, and it therefore depends on α and β. Such a choice, while it may yield an excellent mean-field thermal reference state, will result in convoluted expressions for the mean-field evolution as well as the operator derivatives, and is beyond the scope of current work.

In our implementation, we use zero-temperature ground-state HF, which may or may not preserve the ground-state symmetries, to construct the energy eigenvalues in H_0. Using an α- and β-independent H_0 is analogous to the imaginary-time interaction picture formalism and simplifies our implementation since, (i) an α- and β-independent H_0 leads to clean analytical forms for the mean-field evolution operator as well as the α- and β-derivatives of the thermal quasiparticle operators, and (ii) H_0 being diagonal yields a straightforward thermal Bogoliubov transformation in Eq <ref>. We note, however, that this choice is different from the one used in Ref. harsha_thermofield_2019, where we use simply the one-electron Hamiltonian to construct H_0.



§ RESULTS


Armed with the working equations, we now proceed to present results for the application of the thermal CC in the covariant formalism, truncated at singles and doubles (CCSD), to various many-electron systems, viz. the one-dimensional Hubbard model,<cit.> the pairing or the reduced BCS model, as well as chemical systems (atomic Beryllium and molecular H_2).
In order to make correspondence with the canonical ground state limit, we present results for all of these systems with a fixed number of particles on average - at each β grid-point, we evolve the thermal states in α to fix the average number of particles before evolving again in β.
We compare our results with full configuration interaction (FCI) results.

We first apply the CCSD methods to the one-dimensional Hubbard model with periodic boundary conditions. Having already presented results for thermal CI truncated to singles and doubles (CISD) in Ref. harsha_thermofield_2019, this model system seems to be the right place to start comparing thermal CC with thermal HF and CI. The Hamiltonian is given by

    H = -t ∑_⟨ p,q ⟩,  σ(c_p,σ^†  c_q,σ + h.c.) + U  ∑_p n_p,↑  n_p,↓,

where ⟨,⟩ denotes that the sum is carried over sites connected in the lattice, t denotes the strength of the kinetic energy term, U denotes the strength of the on-site Coulomb repulsion, and n_p,σ = c_p,σ^†  c_p,σ is the number operator for lattice site p and spin σ.  The ratio U/t characterizes the correlation strength.




Figure <ref> shows the temperature dependence of the error in internal energy for a two-site Hubbard model at half-filling on average with U/t=1 as computed by thermal HF, covariant CISD, thermal CCSD as well its fixed-reference formulation. Figure <ref> presents the same for a six-site Hubbard model with U/t=2. For both cases, restricted HF (or RHF) eigenvalues and integrals have been used to construct H_0 and H respectively. Thermal CCSD clearly outperforms CISD, especially for larger systems where CI is a less accurate wavefunction ansatz. Moreover, the covariant thermal CISD and CCSD go to the appropriate ground-state restricted CCSD in the limit β→∞ (or θ→ 0, where θ = 1/β is the temperature). We note that while CCSD and CISD are exact in describing the ground state of the two-site Hubbard model, they are not exact at finite temperatures since we are working in the grand canonical ensemble.
The fixed-reference CCSD, on the other hand, performs poorly except for a small window around β = 0. Such a behaviour can be anticipated since this method uses the thermal reference corresponding to β = 0 as its starting point. Accordingly, for all other results that follow, we present only the covariant methods.

Next, we consider the reduced BCS or the pairing model, the Hamiltonian for which is given by

    H = ∑_p ϵ_p N_p - G ∑_p,q P^†_p P_q

where N_p counts the number of electrons, ϵ_p denotes the energy, and P^†_p/P_p respectively creates/annihilates a pair of electrons in the p^th-level, while G quantifies the attractive pair-hopping interaction. Here we choose the energy levels with a uniform spacing of 1 unit, i.e., Δϵ = ϵ_p+1-ϵ_p = 1.
Figure <ref> describes the temperature dependence of the error in internal energy for a six-level pairing model with G = 0.2 (weakly correlated) at average half-filling. Figure <ref> shows the same for G = 0.5 (near critical regime). Again, we use RHF eigenvalues and integrals to construct H_0 and H, and we see that the covariant thermal CCSD improves significantly over both the HF and CISD.
We also recover the zero-temperature ground-state limit for thermal CISD and CCSD.
Figures <ref> and <ref> show similar trends for atomic Beryllium and molecular H_2 at bond length 0.74 in STO-3G basis sets.






In all of our applications so far, H_0 is constructed using RHF energy eigenvalues. Accordingly, in the zero-temperature limit, thermal HF approaches ground-state RHF and thermal CCSD approaches ground-state RCCSD. In strongly-correlated regime, where zero-temperature RCCSD fails to converge, the evolution of thermal CCSD also eventually diverges. Such an issue can be avoided by using unrestricted HF (or UHF) energy eigenvalues and integrals to construct H_0 and H respectively. Figure <ref> shows the performance of RHF and UHF based thermal CCSD internal energies against FCI results for six-site Hubbard model at half-filling on average with U/t=5.

In addition to the internal energy, we can also compute other physical properties and correlation functions at any temperature / chemical potential. In Figure <ref>, we show the z-component spin-spin correlation function

    χ(i,j) = ⟨ S^z(i) S^z(j) ⟩

for the ten-site Hubbard model with U/t = 2. Here, the expectation values are computed using linear-response density matrices and orbital relaxation effects have not been considered. As one would expect, at very high temperature θ, there is no correlation between adjacent spins. As θ is reduced, the correlation appears and becomes maximal in the zero-temperature limit. Again, as with the internal energies, we see that the correlation function approaches the ground state CC in the limit β→∞.



§ CONCLUSIONS


We have demonstrated that the framework of thermofield dynamics can be exploited to formulate a finite-temperature coupled cluster theory.
We use the CCSD approximation to benchmark our method on various many-electron Hamiltonians and find that it performs substantially better than the thermal Hartree-Fock and thermal CI, just as one would expect for their ground-state counterparts.
This improved accuracy comes with the same asymptotic 𝒪(N^6) scaling (N being the number of spin-orbitals or basis functions) as does standard quasiparticle CCSD, with a modestly larger pre-factor, though note that we must solve these equations at each grid point in the imaginary-time evolution.
We also observe that in the zero-temperature limit, thermal HF and CCSD approaches the appropriate ground-state HF and CCSD.

The thermofield based formulation of CC provides the basis for generalization of more sophisticated CC ansatze, e.g. unitary and multi-reference CC, to finite-temperatures, something that will be explored in future work.
Moreover, while CCSD with a perturbative treatment of triple excitations, or CCSD(T), has claimed the status of “gold standard” among ground-state methods for achieving chemical accuracy in weakly correlated systems, it is not clear how a similar notion can be defined for finite-temperature methods. More benchmark studies need to be performed to understand and improve accuracy of thermal CC.


        This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Computational and Theoretical Chemistry Program under Award No. DE-FG02-09ER16053. G.E.S. acknowledges support as a Welch Foundation Chair (No. C-0036).






§ THERMAL CCSD EQUATIONS

When the cluster operator S is truncated to singles and double excitation operators only, i.e.,

    S = s_0 + ∑_pq s_pq a_p^†ã_q^† + 1/4∑_pqrs s_pqrs a_p^† a_q^†ã_s^†ã_r^†

the evolution equations (c.f. Eq. <ref>) for the CC amplitudes take the following form for the α-evolution,

    ∂ s_0/∂α   = 1/2ℛ_N^0 - S_α^0
        
    ∂ s_pq/∂α   = 1/2ℛ_N^pq - S_α^pq,
        
    ∂ s_pqrs/∂α   = 1/2ℛ_N^pqrs - S_α^pqrs,

where ℛ_N denotes the various CC residuals for the α-evolution,

    
    ℛ_N^0    = ⟨Ψ_0 |( e^-S N e^S  - N  ) |Ψ_0 ⟩,
            
        = ∑_a x_a y_a s_aa
    ℛ_N^pq   = ⟨Ψ_0 |ã_q a_p ( e^-S N e^S  -  N ) |Ψ_0 ⟩,
            
        = ( x_p^2 - y_q^2 ) s_pq
            + ∑_a x_a y_a ( s_apaq - s_pa s_aq)
            
    ℛ_N^pqrs   = ⟨Ψ_0 |ã_r ã_s a_q a_p ( e^-S N e^S  -  N ) |Ψ_0 ⟩,
            
        = ( x_p^2 + x_q^2 - y_r^2 - y_s^2 ) s_pqrs
         + 1/2𝒫(p q) 𝒫(r s)
                    ∑_a x_a y_a ( s_as s_pqar - s_pa s_aqrs)


and the operator-derivative terms are given by

    
    S_α^0    = -1/2∑_a x_a y_a s_aa
    
            S_α^pq   = 1/2∑_a x_a y_a ( s_pa s_aq + s_paaq)
            
    
            S_α^pqrs   = - 1/4𝒫(pq) 𝒫(rs) ∑_a x_a y_a ( s_as s_pqar - s_pa s_aqrs),


where x and y are the thermal Bogoliubov parameters and the dummy indices a, b, … are summed over all the spin-orbitals.
Equations for the β-evolution can be obtained in a similar way,

    ∂ s_0/∂β   = 1/2ℛ_H^0 - S_β^0
        
    ∂ s_pq/∂β   = 1/2ℛ_H^pq - S_β^pq,
        
    ∂ s_pqrs/∂β   = 1/2ℛ_H^pqrs - S_β^pqrs,

where ℛ_H denotes the various CC residuals for the β-evolution. Using the thermal Bogoliubov transformation described in Eq. <ref>, a general two-body Hamiltonian,

    H = ∑_p h_pq c_p^† c_p + 1/4∑_pqrs u_pqrs c_p^† c_q^† c_s c_r,

can be expressed in terms of thermal creation/annihilation operators, and takes the form


    
    H    = h_0 + ∑_ab [
                h^(11)_ab( a^†_a ã^†_b + h.c.)
                + h^(20)_ab a_a^† a_b + h^(02)_abã_a^†ã_b
            ]
            + ∑_abcd[
                h^(221)_abcd( a_a^† a_b^†ã_d^†ã_c^† + h.c.)
                + h^(222)_abcd a_a^†ã_b^†ã_d a_c
            .
            
               .
                + h^(31)_abcd( a_a^† a_b^†ã_c^† a_d + h.c.)
                + h^(13)_abcd( a_a^†ã_b^†ã_c^†ã_d + h.c.)
                + h^(40)_abcd a_a^† a_b^† a_d a_c + h^(04)_abcdã_a^†ã_b^†ã_d ã_c
            ]

    where we use h_0, h^(11), etc. to denote the effective matrix elements of the general quasiparticle Hamiltonian (h^(11) is associated with operators that contain a non-tilde and a tilde quasiparticle each, h^(20) with two non-tilde quasiparticle operators, and so on), which are given by
    
    h_0 = ∑_a y_a^2 h_aa + 1/2∑_ab y_a^2 y_b^2 u_abab

    
    h^(11)_ab = x_a y_b f_ab,    h^(20)_ab = x_a x_b f_ab,    h^(02)_ab = -y_a y_b f_ab,
              with  
            f_ab = δ_ab h_ab + ∑_c y_c^2 u_acbc

    
    h^(221)_abcd = 1/4 x_a x_b y_c y_d u_abcd,
              
            h^(222)_abcd   = x_a x_c y_b y_d u_adbc,
              
            h^(31)_abcd = - 1/2 x_a x_b y_c x_d u_abcd,
            
    
            h^(13)_abcd = -1/2 x_a y_b y_c y_d u_adbc,
              
            h^(40)_abcd   = 1/4 x_a x_b x_c x_d u_abcd,
              
            h^(04)_abcd = 1/4 y_a y_b y_c y_d u_abcd.

    In obtaining the above expressions, we have assumed real matrix elements in the Hamiltonian.
    The residuals can then be expressed compactly in terms of the effective Hamiltonian matrix elements,
    
        
    ℛ_H^0    = ⟨Ψ_0 |( e^-S H e^S  - H_0  ) |Ψ_0 ⟩,
                
        = h_0 - ∑_a y_a^2 _a + ∑_ab h^(11)_ab s_ab + ∑_abcd( 2 s_ac s_bd + s_abcd) h^(221)_abcd
    ℛ_H^pq   = ⟨Ψ_0 |ã_q a_p ( e^-S H e^S  -  H_0 ) |Ψ_0 ⟩,
                
       = h^(11)_pq - δ_pq_p x_p y_p + ∑_a ( h^(02)_aq s_pa + h^(20)_ap s_aq)
                - ∑_ab( h^(11)_ab ( s_aq s_pb + s_pabq ) - h^(222)_pqab s_ab )
                
          + ∑_abc( h^(13)_abcq( 2 s_ab s_pc + s_apbc) - h^(31)_abcp( 2 s_ac s_bq + s_abcq) )
                - 2 ∑_abcd h^(221)_abcd (
                    2 s_ac ( s_bq s_pd + s_pbdq ) - s_aq s_bpcd - s_pc s_abdq )
                
    ℛ_H^pqrs   = ⟨Ψ_0 |ã_r ã_s a_q a_p ( e^-S H e^S  -  H_0 ) |Ψ_0 ⟩,
                
       = 𝒫(pq) 𝒫(rs)  [
                    h^(221)_pqrs + ∑_a (
                        1/2( h^(02)_ar s_pqas + h^(20)_ap s_aqrs) + h^(13)_prsa s_qa + h^(31)_pqsas_ar)
                .
                
       .
                    + 1/2∑_ab(
                        h^(04)_abrs ( 2 s_pa s_qb + s_pqab ) + h^(40)_abpq ( 2 s_ar s_bs + s_abrs )
                        - h^(11)_ab ( s_ar s_pqbs + s_pb s_aqrs ) - 2 h^(222)_prab ( s_qabs + s_as s_qb )
                    )
                .
                
       .
                    + ∑_abc(
                        h^(13)_abcr(
                            s_ab s_pqcs + 1/2 s_as( 2 s_pb s_qc + s_pqbc )  - 2 s_pb s_aqcs)
                        + h^(31)_abcq(
                            s_ac s_bprs + 1/2 s_pc ( 2 s_ar s_bs + s_abrs ) - 2 s_ar s_bpcs)
                    )
                .
                
       .
                    - ∑_abcd h^(221)_abcd(
                        2 s_ac ( s_br s_pqds + s_pd s_bqrs )
                        - 1/2 s_ar s_bs ( 2 s_pc s_qd + s_pqcd ) + 4 s_ar s_pc s_bqds.
                .
                
       .
                    .
                        - 1/4 s_abrs ( 2 s_pc s_qd + s_pqcd ) + ( s_abcr s_pqds + s_apcd s_bqrs )
                        + 2 s_bqdr s_apcs)
                ],

    

where we have used the Baker-Campbell-Hausdorff expansion to simplify the similarity transformation, i.e.

    e^-S H e^S = H + [H,S] + 1/2! [[H,S],S] + ….

For a two-body Hamiltonian, with S truncated to at most double quasiparticle excitations, this expansion truncates at fourth order. Diagrammatic expressions for these equations can also be formulated along similar lines as Bogoliubov coupled cluster methods.<cit.>
The operator-derivative terms in the β-evolution are given by

    
    S_β^0    = 1/2∑_a _a x_a y_a s_aa
    
            S_β^pq   = -1/2∑_a _a x_a y_a ( s_pa s_aq + s_paaq)
            
    
            S_β^pqrs   = 1/4𝒫(pq) 𝒫(rs) ∑_a _a x_a y_a ( s_as s_pqar - s_qa s_aprs) .






§ FIXED-REFERENCE FORMULATION

In the fixed-reference formalism, which is quite analogous to the Schrödinger picture approach, we choose the state |𝕀⟩ as the zeroth order approximation to the thermal state. However, as we have discussed in Ref. , one can chose any value of chemical potential α_0 and temperature β_0 to construct the thermal reference state. It is merely a matter of comfort to use α_0 = 0 = β_0 as the corresponding initial conditions for the cluster amplitudes are trivial. For |𝕀⟩ as our choice of reference, it is convenient to redefine the thermal state as

    |ψ (α, β) ⟩ = e^α N-β H|𝕀⟩,

so that the thermal expectation value of any physical quantity A becomes

    
    ⟨ A ⟩   = ⟨𝕀|  A   e^-β H |𝕀⟩/⟨𝕀|  e^-β H |𝕀⟩, 
       = ⟨𝕀|  A  |ψ (α, β) ⟩/⟨𝕀|ψ (α, β) ⟩,


and a better bra is no longer required. Correspondingly, the governing imaginary time Schrödinger equations become

    
    
    ∂/∂β|ψ (α, β) ⟩   = -H |ψ (α, β) ⟩.
            
    ∂/∂α|ψ (α, β) ⟩   = N |ψ (α, β) ⟩.


The thermal state for a given chemical potential α at inverse temperature β can then be written as an exponential coupled cluster wavefunction

    |Ψ (α, β) ⟩ = e^T (α,β) |𝕀⟩,

where the cluster operator T (α,β) builds correlation atop |𝕀⟩. With this CC wavefunction ansatz, the finite-temperature expectation value of any physical quantity A becomes

    ⟨ A ⟩ = ⟨𝕀| A e^T|𝕀⟩ e^-t_0 = ⟨𝕀| e^-T A e^T|𝕀⟩.

The state |𝕀⟩ is annihilated by thermal quasiparticle operators a_p and ã_p corresponding to x_p = y_p = 1/√(2) in the Bogoliubov transformation in Eq. <ref>. We will refer these field operators as

    d_p,  d^†_p, d̃_p, d̃^†_p.

Therefore, the cluster operator T can be expressed as

    T = t_0
        + ∑_p,q t_pq d^†_p d̃^†_p
        + 1/(2!)^2∑_p,q,r,s t_pqrs d^†_p d^†_q d̃^†_s d̃^†_r
        + …,

where the α- and the β-dependence is carried by the cluster amplitudes.
These cluster amplitudes are found by integrating the imaginary time Schrödinger Eq. <ref>, which, upon substituting the wavefunction ansatz of Eq. <ref>, gives the following working equation

    
    
    ∂ T/∂α |𝕀⟩   =
            e^-T N e^T |𝕀⟩,
            
    ∂ T/∂β |𝕀⟩   =
            - e^-T H e^T |𝕀⟩.


Like conventional ground-state CC, Eq. <ref> can be left-projected with the ground and excited slater determinants to yield the evolution equations for the amplitudes,

    
    
    ∂ t_0/∂β   = - 1/𝒵_𝕀⟨𝕀|  e^-T H e^T |𝕀⟩,
    ∂ t_pq/∂β   = - 1/𝒵_𝕀⟨𝕀| d̃_q d_p |  e^-T H e^T |𝕀⟩,
    ∂ t_pqrs/∂β   = - 1/𝒵_𝕀⟨𝕀| d̃_r d̃_s d_q d_p   e^-T H e^T |𝕀⟩,


and so on, where 𝒵_𝕀 = ⟨𝕀|𝕀⟩. Similar equations can be derived for evolution along α. These equations can be integrated starting from β = 0 (or any other value), where the initial values of the amplitudes are known, to the required inverse temperature and chemical potential. Here, since |𝕀⟩ is exact at α, β = 0, we have the initial conditions,

    t_0 = 0,   t_pq = 0,   t_pqrs = 0,  ….




