Thermofield Theory for Finite-Temperature Quantum Chemistry

TL;DR

This paper introduces a thermofield dynamics framework for finite-temperature quantum chemistry, enabling the development of thermal wave function methods and applying them to the Hubbard model with promising results.

Contribution

It extends thermofield dynamics to quantum chemistry, proposing new finite-temperature correlated wave function methods and two approaches to imaginary time Schrödinger equations.

Findings

Derived thermal configuration interaction theories

Applied methods to the Hubbard model

Compared results with exact benchmarks

Abstract

Thermofield dynamics has proven to be a very useful theory in high-energy physics, particularly since it permits the treatment of both time- and temperature-dependence on an equal footing. We here show that it also has an excellent potential for studying thermal properties of electronic systems in physics and chemistry. We describe a general framework for constructing finite temperature correlated wave function methods typical of ground state methods. We then introduce two distinct approaches to the resulting imaginary time Schrodinger equation, which we refer to as fixed-reference and covariant methods. As an example, we derive the two corresponding versions of thermal configuration interaction theory, and apply them to the Hubbard model, while comparing with exact benchmark results.

Figures

Equations

$$⟨A{⟩}_{\beta }$$

$$=\frac{1}{\mathcal{Z}}\sum _m⟨m\left|A{\mathrm{e}}^{-\beta H}\right|m⟩,$$

$$\hat{\rho }={\mathrm{e}}^{-\beta H},\mathcal{Z}=\sum _m⟨m\left|\hat{\rho }\right|m⟩.$$

$$⟨\hat{A}{⟩}_{\beta }=\frac{⟨\Psi (\beta )|\hat{A}|\Psi (\beta )⟩}{⟨\Psi (\beta )|\Psi (\beta )⟩}.$$

$$c_k,c_k^{†};{\tilde{c}}_k,{\tilde{c}}_k^{†}$$

$$[c_k,c_k^{†}{]}_{\mp }$$

$$[c_k,{\tilde{c}}_k{]}_{\mp }$$

$$[A,B{]}_{-\eta }=AB-\eta BA,$$

$$\widetilde{\left({\tilde{c}}_k\right)}$$

$$\widetilde{\left(\alpha c_k+\delta c_q^{†}\right)}$$

$$\widetilde{\left(c_k{c}_q^{†}\right)}$$

$$-\mathrm{i}\frac{\partial }{\partial t}|\psi ⟩=\tilde{H}|\psi ⟩.$$

$$|\Psi (\beta )⟩={\mathrm{e}}^{-\beta H/2}\sum _m|m⟩\otimes |\tilde{m}⟩,$$

$$\mathcal{Z}=⟨\Psi (\beta )|\Psi (\beta )⟩.$$

$$|\mathbb{I}⟩=|\Psi (\beta =0)⟩=\sum _m|m⟩\otimes |\tilde{m}⟩.$$

$$⟨\hat{A}{⟩}_{\beta }$$

$$=\frac{1}{\mathcal{Z}}\sum _{m,n}⟨m,\tilde{m}|{\mathrm{e}}^{-\beta H/2}\hat{A}{\mathrm{e}}^{-\beta H/2}|n,\tilde{n}⟩$$

$$=\frac{1}{\mathcal{Z}}\sum _m⟨m|{\mathrm{e}}^{-\beta H/2}\hat{A}{\mathrm{e}}^{-\beta H/2}|m⟩$$

$$=\frac{1}{\mathcal{Z}}\mathrm{Tr}\left({\mathrm{e}}^{-\beta H/2}\hat{A}{\mathrm{e}}^{-\beta H/2}\right)$$

$$=\frac{1}{\mathcal{Z}}\mathrm{Tr}\left({\mathrm{e}}^{-\beta H}\hat{A}\right)$$

$$H_{\text{th}}=H-\tilde{H},$$

$$\mathrm{i}\frac{\partial }{\partial t}|\Psi (\beta )⟩=H_{\text{th}}|\Psi (\beta )⟩.$$

$$\frac{\partial }{\partial \beta }|\Psi (\beta )⟩=-\frac{1}{2}H|\Psi (\beta )⟩,$$

$$|\Psi (\alpha ,\beta )⟩$$

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

$${|\Psi (\beta )⟩}_{\text{mf}}=|0(\beta )⟩={\mathrm{e}}^{-\beta H_0/2}\sum _m|m⟩\otimes |\tilde{m}⟩,$$

$$\frac{\partial }{\partial \beta }|0(\beta )⟩=-\frac{1}{2}{H}_0|0(\beta )⟩.$$

$$|0(\beta )⟩=\sum _m{\mathrm{e}}^{-\beta E_m/2}|m⟩\otimes |\tilde{m}⟩.$$

$${\mathcal{Z}}_0=⟨0(\beta )|0(\beta )⟩.$$

$$\left\lvert 0(\beta)\right\rangle=\mathrm{e}^{-\beta H_{0}/2}\prod_{p\in\textrm{levels}}\Big{(}|0\rangle_{p}\otimes|\tilde{0}\rangle_{p}+|1\rangle_{p}\otimes|\tilde{1}\rangle_{p}\Big{)}$$

$$\left[\begin{array}{c}{a}_p(\beta )\\ {\tilde{a}}_p^{†}(\beta )\end{array}\right]=\left[\begin{array}{cc}{w}_p& -z_p\\ z_p& w_p\end{array}\right]\left[\begin{array}{c}{c}_p\\ {\tilde{c}}_p^{†}\end{array}\right],$$

$$a_p|0(\beta )⟩$$

$$⟨0(\beta )|a_p^{†}$$

$$H=\sum _{p,q}{h}_{pq}{c}_p^{†}{c}_q+\frac{1}{4}\sum _{p,q}{v}_{pqrs}{c}_p^{†}{c}_q^{†}{c}_s{c}_r,$$

$$H$$

Full text

I Introduction

The Schrödinger equation for an interacting many-electron system can almost never be solved exactly due to the exponential scaling of the Hilbert space.Dirac (1929) Over the years, several methods have been proposed to compute approximate solutions and properties for such systems. These methods (perhaps with the exception of density functional theory) fall into two broad categories: the path integral approach and the wave function approach. The path integral formulation, generally used in condensed matter physics, is particularly useful when thermodynamically large systems can be approximated using an effective mean-field theory in the presence of perturbing interactions. On the other hand, wave function methods (such as configuration interaction (CI), coupled clusterCrawford and Schaefer (2000); Bartlett and Musiał (2007) (CC), or matrix product statesÖstlund and Rommer (1995); Dukelsky et al. (1998) (MPS)) are often preferred for the study of finite many-body systems such as atoms and molecules – systems generally of interest in chemistry.

While the path integral formulation can be naturally extended to study the thermal behaviour of a system via Matsubara’s imaginary time formalism,Matsubara (1955) the usage of typical wave function methods breaks down at non-zero temperatures. This is because at non-zero temperatures, a quantum system is described by an ensemble density matrix, and one would need to solve for the entire spectrum of the Hamiltonian, as opposed to finding just the ground state energy and the corresponding wave function at zero temperature. There are a number of reasons for developing finite temperature quantum chemistry, a subject that has recently attracted considerable interest.Verstraete, García-Ripoll, and Cirac (2004); Feiguin and White (2005); White (2009); Stoudenmire and White (2010); Pižorn et al. (2014); Hermes and Hirata (2015); Czarnik, Rams, and Dziarmaga (2016); Santra and Schirmer (2017); Zgid and Gull (2017); Claes and Clark (2017); White and Chan (2018); Hummel (2018); Hirata and Jha (2018) The finite temperature variants of the coupled cluster methodHummel (2018); White and Chan (2018) are of particular relevance to our present work. These techniques are based on an extension of the time-dependent CC method to imaginary time, and are framed loosely along similar lines as the Thermal Cluster Cumulant theory proposed by Mukherjee et.al.,Sanyal, Mandal, and Mukherjee (1992); Sanyal et al. (1993); Mandal, Sanyal, and Mukherjee (1998); Mandal et al. (2003) which uses a finite-temperature generalization of Wick’s theorem to introduce a thermal normal ordering. The interaction-picture Boltzmann operator is then compactly expressed as a thermal normal ordered exponential of some cluster operators and a number. Such an exponential ansatz is naturally size-extensive and forms a thermal analogue of ground-state CC.

Umezawa et. al. proposed an alternative approach, known as thermofield dynamics,Matsumoto et al. (1983); Semenoff and Umezawa (1983); Umezawa (1984); Evans et al. (1992) which provides a real-time approach as opposed to Matsubara’s imaginary time formalism. Thermofield dynamics (TFD) can be conveniently described in an operator or wave function formulation, and was originally proposed as a method to study time-dependent and non-equilibrium phenomena in many-body quantum systems, something which cannot be accomplished easily within the original Matsubara formalism, but is addressed within the Keldysh formalism.Keldysh (1965) TFD has been used widely in high energy physicsMaldacena (2003); Shenker and Stanford (2014); Israel (1976); Yang (2018) as well as for the study of time-dependent and open quantum systems.Suzuki (1985); de Vega and Bañuls (2015); Borrelli and Gelin (2016); Chen and Zhao (2017); Borrelli and Gelin (2017)

The wavefunction prescription of TFD has tremendous potential for studying thermal properties of quantum many body systems in physics and chemistry. However, its application to compute equilibrium thermal properties has been relatively unexplored and generally limited to the mean-field level.Hatsuda (1989); Walet and Klein (1990) Nonetheless, there are some methods such as Ancilla Density Matrix Renormalization GroupFeiguin and White (2005); Nocera and Alvarez (2016) which try to harness the novel features of TFD. In this paper, we present a general framework based on TFD to construct wave function methods, extending traditional ground-state methods, to compute thermal averages of physical quantities. Our TFD-based approach provides an alternative to existing techniques.Verstraete, García-Ripoll, and Cirac (2004); Feiguin and White (2005); White (2009); Stoudenmire and White (2010); Pižorn et al. (2014); Hermes and Hirata (2015); Czarnik, Rams, and Dziarmaga (2016); Santra and Schirmer (2017); Zgid and Gull (2017); Claes and Clark (2017); White and Chan (2018); Hummel (2018); Hirata and Jha (2018); Sanyal, Mandal, and Mukherjee (1992); Sanyal et al. (1993); Mandal, Sanyal, and Mukherjee (1998); Mandal et al. (2003) While it shares similarities with many of these techniques (albeit from a different perspective), detailed analysis needs to be carried out to understand the precise distinctions between the different approaches.

II Thermofield Dynamics

At non-zero temperatures, the expectation value of an operator $A$ (i.e. its thermal average) must be evaluated as an ensemble average

[TABLE]

where $\beta$ is the inverse temperature, $H$ is the Hamiltonian, $\{\lvert m\rangle\}$ forms a complete orthonormal basis and the partition function $\mathcal{Z}$ is defined as the trace of $\hat{\rho}$, the density operator

[TABLE]

One can choose to work within the grand canonical or the canonical ensemble by appropriately defining the Hamiltonian $H$ with or without a chemical potential term. In this paper, unless mentioned otherwise, we shall work in the grand canonical ensemble.

The central idea in TFD is to express the ensemble average of an operator $A$ in Eq. 2 as an expectation value over just one state $\lvert\Psi(\beta)\rangle$, known as the thermofield double state, thermal vacuum or just thermal state:

[TABLE]

This idea is realized by introducing a copy of the original Hilbert space $\mathcal{H}$, known as the auxiliary or tilde-conjugate space $\tilde{\mathcal{H}}$, such that:

For every state $\lvert\psi\rangle$ in $\mathcal{H}$, there is a copy $\lvert\tilde{\psi}\rangle$ in $\tilde{\mathcal{H}}$ and likewise for operators. 2. 2.

The tilde operators obey similar (anti-) commutation rules as their un-tilde counterparts. For instance, for a bosonic (fermionic) spin-orbital $k$, we have the extended set of field operators

[TABLE]

following the (anti-) commutation rules given by

[TABLE]

where the convention for commutator / anti-commutator is defined as

[TABLE]

such that $\eta=-1$ in Eq. 5 produces the anticommutation rules for fermions and $\eta=+1$ produces the commutation rules for bosons. 3. 3.

A tilde conjugation operation transforms operators between $\mathcal{H}$ and $\tilde{\mathcal{H}}$ with the following general rules:

[TABLE]

where $\alpha^{\star},\delta^{\star}$ are complex conjugates of $\alpha,\delta$ respectively, and $\eta=\pm 1$ for bosons / fermions in Eq. 7. With these conjugation rules, a Hamiltonian for the tilde system can be defined, generally denoted by $\tilde{H}$. 4. 4.

The time-dependent Schrödinger equation in $\tilde{\mathcal{H}}$ becomes ($\hbar=1$)

[TABLE] 5. 5.

Operators in the physical space do not act on states in the tilde space, and vice versa.

The conjugation and doubling of the Hilbert space has well-justified connections with Hopf algebra.Celeghini et al. (1998); Floquet, Trindade, and Vianna (2017)

In the expanded space, the thermal state in Eq. 4 can be expressed as

[TABLE]

where $\{|m\rangle\}$ forms an orthonormal basis in $\mathcal{H}$ and $H$ is the Hamiltonian of the system (and so does not act on states in the tilde space). The norm of the thermal state gives the partition function

[TABLE]

Note that at infinite temperature or $\beta=0$, the thermal state is merely given by

[TABLE]

Since the Hamiltonian is no longer relevant, the state $|\mathbb{I}\rangle$ only depends on the structure of the Hilbert space and can be computed exactly. Consequently, $|\mathbb{I}\rangle$ is analogous to the identity operator and is independent of the choice of basis $\{|m\rangle\}$.

The $|\tilde{m}\rangle$ states in Eq. 11 perform the role of a tracer. That is, for some operator $\hat{A}$ which acts only on the physical states,

[TABLE]

where we have used the shorthand notation $|m,\tilde{m}\rangle=|m\rangle\otimes|\tilde{m}\rangle$.

Conventional applications of TFD have been mostly centered around the study of dynamics of quantum systems at finite temperatures which is governed by a new Hamiltonian

[TABLE]

which drives the real-time dynamics of the thermal state through the Schrödinger equation

[TABLE]

This real-time Schrödinger equation does not provide any prescription to compute the thermal state $|\Psi(\beta)\rangle$ and hence cannot be used to compute the partition function and other equilibrium thermal properties for the original system. Instead, we realize that $|\Psi(\beta)\rangle$, while not an eigenstate of the physical Hamiltonian $H$, obeys an imaginary time Schrödinger equation, given by

[TABLE]

which allows us to construct the thermal state at any $\beta$ by integrating Eq. 17 starting from temperature $\beta_{0}$ where the thermal state $|\Psi(\beta_{0})\rangle$ is known exactly. Similarly, if an explicit dependence on the chemical potential $\mu$ (or rather $\alpha=\beta\mu$) is considered, an equation for the evolution in $\alpha$ can be established, i.e.

[TABLE]

In what follows, we explicitly describe the $\beta$ evolution and provide a prescription to extend the formalism to the $\alpha$ evolution in accordance with Eq. 19. For the sake of simplicity, we will refer to $\alpha$ as the chemical potential.

II.1 Thermal Mean-Field Theory

While Eq. 11 is formally exact, it is in practice impossible to determine the thermal state exactly except for the simplest model systems (it is, after all, a highly non-trivial problem to compute a single state and Eq. 11 requires us to compute every state). Just as in ground state calculations, practical applications require a systematic way of approximating $\left\lvert\Psi(\beta)\right\rangle$.

The most elementary approximation one can invoke is the mean-field approach, wherein the Hamiltonian $H$ is approximated by a one-body mean-field Hamiltonian (e.g. the Fock operator), and correspondingly, the mean-field thermal vacuum takes the form

[TABLE]

where $H_{0}$ is the mean-field Hamiltonian. This mean-field state satisfies the imaginary time Schrödinger equation not for $H$ but for $H_{0}$:

[TABLE]

If we choose the states $\{\lvert m\rangle\}$ to be the eigenstates of the mean-field Hamiltonian $H_{0}$, which can generally be computed without much computational effort, the thermal vacuum takes the form

[TABLE]

The norm of $\left\lvert 0(\beta)\right\rangle$ gives the mean-field partition function

[TABLE]

Any quantity of interest, such as the energy, can then be approximately evaluated as an expectation value over $\left\lvert 0(\beta)\right\rangle$.

For fermions, working within the grand canonical formulation, the thermal state can also be expressed in terms of the single-particle Fock states:

[TABLE]

where $|0\rangle_{p}$ and $|1\rangle_{p}$ respectively denote that the single-particle level $p$ is empty or occupied. If the levels $p$ are chosen to be the eigenstates of $H_{0}$, a simple thermal Bogoliubov transformation can be defined to construct the thermal field operators $\{a_{p}(\beta),a_{p}^{\dagger}(\beta),\tilde{a}_{p}(\beta),\tilde{a}_{p}^{\dagger}(\beta)\}$ that create or annihilate quasi-particle excitations on to the thermal vacuum, i.e.

[TABLE]

such that

[TABLE]

It is easy to show that $w_{p}=\sqrt{1-f_{F}(\epsilon_{p})}$ and $z_{p}=\sqrt{f_{F}(\epsilon_{p})}$, where $f_{F}$ is the Fermi-Dirac distribution function. Here and in the following, we have dropped the explicit $\beta$-dependence in the operators $a_{p}$, $a_{p}^{\dagger}$, and their tilde counterparts.

While a similar transformation can be defined for bosons, henceforth we shall confine ourselves to the many-electron problem in this paper. Furthermore, while such a transformation can also be built within the canonical ensemble, it is generally non-linear and introduces a new algebra for the thermal operators (see, e.g., the transformation for $SU(2)$ spin operators in Ref. Suzuki, 1986).

The Bogoliubov transformation in Eq. 25 allows us to form a physical intuition about the thermal operators: in the low $T$ (or high $\beta$) limit, the annihilation of a particle in the original space $\mathcal{H}$ is equivalent to creating a particle in the tilde-space $\tilde{\mathcal{H}}$ and vice-versa, whereas in the high $T$ limit, these operators are completely mixed.

A general 2-body Hamiltonian, written in terms of the spin-orbital creation and annihilation operators as

[TABLE]

can be expressed in terms of the thermal operators and takes the form

[TABLE]

Notice that the foregoing expression for the Hamiltonian is normal-ordered with respect to the thermal vacuum. The ordering we follow here is $a^{\dagger}\rightarrow\tilde{a}^{\dagger}\rightarrow\tilde{a}\rightarrow a$. The thermal average of $H$ at the mean-field level is simply found by extracting the scalar component from its normal ordered expression,

[TABLE]

which recovers the standard thermal Hartree-FockMermin (1963); Sokoloff (1967) expression for the energy. One can find the appropriate one-electron basis or molecular orbitals by variationally minimizing the appropriate free energy i.e. the Helmholtz free energy in the canonical ensemble and the grand potential in the grand canonical ensemble.

III Framework for Correlated Methods

We have noted that a practical wave function-based framework for the study of the thermal properties of electronic systems is highly desirable, and have seen that thermofield dynamics allows such a framework at the mean-field level. Here, we wish to include correlation atop thermal mean-field.

Correlated methods frequently use a mean-field reference as the starting point, and we wish to do so here as well, but we face an additional choice which we wish to explore. Recall that the thermal state $|\Psi(\beta)\rangle$ is obtained not from an eigenvalue problem but from an imaginary-time Schrödinger equation. We can choose as our reference the mean-field thermal state corresponding to a fixed temperature $\beta_{0}$, or we can instead use as our reference the mean-field thermal state corresponding to the temperature of interest. By analogy with similar frameworks for coordinates in fluid dynamics as well as general relativity, we call the former approach a fixed-reference formulation and the latter a covariant formulation. Mathematically, the exact thermal vacuum $|\Psi(\beta)\rangle$ is represented as

[TABLE]

in the fixed-reference approach, where $\hat{\Gamma}$ is a wave operator which builds correlation on the reference, and as

[TABLE]

in the covariant approach.

On the one hand, the covariant approach would seem to be more sensible, as less is demanded of the wave operator $\hat{\Gamma}$; on the other hand, the fixed-reference approach has the advantage that the quasiparticle creation and annihilation operators given by the thermal Bogoliubov transformation of Eq. 25 are not themselves temperature-dependent, which considerably simplifies the formulation of correlated methods. In principle any inverse temperature $\beta_{0}$ can be used in the fixed-reference case, but in practice the most convenient choice is $\beta_{0}=0$ for which the mean-field thermal state is exact and the wave operator $\hat{\Gamma}(0,0)$ is simply the identity operator.

Beyond deciding between the fixed-reference and covariant cases, we have a second decision to make. The mathematical structure of the thermal averages expressed in the TFD formalism exhibits a great degree of flexibility in the way in which the bra and the ket thermal states can be split. This is a direct consequence of the cyclic property of the trace involved, or equivalently the freedom in the choice of the path integral representation of TFD.Das and Kalauni (2016)

Consider the thermal average of some physical quantity $A$ which is written in TFD as

[TABLE]

The numerator in Eq. 32 represents the trace of the operator $A$ along with the density operator $\hat{\rho}$. The cyclic property of the trace allows us to rewrite Eq. 32 as

[TABLE]

where we are free to choose $0\leq\sigma\leq 1$.Matsumoto et al. (1983); Das and Kalauni (2016) The cyclic property of the trace is applicable to an expectation value here because $A$ acts only on the states in $\mathcal{H}$. Note that for $\sigma\neq 1/2$, we would require not one but two thermal states (one each for the bra and the ket). Two values of $\sigma$ are particularly useful: $\sigma=1/2$ for which we require only a single thermal state, and $\sigma=1$ for which one of the two requisite thermal states is simply $\langle\mathbb{I}|$. When this $\sigma$-dependence is needed to properly interpret the equations, we will hereafter indicate it by the inclusion of a subscript which for the sake of economy we will suppress when it is not needed.

We close this section by giving reasonably explicit recipes for the thermal state in the fixed-reference and covariant cases, assuming that we evolve in $\beta$ from $\beta=0$. For the fixed-reference case, with $\beta_{0}=0$, we have

[TABLE]

where we satisfy

[TABLE]

and where $\hat{1}$ is the identity operator. In contrast, for the covariant case we have

[TABLE]

and we satisfy

[TABLE]

Inserting Eq. 37a into Eq. 37b gives us an alternative expression:

[TABLE]

where the $\beta$-dependent effective Hamiltonian is

[TABLE]

and $V=H-H_{0}$ is the correction to the mean-field Hamiltonian $H_{0}$. In Sec. IV, we will show detailed formulations for these two approaches through the example of thermal configuration interaction (CI) theory.

III.1 Chemical potential

The framework described so far in Sec. III is common to both the canonical and the grand canonical ensemble. However, as we have noted, working in the grand canonical ensemble is more convenient. Consequently, the thermal state, and therefore all physical quantities, are a function of two state variables: the inverse temperature $\beta$ and chemical potential $\alpha$. While we have focussed our discussion on the evolution of the thermal state with respect to $\beta$, we should also discuss its evolution with respect to $\alpha$.

Just as with the $\beta$-evolution, we can evolve in $\alpha$ using either a fixed-reference or a covariant approach; further, we also have the freedom in splitting the contribution from $\alpha$ (cf. Eq. 18) between the bra and the ket. In principle, one can independently choose between the fixed-reference and covariant approaches for the $\alpha$- and $\beta$-evolutions, and can have different $\sigma$’s for the $\alpha N$ and $\beta H$ terms of the $e^{\alpha N-\beta H}$. In practice, we prefer to either use the fixed-reference approach for both $\alpha$- and $\beta$-evolutions or the covariant approach for both, and we use the same $\sigma$ as well (so that we partition $\alpha N-\beta H$ rather than the two terms separately).

At $\beta=0$, we pick some chemical potential $\alpha_{0}$ as the initial value. Thereafter, for the fixed-reference case, we have

[TABLE]

which, in addition to Eq. 35 (which now uses $|\mathbb{J}\rangle$ instead of $|\mathbb{I}\rangle$ as the reference), satisfies

[TABLE]

On the other hand, for the covariant formulation, we have

[TABLE]

which, in addition to Eq. 37 (where the mean-field reference $|0_{\sigma}(\alpha,\beta)\rangle$ now also depends on $\alpha$), satisfies

[TABLE]

Because we wish to ultimately make contact with the canonical ensemble at zero temperature, we select the initial chemical potential $\alpha_{0}$ so that the reference state $|\mathbb{J}\rangle$ at $\beta=0$ has the desired average particle number, and for each $\beta$ we evolve in $\alpha$ so that $\langle N\rangle$ is constant.

IV Thermal CI Theory

Here, we specialize both the fixed-reference and the covariant versions of our general theory to the case of thermal CI. As will be clear, these approaches are quite similar to time-dependent CIKlamroth (2003); Huber and Klamroth (2005) at zero temperature. While CI is not size-extensive, a property which is highly desirableNooijen, Shamasundar, and Mukherjee (2005) in the study of many-electron systems, we use it as an example to introduce thermofield-based quantum chemistry methods because of its simplicity. The construction of more sophisticated methods, such as thermal coupled cluster theory, is mathematically more involved and its details will be presented in a follow-up article.Harsha, Henderson, and Scuseria Nevertheless, we here present preliminary numerical results of thermal CC with singles and doubles (CCSD) for the purpose of comparison with CI singles and doubles (CISD).

IV.1 Fixed-reference Thermal CI

In this first approach, we express thermal averages with $\sigma=1$, so that we use the asymmetric expectation value

[TABLE]

where the state $|\psi(\alpha,\beta)\rangle$ is governed by the imaginary-time Schrödinger equation

[TABLE]

and is written as

[TABLE]

where all the $\alpha$- and $\beta$-dependence is carried by the expansion coefficients. Notice that the wave operator in Eq. 46 is written in intermediate normalization with an exponential normalization constant. Since the norm of the thermal state represents the partition function which decays exponentially, an exponential parametrization provides numerical stability. The ($\alpha$- and $\beta$-independent) quasiparticle operators $b$, $\tilde{b}^{\dagger}$, and so on are defined by Eq. 25 with

[TABLE]

where $n_{f}$ is the filling fraction at $\beta=0$ for a given basis set. It is interesting to note that the wave operator $\hat{\Gamma}$ in Eq. 46 is composed of terms that contain equal number of tilde and non-tilde quasiparticle creation operators. This is because, by virtue of the thermal Bogoliubov transformation in Eq. 25, the difference in the total number of tilde and non-tilde quasiparticles is a symmetry of the Hamiltonian, i.e.

[TABLE]

and therefore, only terms which contain an equal number of tilde and non-tilde quasiparticle operators would have a non-trivial contribution to the thermal state.

Substituting this ansatz into Eq. 45, we obtain the following set of working equations for $\beta$-evolution:

[TABLE]

while for $\alpha$-evolution, we get

[TABLE]

and similarly for higher order terms, with the infinite-temperature partition function $\mathcal{Z}_{\mathbb{J}}$ being given by

[TABLE]

We can integrate these equations from $(\alpha,\beta)=(0,0)$, with the initial conditions

[TABLE]

to obtain the CI coefficients at any given thermodynamic coordinate $(\alpha,\beta)$.

IV.2 Covariant CI

In the covariant formulation, we wish to work with an $\alpha$- and $\beta-$dependent reference and the expectation values are best described with the $\sigma=1/2$ formulation of TFD and take the form of a symmetric expectation value

[TABLE]

where the thermal state $|\Psi\rangle$ is the same as the state described in Eq. 11 and is governed by the imaginary time Schrödinger equation described in Eq. 17 and 19..

We parametrize the thermal state as

[TABLE]

where unlike in Eq. 46 both the coefficients and the field operators are $\alpha$- and $\beta$-dependent. Because the operators $a_{p}$, $\tilde{a}_{p}^{\dagger}$, and so forth are $\alpha$- and $\beta$-dependent, they have non-trivial $\alpha$ and $\beta$ derivatives. For example, from the thermal Bogoliubov transformation of Eq. 25, we see that

[TABLE]

and so on.

Substituting the wave function ansatz of Eq. 52 into Eq. 17, we find

[TABLE]

and for Eq. 19, we get

[TABLE]

and so forth, where $\partial_{\mathrm{op}}\hat{S}/\partial x$ (with $x=\alpha,\beta$) denotes the derivative of only the operator part of the wave operator $\hat{S}$, and

[TABLE]

while $\mathcal{Z}_{0}$ is the mean-field partition function. As with the fixed-reference case, we can integrate Eq. 54 and 55 from $(\alpha,\beta)=(0,0)$ with the initial conditions

[TABLE]

Our primitive implementation of both the fixed-reference and the covariant formulations of thermal CI have been made publicly accessible through a GitHub repositoryHarsha where we make use of drudgeZhao and Scuseria (computer algebra system) to generate equations and codes.

V Results and Conclusions

Just as for ground-state CI, the theories described above can be truncated at various levels of excitations to yield results with respective levels of accuracy. Here, we truncate after single (CIS) or double (CISD) excitations to study the temperature-dependence of the internal energy of the Hubbard modelHubbard J. and Flowers Brian Hilton (1963) with periodic boundary condition in the grand canonical ensemble. The Hamiltonian for the Hubbard model 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 $\hat{n}_{p,\sigma}=c_{p,\sigma}^{\dagger}\,c_{p,\sigma}$ is the number operator for lattice site $p$ and spin $\sigma$. The correlation strength in the Hubbard model is generally characterized by the ratio $U/t$.

Recall that our formulation has been given in the grand canonical ensemble. Our current implementation treats the inverse temperature $\beta$ and the chemical potential $\alpha$ as independent variables and evolves the differential equations of Eq. 17 and Eq. 19 in the $\alpha$-$\beta$ plane in such a way as to maintain the desired average number of particles. We have omitted explicit orbital optimization of thermal HF at each $\alpha$, $\beta$ value, and the orbitals used in thermal HF and covariant thermal CI are all the same. We choose to work in the basis in which the Fock operator is diagonal. This particularly simplifies the thermal Bogoliubov transformation to the form described in Eq. 25.

We begin with the simple two-site model with (on average) two electrons and $U/t=1$, for which both the errors in the internal energy and its value, calculated using the covariant CI method, are displayed in Fig. 1. While thermal mean-field is exact at $\beta=0$, we quickly see large deviations from the exact result. Thermal CIS does not introduce any significant improvement over the Hartree-Fock internal energy, however the thermal CISD model is fairly accurate. Thermal CCSD further reduces the errors for moderate to large values of $\beta$, although it is not as accurate near $\beta=0$. This is likely due to the fact that in our current implementation, thermal CC evaluates expectation values as

[TABLE]

while it should, in principle, include

[TABLE]

where by $\langle\Psi_{\textrm{L}}|$, we refer to a more elaborate, correlated bra state. On the other hand, thermal CI, by definition, computes thermal averages as a symmetric expectation value, with both the bra and the ket being correlated wavefunctions. In the zero-temperature limit ($\beta\to\infty$) thermal HF reduces to the canonical zero-temperature ground-state restricted HF (RHF), and thermal CI reduces to the canonical zero-temperature ground-state restricted CI; for this reason, the zero-temperature limit of thermal CISD is exact (because CISD is exact for two-particle systems). The error for intermediate values of $\beta$ with thermal CISD is because we work in the grand canonical ensemble for which the exact thermal CI requires also triple and quadruple excitations.

Figure 2 shows results for the six-site model at $U/t=2$ and average half filling; again we use the covariant formulation. The story is similar to the simpler two-site model: there are large errors at thermal HF and CIS level which are substantially reduced by thermal CI and CC. For large $\beta$, thermal CISD reduces to the standard zero-temperature canonical ensemble CISD.

As we have seen, covariant CI works well over the entire temperature range, particularly once we include double excitations. In contrast, the fixed-reference truncated CI does not perform well as we go away from $\beta=0$, and soon becomes inferior to even thermal HF. This is because as we evolve away from the initial $\beta$, the mean-field contribution to the internal energy (or any other property of interest) remains fixed and all the temperature dependence and correlation effects need to be taken care of by the correlator or the wave operator. Analogous to the reference-dependent behaviour of zero temperature ground-state CI, such a method would perform poorly compared to one with a more optimized mean-field reference at each $\beta$-point. This is emphasized in Fig. 3 which shows results for the six-site model at $U/t=2$ with an average of two electrons. Lastly, we note that the performance of any thermal wavefunction method is more or less equivalent to its ground-state symmetry adapted counterpart. Accordingly, the example we implement here, i.e., CI, does not perform equally well for stronger correlations, a feature shared by most ground-state methods based on the RHF reference (see, e.g., Refs. Henderson et al., 2014 and Degroote et al., 2016). Hence, we have included results only for weakly correlated Hubbard systems.

Conclusions

By integrating an imaginary time Schrödinger equation for the thermal state, thermofield theory provides an honest wave function method from which the thermal density matrix and any thermal property may be computed. This technique thus allows for the straightforward extension of ground-state methods to the computation of properties at non-zero temperatures. Moreover, the doubling of the Hilbert space due to the introduction of the fictitious tilde space does not affect the overall computational scaling of any method. In comparison with the ground state method, the cost increases only by a factor proportional to the number of points in the $\beta$-$\mu$ evolution. These features make this framework ideal for application to many-body quantum systems in physics and chemistry and provide a deterministic alternative to the more common stochastic Monte-Carlo methods. Our early results show that even thermal CI considerably improves upon thermal mean-field theory, provided at least that one uses the covariant approach for the evolution of the thermal state. While the fixed-reference formulation is significantly easier to derive and implement, its performance is simply inadequate and we do not recommend its use.

We close by emphasizing that it was not our intention here to delve into a detailed comparison between our TFD-based approach and other finite-temperature methods. Such a comparison would be of value and, as with the details of our TFD-based coupled cluster methodology, will appear in due time. But while the story is incomplete, this early work provides both the desire and the means to generalize more sophisticated wavefunction techniques to the case of finite temperature.

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