Thermofield theory for finite-temperature coupled cluster

TL;DR

This paper develops a thermofield coupled cluster approach to accurately compute properties of many-electron systems at finite temperatures, extending previous methods with a thermal quasiparticle framework and benchmarking on model and real systems.

Contribution

It introduces a thermofield-based coupled cluster method with an exponential ansatz for thermal states, enabling finite-temperature electronic structure calculations.

Findings

Accurately computed properties of model and real systems at finite temperature.

Validated the method against exact results for benchmark systems.

Extended the thermofield framework to coupled cluster theory.

Abstract

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.

Equations

$$⟨A{⟩}_{\beta }=\frac{1}{\mathcal{Z}}\mathrm{Tr}(A\rho ),$$

$$\rho =e^{-\beta H^{\prime }},\mathcal{Z}=\mathrm{Tr}(\rho ),$$

$$|\Psi ⟩=e^T|\Phi ⟩,$$

$$T$$

$$T_1$$

$$T_2$$

$$\bar{H}|\Phi ⟩=E|\Phi ⟩,\bar{H}=e^{-T}H{e}^T$$

$$E$$

$$0$$

$$E=⟨\Phi |(1+Z)e^{-T}H{e}^T|\Phi ⟩$$

$$Z$$

$$Z_1$$

$$Z_2$$

$$⟨{\Psi }^{\prime }|=⟨\Phi |(1+Z)e^{-T}=⟨\Phi |(1+W)e^{w_0},$$

$$⟨A⟩=\frac{⟨\Psi (\alpha ,\beta )|A|\Psi (\alpha ,\beta )⟩}{⟨\Psi (\alpha ,\beta )|\Psi (\alpha ,\beta )⟩}.$$

$$|\Psi (\alpha ,\beta )⟩=e^{(\alpha N-\beta H)/2}|\mathbb{I}⟩,|\mathbb{I}⟩=\sum _m|m,\tilde{m}⟩,$$

$$\frac{\partial }{\partial \beta }|\Psi (\alpha ,\beta )⟩$$

$$\frac{\partial }{\partial \alpha }|\Psi (\alpha ,\beta )⟩$$

$$|\mathbb{I}⟩=\prod _{p\in \mathrm{levels}}\left(|0,\tilde{0}{⟩}_p+|1,\tilde{1}{⟩}_p\right),$$

$$H_0=\sum _p{ϵ}_p{c}_p^{†}{c}_p$$

$$|0(\alpha ,\beta )⟩=\prod _{p\in \mathrm{levels}}\left(x_p|0,\tilde{0}{⟩}_p+y_p|1,\tilde{1}⟩\right),$$

$$x_p$$

$$y_p$$

$$\left[\begin{array}{c}{a}_p\\ {\tilde{a}}_p^{†}\end{array}\right]=\left[\begin{array}{cc}{x}_p& -y_p\\ y_p& x_p\end{array}\right]\left[\begin{array}{c}{c}_p\\ {\tilde{c}}_p^{†}\end{array}\right],$$

$$a_p|0(\alpha ,\beta )⟩=0={\tilde{a}}_p|0(\alpha ,\beta )⟩.$$

$$⟨A{⟩}_{\mathrm{CC}}=\frac{⟨{\Psi }^{\prime }|A|\Psi ⟩}{⟨{\Psi }^{\prime }|\Psi ⟩},$$

$$|\Psi ⟩$$

$$S$$

$$⟨{\Psi }^{\prime }|=⟨{\Psi }_0|(1+Z)e^{z_0}{e}^{-S},$$

$$⟨{\Psi }^{\prime }|$$

$$W$$

$$⟨{\Psi }_0|\left(\frac{\partial W}{\partial \alpha }+(1+W)\frac{\partial w_0}{\partial \alpha }\right)$$

$$⟨{\Psi }_0|\left(\frac{\partial W}{\partial \beta }+(1+W)\frac{\partial w_0}{\partial \beta }\right)$$

$$H_{\text{CI}}$$

$$N_{\text{CI}}$$

Full text

I 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 $\mathrm{eV}$) , which is much larger than room temperature ($\sim 25\mathrm{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,Dirac (1929) 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 theoryCrawford and Schaefer (2000); Bartlett and Musiał (2007) (CC), Monte Carlo methods,Hammond et al. (1994); Foulkes et al. (2001); Zhang and Krakauer (2003); Zhang (2004); Al-Saidi et al. (2006) among many more.White (1993); Östlund and Rommer (1995); Georges et al. (1996)

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,Lee et al. (2006) ultra-cold chemistry,Balakrishnan (2016); Bohn et al. (2017) geochemical processes which generally involve very high temperatures and pressure,Guillot (1999) 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 $\beta$, the ensemble of choice is generally canonical or grand-canonical and the expectation value of an observable $A$ is defined by

[TABLE]

where $\rho$ is the thermal density matrix, and $\mathcal{Z}$ is the partition function, given by

[TABLE]

with $H^{\prime}=H$ for the canonical ensemble and $H^{\prime}=H-\mu N$ for the grand canonical ensemble, where $\mu$ is the chemical potential. In this paper, we shall work explicitly with the grand canonical ensemble.

Exactly computing $\rho$ (or equivalently $\mathcal{Z}$) 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.Mermin (1963); Sokoloff (1967); Lichtenstein et al. (2001); Verstraete et al. (2004); Feiguin and White (2005); White (2009); Stoudenmire and White (2010); Pittalis et al. (2011); Pižorn et al. (2014); Hermes and Hirata (2015); Czarnik et al. (2016); Santra and Schirmer (2017); Zgid and Gull (2017); Claes and Clark (2017); White and Chan (2018); Hummel (2018) Most of these methods can be broadly categorized into deterministic methods such as diagrammatic perturbation theory based on the Matsubara formalism,Matsubara (1955) 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 groupVerstraete et al. (2004); Feiguin and White (2005) and finite-temperature perturbation theories,He et al. (2014); Santra and Schirmer (2017) 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 methodSanyal et al. (1992, 1993); Mandal et al. (1998, 2003) was proposed by Mukherjee et. al. and has been further elaborated recently in independent works by White *et al.*White and Chan (2018) and Hummel.Hummel (2018) 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),Matsumoto et al. (1983); Semenoff and Umezawa (1983); Umezawa (1984); Evans et al. (1992) using a framework we recently explored in Ref. 42. 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.Suzuki (1985); Hatsuda (1989); Walet and Klein (1990); de Vega and Bañuls (2015); Borrelli and Gelin (2016); Nocera and Alvarez (2016); Chen and Zhao (2017); Borrelli and Gelin (2017); Wu and Hsieh (2018) A connection between thermofield formalism and thermal cluster-cumulant theory was also explored by Mukherjee et. al. in Ref. 35. Here, we parametrize this so-called thermal wavefunction as an exponential ansatz to integrate the evolution equation.

II 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ümmelCoester (1958); Coester and Kümmel (1960) and later reformulated for electronic structure theory by Čižek and Paldus,Čižek and Paldus (1971) it uses an exponential wavefunction ansatz

[TABLE]

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

[TABLE]

where we have followed the standard notation for labeling orbital indices, i.e., occupied orbitals are denoted by indices $i,j,k,\ldots$, while unoccupied orbitals are denoted by $a,b,c,\ldots$. Since the wavefunction $|\Psi\rangle$ is constructed as excitations on a single Slater determinant $|\Phi\rangle$, 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

[TABLE]

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

[TABLE]

where $\langle\Phi_{i}^{a}|$ and $\langle\Phi_{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,Arponen (1983); Helgaker and Jørgensen (1988) i.e. one makes the energy functional

[TABLE]

stationary with respect to $T$ and $Z$, where

[TABLE]

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

[TABLE]

where $w_{0}$ is a constant and $W$ has the same operator-form as $Z$.

A similar formulation known as equation of motion CCRowe (1968) (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 ($\mathcal{O}(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.

III 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 $|\Psi(\alpha,\beta)\rangle$, with $\alpha=\beta\mu$, 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. 1 can be replaced by an expectation value over this wavefunction, i.e.

[TABLE]

That the thermal state $|\Psi(\alpha,\beta)\rangle$ cannot be a pure state in the physical Hilbert space $\mathcal{H}$ is easily established since the density matrix $\rho$ 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 $\tilde{\mathcal{H}}$, and the thermal state is then defined in the doubled space $\mathcal{H}\otimes\tilde{\mathcal{H}}$ as

[TABLE]

where $H$ is the Hamiltonian, $\{|m\rangle\}$ is some orthonormal basis in the Hilbert space $\mathcal{H}$ and $|\tilde{m}\rangle$ is the tilde-state corresponding to $|m\rangle$, while we have used the shorthand notation $|m,\tilde{m}\rangle=|m\rangle\otimes|\tilde{m}\rangle$. The state $|\mathbb{I}\rangle$ is the state with maximal entanglement between $\mathcal{H}$ and $\tilde{\mathcal{H}}$ and is the exact thermal state at infinite temperature (i.e. $\beta=0$) and $\alpha=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 $\tilde{\mathcal{H}}$ can be found in Ref. Harsha et al., 2019 and the references therein.

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

[TABLE]

One can solve for $|\Psi(\alpha,\beta)\rangle$ by integrating Eq. 14, generally starting from $(\alpha,\beta)=(0,0)$, where the initial thermal state is known exactly. Exactly evolving the thermal state $|\Psi(\alpha,\beta)\rangle$ is equivalent to the exact diagonalization of the density matrix $\rho$. 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 $|\mathbb{I}\rangle$ can be expressed in terms of single-particle Fock states,

[TABLE]

where $|0\rangle_{p}$ and $|1\rangle_{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.

[TABLE]

the normalized mean-field thermal state can be written as

[TABLE]

where $x_{p}$ and $y_{p}$ are related to the Fermi-Dirac statistics,

[TABLE]

with $x_{p}^{2}+y_{p}^{2}=1$. The mean-field thermal state in Eq. 16 allows us to introduce a thermal Bogoliubov transformation,

[TABLE]

such that

[TABLE]

Correlated methods can be built with either $|\mathbb{I}\rangle$ or $|0(\alpha,\beta)\rangle$ as the reference, while integrating Eq. 14. As we have discussed in Ref 42, the former choice, which we call the fixed-reference formalism, performs well only in the vicinity of $\beta=0$, while the latter, called the covariant formalism, yields accurate results for the entire range of $\beta$. In the next section, we present details to integrate Eq. 14 through the CC ansatz in the covariant formalism. We have also included a short discussion on the fixed-reference approach in Appendix B. For brevity of notation, we will henceforth refer the mean-field thermal state $|0(\alpha,\beta)\rangle$ by $|\Psi_{0}\rangle$.

IV Thermal Coupled Cluster

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

[TABLE]

where both the ket $|\Psi\rangle$ and the bra $\langle\Psi^{\prime}|$ states are approximations to the same thermal state, and consequently evolve according to Eq. 14 and its adjoint respectively.

Given that the Bogoliubov transformation in Eq. 18 is BCS-like, i.e. orbitals $p$ (in $\mathcal{H}$) and $\tilde{p}$ (in $\mathcal{\tilde{H}}$) in the mean-field thermal state (Eq. 16) 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 operatorsHenderson et al. (2014) acting on an $\alpha$- and $\beta$-dependent mean-field thermal reference, $|\Psi_{0}\rangle$, defined in Eq. 16,

[TABLE]

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

[TABLE]

and as explained in Eq. 11, it can be expressed as an effective CI wavefunction

[TABLE]

For both the bra and the ket states, the reference $|\Psi_{0}\rangle$ evolves continuously as we evolve the Schrödinger Eq. 14 (hence the name ‘covariant’). Accordingly, both the amplitudes ($s_{pq}$, $w_{pq}$, etc.) and the quasiparticle operators $a^{\dagger},\tilde{a}^{\dagger}$ carry $\alpha$- and $\beta$-dependence. The $\alpha$- and $\beta$-evolution of $\langle\Psi^{\prime}|$, a thermal CI wavefunction, is governed by

[TABLE]

where $H_{\textrm{CI}}$ and $N_{\mathrm{CI}}$ are effective CI Hamiltonian and Number operators respectively, and are given by

[TABLE]

A detailed discussion on thermal CI and the derivation of these equations can be found in Ref. Harsha et al., 2019.

For the evolution of $|\Psi\rangle$, substituting the CC ansatz from Eq. 20a into the Schrödinger Eq. 14 gives

[TABLE]

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

IV.1 Wilcox identity

The process of reducing Eq. 24 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.

[TABLE]

where $x=\alpha,\beta$. The derivative of the exponential cluster operator is appropriately performed by making use of the Wilcox identity,Wilcox (1967) which states that the derivative of the exponential of an operator $M$ with respect to some parameter $\lambda$ can be evaluated as

[TABLE]

With this, the left-hand side of Eq. 24 becomes

[TABLE]

where we have used the shorthand $\partial_{x}$ for $\partial/\partial x$, and made use of the Baker-Campbell-Hausdorff expansion in going from Eq. 26a to 26b. Finally, breaking the derivative $\partial_{x}S$ into the amplitude ($\partial_{\mathrm{amp}}S$) and operator ($\partial_{\mathrm{op}}S$) derivatives,

[TABLE]

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

[TABLE]

where the integration over $y$ in the second equation can be carried out analytically, as explained in Eq. 26b. With these details, Eq. 24 can be further simplified as

[TABLE]

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 $\alpha$-$\beta$ space. Complete expressions for the CCSD evolution equations are included in Appendix A.

V Implementation

V.1 Integration

An explicit Runge-Kutta method of order (4)5 with step-size controlDormand and Prince (1980) 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.

V.2 Choice of $H_{0}$

As mentioned earlier, for the covariant thermal CC, a mean-field Hamiltonian $H_{0}=\sum_{p}\epsilon_{p}a_{p}^{\dagger}a_{p}$ is used to construct the mean-field thermal reference state $|\Psi_{0}\rangle$. There are several alternatives to pick this mean-field Hamiltonian – in thermal HF theory,Mermin (1963); Sokoloff (1967) one uses an $H_{0}$ that optimizes the applicable free energy, and it therefore depends on $\alpha$ and $\beta$. 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 $\alpha$- and $\beta$-independent $H_{0}$ is analogous to the imaginary-time interaction picture formalism and simplifies our implementation since, (i) an $\alpha$- and $\beta$-independent $H_{0}$ leads to clean analytical forms for the mean-field evolution operator as well as the $\alpha$- and $\beta$-derivatives of the thermal quasiparticle operators, and (ii) $H_{0}$ being diagonal yields a straightforward thermal Bogoliubov transformation in Eq 18. We note, however, that this choice is different from the one used in Ref. Harsha et al., 2019, where we use simply the one-electron Hamiltonian to construct $H_{0}$.

VI 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,Hubbard J. and Flowers Brian Hilton (1963) the pairing or the reduced BCS model, as well as chemical systems (atomic Beryllium and molecular $\mathrm{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 $\beta$ grid-point, we evolve the thermal states in $\alpha$ to fix the average number of particles before evolving again in $\beta$. 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 et al., 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

[TABLE]

where $\langle,\rangle$ 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,\sigma}=c_{p,\sigma}^{\dagger}\,c_{p,\sigma}$ is the number operator for lattice site $p$ and spin $\sigma$. The ratio $U/t$ characterizes the correlation strength.

Figure 1a 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 1b 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 $\beta\rightarrow\infty$ (or $\theta\rightarrow 0$, where $\theta=1/\beta$ 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 $\beta=0$. Such a behaviour can be anticipated since this method uses the thermal reference corresponding to $\beta=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

[TABLE]

where $N_{p}$ counts the number of electrons, $\epsilon_{p}$ denotes the energy, and $P^{\dagger}_{p}$/$P_{p}$ respectively creates/annihilates a pair of electrons in the $p^{\mathrm{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., $\Delta\epsilon=\epsilon_{p+1}-\epsilon_{p}=1$. Figure 2a 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 2b 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 3 and 4 show similar trends for atomic Beryllium and molecular H2 at bond length $0.74\textup{\AA}$ 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 5 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 6, we show the $z$-component spin-spin correlation function

[TABLE]

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 $\theta$, there is no correlation between adjacent spins. As $\theta$ 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 $\beta\rightarrow\infty$.

VII 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 $\mathcal{O}(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.

Acknowledgements.

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

Appendix A Thermal CCSD Equations

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

[TABLE]

the evolution equations (c.f. Eq. 24) for the CC amplitudes take the following form for the $\alpha$-evolution,

[TABLE]

where $\mathcal{R}_{N}$ denotes the various CC residuals for the $\alpha$-evolution,

[TABLE]

and the operator-derivative terms are given by

[TABLE]

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

[TABLE]

where $\mathcal{R}_{H}$ denotes the various CC residuals for the $\beta$-evolution. Using the thermal Bogoliubov transformation described in Eq. 18, a general two-body Hamiltonian,

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

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,

[TABLE]

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

[TABLE]

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.Signoracci et al. (2015) The operator-derivative terms in the $\beta$-evolution are given by

[TABLE]

Appendix B Fixed-reference formulation

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

[TABLE]

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

[TABLE]

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

[TABLE]

The thermal state for a given chemical potential $\alpha$ at inverse temperature $\beta$ can then be written as an exponential coupled cluster wavefunction

[TABLE]

where the cluster operator $T(\alpha,\beta)$ builds correlation atop $|\mathbb{I}\rangle$. With this CC wavefunction ansatz, the finite-temperature expectation value of any physical quantity $A$ becomes

[TABLE]

The state $|\mathbb{I}\rangle$ is annihilated by thermal quasiparticle operators $a_{p}$ and $\tilde{a}_{p}$ corresponding to $x_{p}=y_{p}=1/\sqrt{2}$ in the Bogoliubov transformation in Eq. 18. We will refer these field operators as

[TABLE]

Therefore, the cluster operator $T$ can be expressed as

[TABLE]

where the $\alpha$- and the $\beta$-dependence is carried by the cluster amplitudes. These cluster amplitudes are found by integrating the imaginary time Schrödinger Eq. 49, which, upon substituting the wavefunction ansatz of Eq. 52, gives the following working equation

[TABLE]

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

[TABLE]

and so on, where $\mathcal{Z}_{\mathbb{I}}=\langle\mathbb{I}|\mathbb{I}\rangle$. Similar equations can be derived for evolution along $\alpha$. These equations can be integrated starting from $\beta=0$ (or any other value), where the initial values of the amplitudes are known, to the required inverse temperature and chemical potential. Here, since $|\mathbb{I}\rangle$ is exact at $\alpha,\beta=0$, we have the initial conditions,

[TABLE]