Thermofield theory for finite-temperature electronic structure

Gaurav Harsha Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48105 Thomas M. Henderson Department of Physics and Astronomy, Rice University, Houston, Texas 77005 Department of Chemistry, Rice University, Houston, Texas 77005 Gustavo E. Scuseria Department of Physics and Astronomy, Rice University, Houston, Texas 77005 Department of Chemistry, Rice University, Houston, Texas 77005

Abstract

Wave-function methods have offered a robust, systematically improvable means to study ground-state properties in quantum many-body systems. Theories like coupled cluster and their derivatives provide highly accurate approximations to the energy landscape at a reasonable computational cost. Analogs of such methods to study thermal properties, though highly desirable, have been lacking because evaluating thermal properties involve a trace over the entire Hilbert space, which is a formidable task. Besides, excited-state theories are generally not as well studied as ground-state ones. In this mini-review, we present an overview of a finite-temperature wave function formalism based on thermofield dynamics to overcome these difficulties. Thermofield dynamics allows us to map the equilibrium thermal density matrix to a pure state, i.e., a single wave function, albeit in an expanded Hilbert space. Ensemble averages become expectation values over this so-called thermal state. Around this thermal state, we have developed a procedure to generalize ground-state wave function theories to finite temperatures. As explicit examples, we highlight formulations of mean-field, configuration interaction, and coupled cluster theories for thermal properties of fermions in the grand-canonical ensemble. To assess the quality of these approximations, we also show benchmark studies for the one-dimensional Hubbard model, while comparing against exact results. We will see that the thermal methods perform similarly to their ground-state counterparts, while merely adding a pre-factor to the asymptotic computational cost. They also inherit all the properties, good or bad, from the ground-state methods, signifying the robustness of our formalism and the scope for future development.

I Introduction

The design and study of new materials and molecules is becoming heavily reliant on the ability to describe and predict their macroscopic properties from a quantum mechanical treatment of the underlying electronic structure. Many conventional applications involve molecules and materials with optical gaps much larger than the temperature of the environment, e.g., large gap molecules and insulating solids at room temperatures. The properties in such systems are governed solely by the ground state because the thermal energy is incapable of exciting enough electrons into higher energy states to make any practically observable difference. On the other hand, for quantum systems with optical gaps comparable to, or smaller than, the thermal energy, the ground sate is no longer a sufficient description. Accounting for thermal effects become essential for an accurate understanding of the electronic structure under these conditions, which often arise in novel applications based on strongly correlated materials (e.g., high-Tc superconductors and transition metal complexes at room temperature). [1, 2] Other examples include, but are not limited to, chemical reactions driven by hot-electrons [3] and those or in extreme geological environments. [4, 5]

Finite-temperature properties are defined as an average of all expected outcomes of an observable, weighted over an appropriate ensemble. For example, in the grand-canonical ensemble, the expectation value of an observable $A$ at inverse temperature $\beta=1/k_{B}T$ and chemical potential $\mu$ is defined as

\braket{A}(\beta,\mu)=\frac{1}{\mathcal{Z}}\sum_{m}\braket{m}{A}{m}e^{-\beta(E_{m}-\mu N_{m})},

(1)

where $\bra{m}$ and $\ket{m}$ are the Fock-space bra- and ket-eigenvectors of the system Hamiltonian $H$ with energy $E_{m}$ and particle number $N_{m}$ , and $\mathcal{Z}=\sum_{m}e^{-\beta(E_{m}-\mu N_{m})}$ is the partition function. This expectation value can also be expressed as a trace of $A$ over the ensemble density matrix

\braket{A}(\beta,\mu)=\mathrm{Tr}(\rho A),

(2)

where $\rho=e^{-\beta(H-\mu N)}/\mathcal{Z}$ is the Gibbs-Boltzmann density operator, with $H$ and $N$ being the Hamiltonian and the total particle-number operator. A direct application of conventional, state-specific, wave function methods to calculate thermal ensemble averages would be highly impractical because:

1.

The number of states in the Hilbert space that need to be counted in the ensemble grows exponentially with the system size.
2.

Approximating individual eigenstates of the Hamiltonian is a challenging task. While several reliable methods are now available for the ground-state, few excited state methods are equally accurate and efficient.

As a result, the many-body problem becomes far more difficult at finite temperatures. Therefore, beyond mean-field methods, [6, 7] Matsubara perturbation theory [8] (generally employed in Green’s function methods) as well as quantum Monte Carlo have been considered an important ingredient in a successful finite-temperature electronic structure method.

On the other hand, state-specific wave function theories offer more control, are systematically improvable, and often provide accurate results at a computational cost that grows only as a polynomial in system size. For instance, Hartree-Fock (HF) theory is a mean-field method that generally provides a good qualitative result at $\mathcal{O}(n^{3})$ computational scaling ( $n$ being a measure of the system size). Configuration interaction (CI) with singles and doubles (CISD) and second-order Møller-Plesset perturbation theory (MP2) improve upon HF but require increased computational effort ( $\mathcal{O}(n^{6})$ for CISD and $\mathcal{O}(n^{5})$ for MP2). Coupled cluster (CC) theory, [9, 10] truncated at singles and doubles as well (CCSD), generally outperforms HF, CISD, and MP2 and scales as $\mathcal{O}(n^{6})$ . Wave function methods have accumulated decades of experience, and several existing codes and numerical techniques are available to make them even more efficient. [11, 12] Because of such important qualities, the development of new wave function methods to study quantum mechanical systems at finite temperatures is highly desirable. Consequently, finite-temperature wave function methods have seen a growing interest in recent years. Beyond the thermal mean-field theory, several new methods have been proposed, e.g. Ancilla density matrix renormalization group (DMRG), [13, 14, 15, 16] minimally entangled typical thermal states, [17, 18] thermal density functional theory, [19] perturbation theory methods, [20, 21] and CI- or CC-like techniques, [22, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32] to name a few.

In this conference article, we present a short review of thermal wave function methods, with particular emphasis on a new formalism for generalizing existing ground-state wave function theories to finite-temperatures by using the theory of thermofield dynamics. [33, 34, 35] Most of the theory, as well as some key results, discussed here were reported by the authors in Refs. 28, 29, 31, 32.

II Wave function methods at finite temperature

While finite-temperature (or imaginary-time) generalizations of the time-dependent perturbation theory have been around for over half a century, one of the earliest ideas to generalize a wave function theory (coupled cluster, in this context) to finite-temperatures was proposed by Sanyal, Mukherjee and co-authors. In a series of papers, [22, 22, 23, 24] they developed the thermal cluster cumulant (TCC) theory, in which the interaction picture Gibbs-Boltzmann operator is parametrized using a CC-like exponential operator. By using a thermal generalization of the Wick’s theorem and associated normal ordering for this exponential ansatz, they showed that the free energy can be expressed using fully contracted term in the normal ordered expansion.

Recently, different ways to generalize the coupled cluster theory to finite temperatures have been sought separately by Hummel, [27] and White et al. [26] Both of these works rely on re-interpreting the coupled cluster diagrams from the perspective of time-dependent perturbation theory on the imaginary-time axis. While Hummel and White’s approaches differ in inspiration, by using the time-dependent perturbation theory, they end up relying on the thermal Wick’s theorem, much like in Sanyal and Mukherjee’s TCC theory. Consequently, up to difference in implementation, Hummel, Sanyal, and White’s finite-temperature coupled cluster result in similar working equations.

Another formalism for thermal wave function methods, inspired by thermofield dynamics, was proposed by the authors of this review. The thermofield approach relies on an explicit wave function representation of the Gibbs-Boltzmann density operator, which allows a direct generalization of any ground-state wave function ansatzes to finite temperatures.

Using thermofield dynamics or related principles to construct thermal wave function theories is not new. A connection of the TCC with the thermofield dynamics was explored by Sanyal and Mukherjee in Refs. 22, 23. They later noted that the thermofield theory was not a necessary ingredient in TCC. Ancilla DMRG also builds on techniques similar to thermofield dynamics for purification of the ensemble density matrix. Furthermore, a coupled cluster formalism for finite-temperature dynamics, inspired by thermofield dynamics, was also explored by Shushkov et al. in Ref. 30. However, to our knowledge, a general framework was proposed only in Ref. 28.

Refer to caption — Figure 1: Physical interpretation of how a single wave function in an enlarged space in the thermofield framework can encode the thermal ensemble density matrix

III Thermofield Dynamics

Thermofield dynamics was originally proposed to study time-dependent phenomena at finite-temperature using path integrals, much like the Keldysh approach [36] in the Matsubara imaginary time formalism. It provides a way to construct a pure-state representation for the thermal density matrix. That is, we can write down a thermal wave function $\ket{\Psi(\beta,\mu)}$ , such that the trace over an ensemble of states in Eq. 1 can be expressed as an expectation value over the thermal state, i.e.,

\braket{A}=\frac{\braket{\Psi(\beta,\mu)}{A}{\Psi(\beta,\mu)}}{\langle\Psi(\beta,\mu)|\Psi(\beta,\mu)\rangle}.

(3)

Such a purification of the ensemble density matrix is achieved by working in an enlarged space that is made up from

•

the original Hilbert space $\mathcal{H}$ ,
•

and its complex conjugate copy $\tilde{\mathcal{H}}$ , also known as the tilde-conjugate space.

The thermal state is defined as

\ket{\Psi(\beta,\mu)}=\frac{1}{\mathcal{Z}}\sum_{m}e^{-\beta(E_{m}-\mu N_{m})/2}\ket{m}\ket{\tilde{m}},

(4)

where, once again, $\beta=1/k_{B}T$ is the inverse temperature, $\mathcal{Z}$ is the partition function, $\ket{m}$ is an eigenstate of the Hamiltonian $H$ with eigenvalue $E_{m}$ , and $N_{m}$ number of particles, and $\ket{\tilde{m}}$ is the conjugate-space counterpart of $\ket{m}$ . In the grand-canonical ensemble, the thermal state can be explicitly written as


$\displaystyle\ket{\Psi(\beta,\mu)}$	$\displaystyle=e^{-\beta(H-\mu N)/2}\ket{\mathbb{I}},$	(5a)
$\displaystyle\ket{\mathbb{I}}$	$\displaystyle=\prod_{p}\left(1+c_{p}^{\dagger}\tilde{c}_{p}^{\dagger}\right)\ket{-;-},$	(5b)

where $p$ is a spin-orbital index, $c_{p}^{\dagger}$ creates a particle in the $p$ th physical spin-orbital and $\tilde{c}^{\dagger}_{p}$ creates a particle in the corresponding conjugate spin-orbital. Finally, $\ket{-;-}$ denotes the vacuum for physical and conjugate spaces. The state $\ket{\mathbb{I}}$ is also 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 $\mu=0$ .

A simple physical interpretation for the thermofield approach is that by introducing the tilde-conjugate degrees of freedom, the system and the thermal bath (or the environment) are treated equally and as one big, isolated system. The thermal state can be thought of as the wave function for this combined macro system. Since the dimensionality of a typical quantum system is much smaller than that of the environment, we can use the Schmidt decomposition [37] to factorize the thermal state in the form of Eq. 5. Figure 1 provides a graphical depiction of this interpretation. Thermofield dynamics has already shown great potential in the study many-electron systems in chemical and condensed-matter systems.[38, 39, 40, 41, 42, 43, 44, 45, 46]

IV Finite-temperature wave function methods

By definition, the grand-canonical thermal state in Eq. 5 satisfies the following evolution equations in the inverse temperature $\beta$ and chemical potential $\mu$ :


$\displaystyle\frac{\partial\ket{\Psi(\beta,\mu)}}{\partial\beta}$	$\displaystyle=-\frac{1}{2}H\ket{\Psi(\beta,\mu)},$	(6a)
$\displaystyle\frac{\partial\ket{\Psi(\beta,\mu)}}{\partial\mu}$	$\displaystyle=\frac{\beta}{2}N\ket{\Psi(\beta,\mu)}.$	(6b)

An exact solution of these evolution equations is equivalent to exact diagonalization (or Full CI) in the Fock space. However, since we only need to work with a single wave function, we can use any conventional ground-state wave function approximation to solve and integrate the imaginary time ( $\beta$ ) and chemical potential ( $\mu$ ) evolution equations. Here, we will show explicit constructions for the mean-field, CI and CC approximations.

IV.1 Mean-field theory

In the mean-field approximation, we use an effective one-body Hamiltonian $H_{0}$ to integrate the evolution equations in Eq. 6. As long as $H_{0}$ does not depend on $\beta$ and $\mu$ , the integration can be performed analytically and the mean-field thermal state $\ket{\Phi(\beta,\mu)}$ is given by

\ket{\Phi(\beta,\mu)}=e^{-\beta(H_{0}-\mu N)/2}\ket{\mathbb{I}}.

(7)

For electronic systems, a convenient choice for $H_{0}$ is the sum of Fock operators from the ground-state Hartree-Fock state. In the Hartree-Fock molecular orbital basis, in which the Fock matrix is diagonal and $H_{0}=\sum_{p}\epsilon_{p}c_{p}^{\dagger}c_{p}$ , the normalized mean-field thermal state can be factorized as

\ket{\Phi(\beta,\mu)}=\prod_{p}\left(u_{p}+v_{p}c_{p}^{\dagger}\tilde{c}_{p}^{\dagger}\right)\ket{-;-},

(8)

where $u_{p}=1/\sqrt{1+e^{-\beta(\epsilon_{p}-\mu)}}$ , and $v_{p}=\sqrt{1-u_{p}^{2}}$ . Here, the product runs over all spin-orbitals $p$ . In order to maintain charge neutrality, the chemical potential $\mu$ is adjusted for every $\beta$ grid-point. The mean-field thermal state in Eq. 8 is a Bardeen-Cooper-Schrieffer (BCS) state which reduces to the ground-state Hartree-Fock in the zero-temperature limit where $v_{\text{occupied}}=1=u_{\text{virtual}}$ , and $v_{\text{virtual}}=0=u_{\text{occupied}}$ . The thermal analogue of Cooper pairs are formed between the physical and auxiliary orbital pairs $\{p,\tilde{p}\}$ .

IV.2 Correlated methods

For an accurate description of quantum mechanical systems, we need to go beyond mean-field theory and capture the effects of correlation in the system. Traditional ground-state methods such as perturbation theory, CI, and CC are some of the most popular and effective methods for this task. In general, a correlated approximation to the thermal state is constructed as

\ket{\Psi(\beta,\mu)}=\Gamma(\beta,\mu)\ket{\Phi(\beta,\mu)},

(9)

where the wave operator $\Gamma(\beta,\mu)$ adds correlation atop the mean-field reference. The wave operator $\Gamma(\beta,\mu)$ is constructed using quasiparticle excitation operators for the mean-field reference $\ket{\Phi(\beta,\mu)}$ , in the same way as particle-hole excitations for a Hartree-Fock reference are used in ground-state theories. These thermal quasiparticles operators, $\{a_{p},a_{p}^{\dagger},\tilde{a}_{p},\tilde{a}_{p}^{\dagger}\}$ are defined in such a way that

a_{p}(\beta,\mu)\ket{\Phi(\beta,\mu)}=0=\tilde{a}_{p}(\beta,\mu)\ket{\Phi(\beta,\mu)}.

(10)

For the grand canonical mean-field state, defined in Eq. 8, the fermion and quasiparticle creation and annihilation operators are related via a Bogoliubov transformation:

\begin{bmatrix}a_{p}\\ \tilde{a}_{p}^{\dagger}\end{bmatrix}=\begin{bmatrix}u_{p}&-v_{p}\\ v_{p}&u_{p}\end{bmatrix}\begin{bmatrix}c_{p}\\ \tilde{c}_{p}^{\dagger}\end{bmatrix}.

(11)

The $\beta$ - and $\mu$ -evolution equations for the correlated wave function becomes


$\displaystyle\frac{\partial\Gamma}{\partial\beta}\ket{\Phi}$	$\displaystyle=-\frac{1}{2}(H\Gamma-\Gamma H_{0})\ket{\Phi},$	(12a)
$\displaystyle\frac{\partial\Gamma}{\partial\mu}\ket{\Phi}$	$\displaystyle=\frac{\beta}{2}(N\Gamma-\Gamma N)\ket{\Phi},$	(12b)

where both $\Gamma$ and $\ket{\Phi}$ depend on $\beta$ and $\mu$ , as noted in Eq. 9. Comparing different orders of quasiparticle excitations on the left and right hand side of this equation defines our initial value problem. Note that at $\beta=0$ , mean-field theory is exact, and $\Gamma(\beta=0,\mu)=1$ .

IV.2.1 Configuration interaction theory

For the CI approximation, the wave operator $\Gamma(\beta,\mu)$ is expressed as a linear combination of various quasiparticle excitations. At the level of single and double excitations (CISD), the thermal state is defined as

\ket{\Psi_{CI}(\beta,\mu)}=e^{t_{0}}(1+T)\ket{\Phi(\beta,\mu)},

(13)

where $t_{0}$ is a scalar that keeps track of the norm, and the excitation operator $T$ is defined as

	$\displaystyle T$	$\displaystyle=\sum_{pq}t_{pq}a_{p}^{\dagger}(\beta,\mu)\tilde{a}_{q}^{\dagger}(\beta,\mu)$
		$\displaystyle+\frac{1}{4}\sum_{pqrs}t_{pqrs}a_{p}^{\dagger}(\beta,\mu)a_{q}^{\dagger}(\beta,\mu)\tilde{a}_{s}^{\dagger}(\beta,\mu)\tilde{a}_{r}^{\dagger}(\beta,\mu).$		(14)

While $t_{0}$ , $t_{pq}$ , and $t_{pqrs}$ also depend on $\beta$ and $\mu$ , we have suppressed the explicit dependence for brevity.

IV.2.2 Coupled cluster theory

In coupled cluster theory, the wave operator is expressed as an exponential of the quasiparticle excitations. At the level of singles and doubles approximation (CCSD), the thermal state is defined as

\ket{\Psi_{CC}(\beta,\mu)}=e^{t_{0}+T}\ket{\Phi(\beta,\mu)},

(15)

where $t_{0}$ is a scalar and $T$ has the same form as in Eq. 14. Technical details on the derivation and exact form of evolution equations for the amplitudes $t_{pq}$ and $t_{pqrs}$ in CISD and CCSD can be found in Refs. 28, 29.

IV.2.3 Computational scaling and accuracy

The asymptotic computational scaling of thermal wave function methods is the same as their ground-state counterparts, albeit with a pre-factor proportional to the grids used in integration along $\beta$ - and $\mu$ -directions. Therefore, while finite-temperature mean-field scales as $\mathcal{O}(n^{3})$ , correlated methods such as CISD and CCSD scale as $\mathcal{O}(n^{6})$ , where $n$ is the number of spin-orbitals used to describe the system. Additionally, the scaling of ground-state methods often benefit from a clear distinction between the occupied and unoccupied (or virtual) orbitals, e.g., ground-state CCSD scales as $\mathcal{O}(n_{o}^{2}n_{v}^{4})$ , where $n_{o}$ and $n_{v}$ are the number of occupied and virtual orbitals. At finite temperatures, all orbitals are treated equally, and a distinction based on orbital occupation does not exist. In fact, considering that finite-temperature mean-field is a BCS-style state, thermal CI and CC are equivalent in formulation and scaling to ground-state quasiparticle CI and CC. [47]

While higher order approximations have not been implemented so far, we expect that the thermal CI and CC will become more accurate as quasiparticle excitations at the level of triples, quadruples, etc. are incorporated. However, it is interesting to note that in the grand-canonical ensemble, finite-temperature CI and CC will become exact only if quasiparticle excitations up to the order equal to the number of spin-orbitals are used. For instance, in Ref. 29, we showed that while ground-state CCSD and CISD are exact for two-electron systems, finite-temperature CISD and CISD are not. This again is a property inherited from quasiparticle CI and CC. [47]

IV.3 Physical properties

With the correlated thermal states, the finite-temperature expectation value of a physical observable, represented by an operator $A$ , can be computed using Eq. 3. For instance, the CI approximation to $\braket{A}(\beta,\mu)$ is defined as

\braket{A}_{CI}(\beta,\mu)=\frac{\braket{\Psi_{CI}(\beta,\mu)}{A}{\Psi_{CI}(\beta,\mu)}}{\braket{\Psi_{CI}(\beta,\mu)}{\Psi_{CI}(\beta,\mu)}}.

(16)

On the other hand, for the coupled cluster method, a symmetric expectation value is generally not feasible. Taking inspiration from ground-state CC, we use a linear, or CI, approximation for the bra thermal state, i.e.,


$\displaystyle\bra{\Psi_{CI}(\beta,\mu)}$	$\displaystyle=\bra{\Phi(\beta,\mu)}\left(1+W(\beta,\mu)\right)e^{w_{0}},$	(17a)
	$\displaystyle=\bra{\Phi(\beta,\mu)}\left(1+Z(\beta,\mu)\right)e^{z_{0}-T(\beta,\mu)},$	(17b)

where the operators $W(\beta,\mu)$ and $Z(\beta,\mu)$ are quasiparticle de-excitation operators (excitations on the bra), and are defined as


$\displaystyle W(\beta,\mu)$	$\displaystyle=\sum_{pq}w_{pq}\tilde{a}_{q}(\beta,\mu)a_{p}(\beta,\mu)$
$\displaystyle+\frac{1}{4}\sum_{pqrs}$	$\displaystyle w_{pqrs}\tilde{a}_{r}(\beta,\mu)\tilde{a}_{s}(\beta,\mu)a_{q}(\beta,\mu)a_{p}(\beta,\mu),$	(18a)
$\displaystyle Z(\beta,\mu)$	$\displaystyle=\sum_{pq}z_{pq}\tilde{a}_{q}(\beta,\mu)a_{p}(\beta,\mu)$
$\displaystyle+\frac{1}{4}\sum_{pqrs}$	$\displaystyle z_{pqrs}\tilde{a}_{r}(\beta,\mu)\tilde{a}_{s}(\beta,\mu)a_{q}(\beta,\mu)a_{p}(\beta,\mu).$	(18b)

Once again, the explicit $\beta$ - and $\mu$ -dependence in the $w$ - and $z$ -amplitudes has been suppressed for brevity. Using the fact that $\bra{\Phi(\beta,\mu)}e^{T(\beta,\mu)}=\bra{\Phi(\beta,\mu)}$ , a transformation between $Z$ and $W$ can be derived:


$\displaystyle\bra{\Phi(\beta,\mu)}(1+W(\beta,\mu))e^{w_{0}}$	$\displaystyle=\bra{\Phi(\beta,\mu)}e^{-T(\beta,\mu)}(1+W(\beta,\mu))e^{w_{0}}e^{T(\beta,\mu)}e^{-T(\beta,\mu)},$	(19a)
$\displaystyle\Rightarrow\bra{\Phi(\beta,\mu)}(1+Z(\beta,\mu))e^{z_{0}}$	$\displaystyle=\bra{\Phi(\beta,\mu)}e^{-T(\beta,\mu)}(1+W(\beta,\mu))e^{w_{0}}e^{T(\beta,\mu)}.$	(19b)

Finally, with a linearized approximation to the bra thermal state, we define the finite-temperature CC expectation values as


$\displaystyle\braket{A}_{CC}(\beta,\mu)$	$\displaystyle=\frac{\braket{\Psi_{CI}(\beta,\mu)}{A}{\Psi_{CC}(\beta,\mu)}}{\braket{\Psi_{CI}(\beta,\mu)}{\Psi_{CC}(\beta,\mu)}},$	(20a)
	$\displaystyle=\frac{\braket{\Phi(\beta,\mu)}{\left(1+Z(\beta,\mu)\right)e^{-T(\beta,\mu)}Ae^{T(\beta,\mu)}}{\Phi(\beta,\mu)}}{\braket{\Phi(\beta,\mu)}{\Phi(\beta,\mu)}}.$	(20b)

The grand potential cannot be computed as an operator expectation value. However, it is directly related to the partition function which, in thermofield dynamics, is given by the normalization of the thermal wave function. The finite-temperature coupled cluster grand potential can thus be computed as


$\displaystyle\Omega_{CC}$	$\displaystyle=\Omega_{0}+\Omega_{c},$	(21a)
$\displaystyle\Omega_{0}$	$\displaystyle=-\frac{1}{\beta}\log\braket{\mathbb{I}}{e^{-\beta(H_{0}-\mu N)}}{\mathbb{I}}=-\frac{1}{\beta}\log\left[\prod_{p}\left(1+e^{-\beta(\epsilon_{p}-\mu)}\right)\right],$	(21b)
$\displaystyle\Omega_{c}$	$\displaystyle=-\frac{1}{\beta}\log\left[\frac{\braket{\Phi(\beta,\mu)}{(1+Z(\beta,\mu))e^{z_{0}}e^{-T(\beta,\mu)}e^{t_{0}+T(\beta,\mu)}}{\Phi(\beta,\mu)}}{\braket{\Phi(\beta,\mu)}{\Phi(\beta,\mu)}}\right]=-\frac{1}{\beta}\left(t_{0}+z_{0}\right).$	(21c)

Here, $\Omega_{0}$ and $\Omega_{c}$ denote the the mean-field and correlation contributions to the grand potential. In simplifying Eq. 21b, we recall our definition for the mean-field Hamiltonian, $H_{0}=\sum_{p}\epsilon_{p}c_{p}^{\dagger}c_{p}$ . Other thermodynamic properties can be derived from the grand potential. For example, given the coupled cluster internal energy $E=\braket{H}_{CC}$ and particle number $N=\braket{N}_{CC}$ , the entropy is defined as:

S_{CC}=\frac{1}{T}\left(\braket{H}_{CC}-\mu\braket{N}_{CC}-\Omega_{CC}\right),

(22)

where the internal energy, $E=\braket{H}_{CC}$ , and the number of particles, $N=\braket{N}_{CC}$ , are simply the thermal expectation values of the Hamiltonian and the number operator, We can also define the particle number, entropy, and the internal energy as derivatives of the grand potential:


$\displaystyle N$	$\displaystyle=-\frac{\partial\Omega}{\partial\mu},$	(23a)
$\displaystyle S$	$\displaystyle=\beta^{2}\frac{\partial\Omega}{\partial\beta},$	(23b)
$\displaystyle E$	$\displaystyle=\Omega+\mu N+\frac{1}{\beta}S.$	(23c)

In principle, the derivatives of the grand potential (Eq. 23) and the expectation values (e.g., Eq. 22) will result in properties that are numerically identical only when the thermal bra and ket make the grand potential stationary. [32] For our finite-temperature coupled cluster theory, where we time-evolve the CI-like bra and the CC-like ket separately, the stationarity of the grand potential is not guaranteed. However, in practice, the derivative and expectation value formalisms lead to very similar results (see, for example, the comparison of magnetization properties in the one-dimensional transverse field Ising model in Ref. 32)

V Results

To study the accuracy of various approximations to the thermal state, we use the one-dimensional Hubbard model [48] as a benchmark system, the Hamiltonian for which is defined as

H=-t\sum_{\langle p,q\rangle,\sigma}\left(c_{p,\sigma}^{\dagger}\,c_{q,\sigma}+\textrm{h.c.}\right)+U\,\sum_{p}n_{p,\uparrow}\,n_{p,\downarrow},

(24)

where $\langle,\rangle$ denotes that the sum is carried over connected lattice sites, $t$ denotes the hopping strength, $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 physics of the Hubbard model depends only on the ratio $U/t$ , which characterizes the correlation strength.

The left panel in Fig. 2 shows the total internal energy for the mean-field theory, CISD, CCSD, and exact thermal wave functions for the 6-site Hubbard model with $U/t=2$ at half-filling on average. The right panel shows the errors in internal energy for the approximate methods against the exact result. The error plot clearly show that thermal CISD and CCSD significantly improve the results as compared to the mean-field theory, with CCSD outperforming the other two. The markers in the y-axis of the error plot indicate the errors in the corresponding ground-state theories, which establishes that the thermal wave function methods approach the ground-state methods in the limit $\beta\rightarrow\infty$ (or $T\rightarrow 0$ ).

A common feature in ground-state methods is that the results in the correlated methods depend on the choice of mean-field reference state, e.g., a symmetry adapted restricted Hartree-Fock vs. a broken-symmetry unrestricted Hartree-Fock state. Figure 3 shows a similar behavior for thermal wave function methods. Here, we have used the 6-site Hubbard model with $U/t=5$ at half-filling on average. The symmetry adapted thermal mean-field was constructed using zero-temperature restricted Hartree-Fock eigenvalues and orbitals as our $H_{0}$ , while the broken symmetry mean-field was constructed from zero-temperature unrestricted Hartree-Fock. The coupled cluster results were built upon the respective reference states. The results once again show that thermal CCSD improves significantly over the respective mean-field reference, although the results depend on the choice of mean-field reference. Both the symmetry adapted and broken-symmetry thermal mean-field and CCSD approach their ground-state counterparts as we approach zero temperature.

In Fig. 4, we also compare mean-field, CISD, and CCSD approximations to the grand potential for 6-site Hubbard models at $U/t=2$ and at $U/t=5$ . For the latter, we use both the symmetry adapted and broken symmetry thermal mean-field reference, similar to Fig. 3. For the comparison of grand potentials, we fix the chemical potential ( $\mu=0$ ) rather than fixing the total number of electrons. As expected, we find that both CISD and CCSD improve significantly over the mean-field grand potential. While other mean-field observables such as internal energy, entropy, etc. are, by definition, exact in the limit $\beta\rightarrow 0$ (i.e., infinite temperature), the grand potential $\Omega_{0}$ defined in Eq. 21b is not exact in the high-temperature limit. This is because $\Omega_{0}$ includes contribution only from the mean-field Hamiltonian $H_{0}$ . For $U/t=5$ , it is well known that ground-state CCSD, whether symmetry adapted or symmetry broken, struggles to recover the correlation energy. This ineffectiveness of single-reference methods like CI and CC is also visible at finite temperatures.

At finite temperatures, we are more interested in properties other than internal energy. For example, in the one-dimensional Hubbard model, the spin-spin correlation function, defined as

\chi(i,j)=\langle S^{z}(i)\>S^{z}(j)\rangle

(25)

is an observable quantification of the correlation in the system. For a weakly interacting Hubbard model, the strength of this correlation function decays very rapidly with the distance between two lattice sites. On the other hand, as the interaction strength is increased, the length-scales over which $\chi(i,j)$ decays increases. We use the biorthogonal framework outlined in Eq. 20 to evaluate these correlation functions. Figure 5 shows the trends in $\chi(i,j)$ at various temperatures for the 10-site Hubbard model with $U/t=2$ . A fixed chemical potential was used for this calculation which ensures near-half-filling for all temperature points used. The correlation function for exact and coupled cluster calculations for the ground-state (i.e., at $T=0$ ) is also shown for reference. We can observe that as the temperature is increased, the correlation strength as well as the correlation length decreases, which establishes that thermal CCSD is capable of accurately describing trends in the order parameters as a function of temperature.

V.1 Zero-temperature limit of thermal methods

Figures 2 and 3 show that all the approximate finite-temperature methods converge to their ground-state counterparts in the zero-temperature limit. However, the well behaved zero-temperature convergence of these methods arises from the fact that the thermal mean-field converges to the appropriate reference determinant, i.e.

\lim_{\beta\rightarrow\infty}\ket{\Phi(\beta)}=\ket{\Phi_{0}},

(26)

where $\ket{\Phi_{0}}$ is the ground-state reference Slater determinant. This limit is well behaved only if the ground-state of the mean-field driver Hamiltonian $H_{0}=\sum_{p}\epsilon_{p}c^{\dagger}_{p}c_{p}$ is non-degenerate. Otherwise, the zero-temperature limit of $\ket{\Phi(\beta)}$ is not guaranteed to be a single Slater determinant, and consequently, the correlated methods would not necessarily converge to single-reference ground-state counterparts. For a repulsive Hamiltonian, even if the exact ground-state is degenerate, it is known that the general Hartree-Fock eigenvalues will not be degenerate. [49] We would also like to note that degeneracies caused by spin are generally not a problem, because $H_{0}$ is defined by a zero-temperature Fock operator which tends to lift these degeneracies even for unrestricted or restricted open-shell Hartree-Fock. Consider, for example, the hydrogen atom where we have placed the electron in the $1s\uparrow$ spin state, whose orbital energy is given by

\epsilon_{1s\uparrow}=\braket{1s\uparrow}{h}{1s\uparrow}

(27)

since the self-Coulomb is cancelled by the self-exchange in Hartree-Fock theory. On the other hand, the orbital energy of the $1s\downarrow$ state is raised by the Coulomb interaction with the $1s\uparrow$ orbital:

\epsilon_{1s\downarrow}=\braket{1s\downarrow}{h}{1s\downarrow}+\braket{1s\downarrow,1s\uparrow}{v}{1s\downarrow,1s\uparrow}.

(28)

Here, $h$ and $v$ denote the one- and two-electron interaction terms in the Hamiltonian. Because the two orbitals are not degenerate, the zeroth-order Hamiltonian has a non-degenerate ground state to which, as $\beta\rightarrow\infty$ , the thermal mean-field state converges. In this case, our thermal mean-field will converge to the $\uparrow$ -spin ground state.

Furthermore, in most scenarios with near-degeneracies in the mean-field Hamiltonian, ground-state single-reference coupled cluster breaks down, and the question of a ground-state limit for the finite-temperature method does not arise. In such cases, we observe a divergence in the thermal CC amplitudes as we evolve towards zero temperature.

VI Connections with other finite-temperature CC theories

In Section II, we introduced thermal CC theories by Sanyal et al., Hummel, and White et al., all of which are inspired by thermal Wick’s theorem and imaginary-time-dependent perturbation theory. While the thermofield version is completely different in its formulation, the resulting equations for the finite-temperature CC are identical to the other three methods, provided that other implementation details, such as the mean-field driver $H_{0}$ and the equation for the CC-bra, are consistent. The similarity is not surprising considering the fact that thermofield dynamics offers a unique way of computing thermal traces, much like the thermal Wick’s theorem.

However, an explicit wave function representation of the thermal density matrix renders the thermofield approach far more versatile. For instance, Ancilla DMRG, which is also based on the purification of an ensemble density matrix, is a thermal generalization of the matrix product states. Similarly, wave functions such as multi-reference CI and CC and active space CI, for which diagrammatic equations and Wick’s theorem are not straightforward, can also be extended to finite temperatures using our formalism. This is currently a work in progress. Moreover, the different thermal CC methods discussed above work in the grand canonical ensemble. While the finite-temperature Wick’s theorem and diagrammatic perturbation theory do not have an analog in the canonical ensemble, the wave function in thermofield dynamics is amenable to ground-state number projection techniques, such that a canonical ensemble formulation is straightforward and has been explored by the authors in Ref. 31.

VII Conclusions

In this short review, we have presented recent developments in finite-temperature wave function theories, focussing particularly on the formalism inspired by thermofield dynamics. Within this approach, finite-temperature generalizations of the whole hierarchy of ground-state methods can be constructed. We have presented benchmark calculations on a simple, yet non-trivial, one-dimensional Hubbard model. By comparing errors in internal energy, grand-potential and correlation functions, we confirm that finite-temperature wave functions behave similarly to their ground-state counterparts. For moderately correlated electronic systems, both thermal CI and CC truncated at single and double excitations provide a significant improvement over mean-field. While our implementation is capable of studying larger and more realistic electronic systems, the absence of exact benchmark data makes it difficult to assess the performance of the thermal wave function theories. A comparative study involving different finite-temperature methods would be a natural step forward.

The explicit wave function representation of the finite-temperature Gibbs-Boltzmann density matrix allows the potential of extending this formalism not only to multi-reference methods, but also to emerging algorithms for quantum computers. Application of, say, thermal CC to periodic systems is another venue that deserves more attention. It would largely facilitate the study of temperature-correlation phase diagrams in realistic materials. Due to its size-extensive property, we expect thermal CCSD to perform well, particularly in comparison to mean-field and CISD, as we go to larger systems. The finite-temperature methods inherit most of the features from their ground-state counterparts such as computational scaling, dependence of CC and CI on the reference state, etc. Accordingly, we expect that thermal CC and CI theories would be most effective for weakly correlated systems. The emerging field of finite-temperature wave function theories shows a lot of promise to compete with, as well as to supplement, presently standard methods such as quantum Monte Carlo and Matsubara perturbation theory.

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

References

Lee et al. [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a Mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78, 17 (2006).
Raveau [2011] B. Raveau, Strongly correlated electron systems: From chemistry to physics, Comptes Rendus Chimie Future of sciences, sciences for the future: Chemistry and its interfaces with biology and physics, 14, 856 (2011).
Mukherjee et al. [2013] S. Mukherjee, F. Libisch, N. Large, O. Neumann, L. V. Brown, J. Cheng, J. B. Lassiter, E. A. Carter, P. Nordlander, and N. J. Halas, Hot Electrons Do the Impossible: Plasmon-Induced Dissociation of H2 on Au, Nano Lett. 13, 240 (2013), publisher: American Chemical Society.
Guillot [1999] T. Guillot, Interiors of giant planets inside and outside the solar system, Science 286, 72 (1999).
Wann et al. [2017] E. T. H. Wann, L. Vočadlo, and I. G. Wood, High-temperature ab initio calculations on FeSi and NiSi at conditions relevant to small planetary cores, Phys Chem Minerals 44, 477 (2017).
Mermin [1963] N. D. Mermin, Stability of the thermal Hartree-Fock approximation, Ann. Phys. 21, 99 (1963).
Sokoloff [1967] J. Sokoloff, Some consequences of the thermal Hartree-Fock approximation at zero temperature, Ann. Phys. 45, 186 (1967).
Matsubara [1955] T. Matsubara, A New Approach to Quantum-Statistical Mechanics, Prog. Theor. Phys. 14, 351 (1955).
Crawford and Schaefer [2000] T. D. Crawford and H. F. Schaefer, An Introduction to Coupled Cluster Theory for Computational Chemists, in Reviews in Computational Chemistry, edited by K. B. Lipkowitz and D. B. Boyd (John Wiley & Sons, Inc., 2000) pp. 33–136.
Bartlett and Musiał [2007] R. J. Bartlett and M. Musiał, Coupled-cluster theory in quantum chemistry, Rev. Mod. Phys. 79, 291 (2007).
Sun et al. [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, et al., Pyscf: the python-based simulations of chemistry framework, WIREs Computational Molecular Science 8, e1340 (2018).
Giannozzi et al. [2020] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, et al., Quantum espresso toward the exascale, J. Chem. Phys. 152, 154105 (2020).
Lichtenstein et al. [2001] A. I. Lichtenstein, M. I. Katsnelson, and G. Kotliar, Finite-Temperature Magnetism of Transition Metals: An ab initio Dynamical Mean-Field Theory, Phys. Rev. Lett. 87, 067205 (2001).
Verstraete et al. [2004] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Matrix product density operators: Simulation of finite-temperature and dissipative systems, Phys. Rev. Lett. 93, 207204 (2004).
Feiguin and White [2005] A. E. Feiguin and S. R. White, Finite-temperature density matrix renormalization using an enlarged Hilbert space, Phys. Rev. B 72, 220401 (2005).
Czarnik et al. [2016] P. Czarnik, M. M. Rams, and J. Dziarmaga, Variational tensor network renormalization in imaginary time: Benchmark results in the Hubbard model at finite temperature, Phys. Rev. B 94, 235142 (2016).
White [2009] S. R. White, Minimally Entangled Typical Quantum States at Finite Temperature, Phys. Rev. Lett. 102, 190601 (2009).
Stoudenmire and White [2010] E. M. Stoudenmire and S. R. White, Minimally entangled typical thermal state algorithms, New J. Phys. 12, 055026 (2010).
Pittalis et al. [2011] S. Pittalis, C. R. Proetto, A. Floris, A. Sanna, C. Bersier, K. Burke, and E. K. U. Gross, Exact Conditions in Finite-Temperature Density-Functional Theory, Phys. Rev. Lett. 107, 163001 (2011).
Santra and Schirmer [2017] R. Santra and J. Schirmer, Finite-temperature second-order many-body perturbation theory revisited, Chem. Phys. Electrons and nuclei in motion - correlation and dynamics in molecules (on the occasion of the 70th birthday of Lorenz S. Cederbaum), 482, 355 (2017).
Hirata and Jha [2019] S. Hirata and P. K. Jha, Chapter Two - Converging finite-temperature many-body perturbation theory in the grand canonical ensemble that conserves the average number of electrons, in Annual Reports in Computational Chemistry, Vol. 15, edited by D. A. Dixon (Elsevier, 2019) pp. 17–37.
Sanyal et al. [1992] G. Sanyal, S. H. Mandal, and D. Mukherjee, Thermal averaging in quantum many-body systems: a non-perturbative thermal cluster cumulant approach, Chem. Phys. Lett. 192, 55 (1992).
Sanyal et al. [1993] G. Sanyal, S. H. Mandal, S. Guha, and D. Mukherjee, Systematic nonperturbative approach for thermal averages in quantum many-body systems: The thermal-cluster-cumulant method, Phys. Rev. E 48, 3373 (1993).
Mandal et al. [2003] S. H. Mandal, R. Ghosh, G. Sanyal, and D. Mukherjee, A finite-temperature generalisation of the coupled cluster method: a non-perturbative access to grand partition functions, Int. J. Mod. Phys. B 17, 5367 (2003).
Hermes and Hirata [2015] M. R. Hermes and S. Hirata, Finite-temperature coupled-cluster, many-body perturbation, and restricted and unrestricted Hartree–Fock study on one-dimensional solids: Luttinger liquids, Peierls transitions, and spin- and charge-density waves, J. Chem. Phys. 143, 102818 (2015).
White and Chan [2018] A. F. White and G. K.-L. Chan, A Time-Dependent Formulation of Coupled-Cluster Theory for Many-Fermion Systems at Finite Temperature, J. Chem. Theory Comput. 14, 5690 (2018).
Hummel [2018] F. Hummel, Finite Temperature Coupled Cluster Theories for Extended Systems, J. Chem. Theory Comput. 14, 6505 (2018), pMID: 30354112, https://doi.org/10.1021/acs.jctc.8b00793 .
Harsha et al. [2019a] G. Harsha, T. M. Henderson, and G. E. Scuseria, Thermofield theory for finite-temperature quantum chemistry, J. Chem. Phys. 150, 154109 (2019a).
Harsha et al. [2019b] G. Harsha, T. M. Henderson, and G. E. Scuseria, Thermofield theory for finite-temperature coupled cluster, J. Chem. Theory Comput. 15, 6127 (2019b).
Shushkov and Miller [2019] P. Shushkov and T. F. Miller, Real-time density-matrix coupled-cluster approach for closed and open systems at finite temperature, J. Chem. Phys. 151, 134107 (2019).
Harsha et al. [2020] G. Harsha, T. M. Henderson, and G. E. Scuseria, Wave function methods for canonical ensemble thermal averages in correlated many-fermion systems, J. Chem. Phys. 153, 124115 (2020).
Harsha et al. [2022] G. Harsha, Y. Xu, T. M. Henderson, and G. E. Scuseria, Thermal coupled cluster theory for SU(2) systems, Phys. Rev. B 105, 045125 (2022), publisher: American Physical Society.
Matsumoto et al. [1983] H. Matsumoto, Y. Nakano, H. Umezawa, F. Mancini, and M. Marinaro, Thermo Field Dynamics in Interaction Representation, Prog. Theor. Phys. 70, 599 (1983).
Semenoff and Umezawa [1983] G. W. Semenoff and H. Umezawa, Functional methods in thermofield dynamics: A real-time perturbation theory for quantum statistical mechanics, Nucl. Phys. B 220, 196 (1983).
Umezawa [1984] H. Umezawa, Methods of Quantum Field Theory in Condensed Matter Physics— New Perspectives, Extensions and Applications —, Prog. Theor. Phys. 80, 26 (1984).
Keldysh [1965] L. P. Keldysh, Diagram technique for nonequilibrium processes, Sov. Phys. JETP 20, 1018 (1965).
Schmidt [1907] E. Schmidt, Zur Theorie der linearen und nichtlinearen Integralgleichungen, Math. Ann. 63, 433 (1907).
Suzuki [1985] M. Suzuki, Thermo Field Dynamics in Equilibrium and Non-Equilibrium Interacting Quantum Systems, J. Phys. Soc. Jpn. 54, 4483 (1985).
Hatsuda [1989] T. Hatsuda, Mean field theory and boson expansion at finite temperature on the basis of the thermo field dynamics, Nucl. Phys. A 492, 187 (1989).
Walet and Klein [1990] N. R. Walet and A. Klein, Thermal boson expansions and dynamical symmetry, Nucl. Phys. A 510, 261 (1990).
de Vega and Bañuls [2015] I. de Vega and M.-C. Bañuls, Thermofield-based chain-mapping approach for open quantum systems, Phys. Rev. A 92, 052116 (2015).
Borrelli and Gelin [2016] R. Borrelli and M. F. Gelin, Quantum electron-vibrational dynamics at finite temperature: Thermo field dynamics approach, J. Chem. Phys. 145, 224101 (2016).
Nocera and Alvarez [2016] A. Nocera and G. Alvarez, Symmetry-conserving purification of quantum states within the density matrix renormalization group, Phys. Rev. B 93, 045137 (2016).
Chen and Zhao [2017] L. Chen and Y. Zhao, Finite temperature dynamics of a Holstein polaron: The thermo-field dynamics approach, J. Chem. Phys. 147, 214102 (2017).
Borrelli and Gelin [2017] R. Borrelli and M. F. Gelin, Simulation of Quantum Dynamics of Excitonic Systems at Finite Temperature: an efficient method based on Thermo Field Dynamics, Sci. Rep. 7, 9127 (2017).
Wu and Hsieh [2019] J. Wu and T. H. Hsieh, Variational Thermal Quantum Simulation via Thermofield Double States, Phys. Rev. Lett. 123, 220502 (2019), publisher: American Physical Society.
Henderson et al. [2014] T. M. Henderson, G. E. Scuseria, J. Dukelsky, A. Signoracci, and T. Duguet, Quasiparticle coupled cluster theory for pairing interactions, Phys. Rev. C 89, 054305 (2014).
Hubbard J. and Flowers Brian Hilton [1963] Hubbard J. and Flowers Brian Hilton, Electron correlations in narrow energy bands, P. R. Soc. Lond. A Mat. 276, 238 (1963).
Bach et al. [1994] V. Bach, E. H. Lieb, M. Loss, and J. P. Solovej, There are no unfilled shells in unrestricted Hartree-Fock theory, Phys. Rev. Lett. 72, 2981 (1994), publisher: American Physical Society.