Spin-Projected Generalized Hartree-Fock as a Polynomial of Particle-Hole Excitations

TL;DR

This paper demonstrates that spin-projected generalized Hartree-Fock wave functions can be expressed as polynomials of excitation operators, enabling potential integration with coupled cluster methods to better describe both weak and strong electron correlations.

Contribution

It extends previous work by showing that spin-projected generalized Hartree-Fock wave functions are polynomials of excitation operators, providing a foundation for combined methods with coupled cluster theory.

Findings

Wave functions are expressible as low-order excitation polynomials.

Preliminary applications to the Hubbard model show promise.

Method could improve accuracy across correlation regimes.

Abstract

The past several years have seen renewed interest in the use of symmetry-projected Hartree-Fock for the description of strong correlations. Unfortunately, these symmetry-projected mean-field methods do not adequately account for dynamic correlation. Presumably, this shortcoming could be addressed if one could combine symmetry-projected Hartree-Fock with a many-body method such as coupled cluster theory, but this is by no means straightforward because the two techniques are formulated in very different ways. However, we have recently shown that the singlet $S^2$-projected unrestricted Hartree-Fock wave function can in fact be written in a coupled cluster-like wave function: that is, the spin-projected unrestricted Hartree-Fock wave function can be written as a polynomial of a double-excitation operator acting on some closed-shell reference determinant. Here, we extend this result and show that the spin-projected generalized Hartree-Fock wave function (which has both $S^2$ and $S_z$ projection) is likewise a polynomial of low-order excitation operators acting on a closed-shell determinant, and provide a closed-form expression for the resulting polynomial coefficients. We include a few preliminary applications of the combination of this spin-projected Hartree-Fock and coupled cluster theory to the Hubbard Hamiltonian, and comment on generalizations of the methodology. Results here are not for production level, but a similarity transformed theory that combines the two offers the promise of being accurate for both weak and strong correlation, and particularly may offer significant improvements in the intermediate correlation regime where neither projected Hartree-Fock nor coupled cluster is particularly accurate.

Tables

Table 1. Table 1: Coefficients λ 𝜆 \lambda for the SGHF wave function and for the analogous coupled cluster wave function. The excitation level is N = 2 ​ i + 3 ​ j + 4 ​ k 𝑁 2 𝑖 3 𝑗 4 𝑘 N=2i+3j+4k . The coefficients in coupled cluster theory are simply 1 / ( i ! ​ j ! ​ k ! ) 1 𝑖 𝑗 𝑘 1/(i!\,j!\,k!) , and we call them λ i ​ j ​ k CC superscript subscript 𝜆 𝑖 𝑗 𝑘 CC \lambda_{ijk}^{\mathrm{CC}} .

$N$	$i$	$j$	$k$	${\lambda }_{ijk}$	${\lambda }_{ijk}^{\mathrm{CC}}$
2	1	0	0	1	1
3	0	1	0	1	1
4	0	0	1	1	1
	2	0	0	3/10	1/2
5	1	1	0	3/5	1
6	1	0	1	3/7	1
	0	2	0	6/35	1/2
	3	0	0	3/70	1/6
7	0	1	1	3/14	1
	2	1	0	9/70	1/2
8	0	0	2	5/84	1/2
	2	0	1	1/14	1/2
	1	2	0	1/14	1/2
	4	0	0	1/280	1/24

Equations

$$|\mathrm{GHF}⟩={\mathrm{e}}^{T_1+U_0+V_{+}+W_{-}}|\mathrm{RHF}⟩$$

$$T_1$$

$$U_0$$

$$V_{+}$$

$$W_{-}$$

$$E_{ai}$$

$$S_{ai}^0$$

$$S_{ai}^{+}$$

$$S_{ai}^{-}$$

$$|\mathrm{SGHF}⟩$$

$$=\frac{1}{2}{P}_{s_z=0}{\int }_0^{\pi }\mathrm{d}\beta \sin (\beta ){\mathrm{e}}^{\mathrm{i}\beta S_y}{P}_{s_z=0}|\mathrm{GHF}⟩,$$

$$R(\Omega )$$

$$P_{s_z=0}|\mathrm{GHF}⟩$$

$$={\mathrm{e}}^{U_0}{P}_{s_z=0}\sum _{m,n}\frac{1}{m!}\frac{1}{n!}{V}_{+}^m{W}_{-}^n|\mathrm{RHF}⟩$$

$$={\mathrm{e}}^{U_0}\sum _n\frac{1}{(n!{)}^2}{\left(V_{+}{W}_{-}\right)}^n|\mathrm{RHF}⟩,$$

$${\mathrm{e}}^{\mathrm{i}\beta S_y}\mathcal{O}|\mathrm{RHF}⟩$$

$$=\tilde{\mathcal{O}}|\mathrm{RHF}⟩,$$

$$\tilde{\mathcal{O}}={\mathrm{e}}^{\mathrm{i}\beta S_y}\mathcal{O}{\mathrm{e}}^{-\mathrm{i}\beta S_y}.$$

$$\tilde{U}$$

$$\tilde{V}$$

$$\tilde{W}$$

$${\Theta }_k=\sum {\theta }_i^a{S}_{ai}^k$$

$$|\mathrm{SGHF}⟩=\sum _N{C}_N|\mathrm{RHF}⟩$$

$$C_N$$

$$\times P_{s_z=0}{\int }_0^{\pi }\mathrm{d}\beta \sin (\beta ){\tilde{U}}^{N-2M}{\tilde{V}}^M{\tilde{W}}^M$$

$$\frac{{\tilde{U}}^{N-2M}}{(N-2M)!}$$

$$\frac{{\tilde{V}}^M}{M!}$$

$$\frac{{\tilde{W}}^M}{M!}$$

$$C_N$$

$$\times {\delta }_{i+k+m,j+l+n}\frac{U_0^{N-2M-i-j}}{(N-2M-i-j)!}\frac{U_{+}^i}{i!}\frac{U_{-}^j}{j!}\frac{V_0^{M-k-l}}{(M-k-l)!}\frac{V_{+}^k}{k!}\frac{V_{-}^l}{l!}\frac{W_0^{M-m-n}}{(M-m-n)!}\frac{W_{+}^m}{m!}\frac{W_{-}^n}{n!},$$

$$\mathcal{I}(p,q,r)={\int }_{-1}^1\mathrm{d}x{x}^p(1+x{)}^q(1-x{)}^r.$$

$$\sin (\beta {)}^{2M+i+j-k-l-m-n}=\sin (\beta {)}^{2(M+i-l-n)}$$

$${\left[\sin (\beta {)}^2\right]}^{M+i-l-n}$$

$$\times {\left[1-\cos (\beta )\right]}^{M+i-l-n}.$$

$$C_N=\sum _{2i+3j+4k=N}{\lambda }_{ijk}{C}_2^i{C}_3^j{K}_4^k;$$

$$|\mathrm{SGHF}⟩=\sum _{ijk}{\lambda }_{ijk}{C}_2^i{C}_3^j{K}_4^k|\mathrm{RHF}⟩$$

$$C_2$$

$$C_3$$

$$+V_0\left(W_{+}{U}_{-}-W_{-}{U}_{+}\right)$$

$$\displaystyle\quad+W_{0}\,\left(U_{+}\,V_{-}-U_{-}\,V_{+}\right)\big{]},$$

$$K_4$$

Full text

I Introduction

Symmetry-projected Hartree-Fock Löwdin (1955); Ring and Schuck (1980); Blaizot and Ripka (1985) has been used for many years in nuclear physics as a tool for the description of strong correlation (see, e.g., Ref. Schmid, 2004) and for the past several years spin-projected Hartree-Fock Jiménez-Hoyos et al. (2012) has seen increasing use in the computational quantum chemistry community as well.Rivero et al. (2013a, b); Jiménez-Hoyos et al. (2013); Stein et al. (2014); Jiménez-Hoyos et al. (2014); Rodríguez-Guzmán et al. (2014) The idea is beguilingly simple. In the presence of strong correlations, the Hartree-Fock method predicts its own failure by breaking some or all of the symmetries of the Hamiltonian. Thus, for example, while restricted Hartree-Fock (RHF) respects both $S^{2}$ and $S_{z}$ spin symmetries, the unrestricted Hartree-Fock (UHF) wave function allows $S^{2}$ spin symmetry to break, the generalized Hartree-Fock (GHF) wave function allows both $S^{2}$ and $S_{z}$ spin symmetries to break, and the Bardeen-Cooper-Schrieffer (BCS) and Hartree-Fock-Bogoliubov (HFB) wave functions allow particle number symmetry to break (and, in the case of HFB, may also allow for spin symmetry breaking). By breaking symmetry, Hartree-Fock gains variational flexibility and provides a better energy at the cost of a perhaps less physical wave function.

One could, however, retain only the component of the broken-symmetry determinant which has good symmetries, i.e. one could project the broken-symmetry Hartree-Fock determinant back onto symmetry eigenstates. The resulting wave function is no longer a single determinant, though it can be expressed as a linear combination of non-orthogonal broken-symmetry determinants. This projected Hartree-Fock (PHF) wave function has the convenient properties that it retains the energetic benefits of broken-symmetry Hartree-Fock but has good symmetries; it is also parameterized in terms of a single determinant, which is computationally convenient. In practice a variation after projection approach is preferred, in which the determinant is chosen to variationally optimize the projected energy (as opposed to a projection after variation approach, in which the determinant is chosen to minimize the mean-field energy and is subsequently projected). An important benefit of this variation after projection technique is that the underlying determinant always breaks symmetry, resulting in smooth potential energy curves.

While PHF has had some significant successes for the description of strong correlation, it is much less successful at capturing dynamic correlation. Presumably, the fact that the wave function is based on an underlying single determinant is too severe a limitation. Adding dynamic correlation to PHF has been a long-standing goal, and one which is beginning to come to fruition.Tshuchimochi and Van Voorhis (2014); Tsuchimochi and Ten-no (2016a, b); Wahlen-Strothman et al. (2016); Hermes et al. (2017); Qiu et al. (2017)

It is worth taking a moment to understand why adding dynamic correlation to PHF is so challenging. The basic problem is that the PHF wave function is not readily expressible as a (hopefully short) combination of single determinants in some form of active space. On the one hand, this is a real strength of PHF: one need not define an active space at all (but can if one wishes). On the other hand, the complicated nature of the PHF wave function makes it challenging to use PHF in combination with traditional multi-reference perturbation theory or coupled-cluster theory.

Fortunately, an alternative exists.Piecuch et al. (1996); Degroote et al. (2016); Qiu et al. (2016, 2017) One can parameterize the broken-symmetry determinant as a Thouless transformationThouless (1960) of a symmetry-adapted reference determinant; this permits the symmetry projection to be done analytically, resulting in an expression for the PHF wave function in terms of particle-hole excitations out of a reference. Written this way, one can hope to combine PHF with more traditional wave function techniques which use a similar form of the wave function.

There is, of course, no such thing as a free lunch. In this case, the important limitation is that different symmetries being projected yield different expressions for the PHF wave function. Moreover, different quantum numbers of the same symmetry may do likewise. So far as we are aware, only the particle-hole representations of the singlet spin-projected UHF (SUHF) wave function,Piecuch et al. (1996); Qiu et al. (2017) the parity projected wave function in the Lipkin Hamiltonian,Wahlen-Strothman et al. (2016) and the number projected HFB wave functionRing and Schuck (1980); Dukelsky et al. (2003); Degroote et al. (2016) are known, at least among symmetries relevant for molecular systems. We should note that the literature has focused on the number projected BCS (PBCS) case, but the number projection of HFB would be done in precisely the same way.

This manuscript adds singlet spin-projected generalized Hartree-Fock (SGHF) to the list. By allowing all spin symmetries to break, GHF has the attractive property that it is always size consistent and size extensiveJiménez-Hoyos et al. (2011) (though this is not true of SGHF or indeed of any PHF state, the extensive component of which is degenerate with the underlying broken-symmetry determinant). Although GHF solutions are, on the whole, rare, SGHF is almost always different from (and variationally superior to) SUHF. Presumably, by being variationally superior to SUHF, SGHF lends itself that little bit more readily to combination with methods which seek to add dynamic correlation. To this end, we show a few preliminary calculations which combine SGHF with variational coupled cluster theory and compare to the combination of SUHF with variational coupled cluster. We show results only for the Hubbard Hamiltonian,Hubbard (1963) where tuning the strength of the interaction allows us to proceed smoothly from weak to strong correlation and where our previous results Qiu et al. (2017) show that the combination of SUHF with variational coupled cluster benefits enormously from the inclusion of triple excitations which, as we shall see, are built into the SGHF wave function already. We will also discuss the combination of spin projection and number projection, to indicate how it might be done. We emphasize, however, that the main point of this paper is to introduce the particle-hole representation of the SGHF wave function.

II The SGHF Wave Function

Here, we will derive an expression for the SGHF wave function in terms of particle-hole excitations out of some RHF reference determinant. We will follow the basic scheme explained in Ref. Qiu et al., 2017. The derivation is a little tedious, but our main result is contained in Eqns. 18 - 20.

We will start by writing a representation of the GHF determinant as a particle-hole Thouless transformation out of the RHF reference:

[TABLE]

where we have separated the Thouless transformation into different components. Specifically, we have defined

[TABLE]

where the excitation operators are

[TABLE]

and where $i$ indexes spatial orbitals occupied in $|\mathrm{RHF}\rangle$ and $a$ indexes spatial orbitals unoccupied in $|\mathrm{RHF}\rangle$. Note that $T_{1}$ preserves the spin quantum numbers $s$ and $s_{z}$, while $U_{0}$ turns a singlet state into a triplet state with $s_{z}=0$, and $V_{+}$ and $W_{-}$ not only turn a singlet into a triplet but also increase or decrease $s_{z}$ by 1, respectively. We use the notation $S_{z}$ to mean the spin operator and $s_{z}$ to denote its eigenvalues.

We should point out two limitations of the particle-hole Thouless transformation we have adopted. First, and more seriously, we can only define it for states $|\mathrm{RHF}\rangle$ and $|\mathrm{GHF}\rangle$ which are not orthogonal, and it becomes numerically ill-conditioned if $|\mathrm{GHF}\rangle$ and $|\mathrm{RHF}\rangle$ have too little overlap. Second, we have written the GHF determinant in a way which does not preserve its normalization. Since the spin projection changes the normalization in any event, this is not a real limitation, but the reader should be aware that neither the GHF nor the SGHF wave functions are properly normalized in this work (though they are normalizable).

On the other hand, the particle-hole Thouless transformation has a few significant advantages. First, the various operators all commute because they are all pure excitation operators. Second, the coefficients of a given excitation level are easily computed.

The SGHF wave function can be extracted from the GHF determinant by finding the component which is rotationally invariant in spin space; this totally symmetric component has $s=s_{z}=0$ and can be computed by Jiménez-Hoyos et al. (2012)

[TABLE]

Here, $P_{s_{z}=0}$ is the projector onto states with $s_{z}=0$. Because the $T_{1}$ component of the Thouless transformation carrying RHF to GHF commutes with the spin rotation operator $R(\Omega)$, it is a spectator in the various projections and integrations of the foregoing equation, and we will omit it for simplicity.

Before we can go further, we need the portion of the GHF wave function which has $s_{z}=0$. This is readily obtained. Noting that the projection operator commutes with $U_{0}$, we can write

[TABLE]

where we have simply used the fact that the projection onto $s_{z}=0$ requires us to increase the $z$-component of spin using $V_{+}$ precisely as many times as we decrease it using $W_{-}$. This is the same logic as used to derive the number projection in number-projected BCS.Degroote et al. (2016)

Now we need to act the rotation operator $\exp(\mathrm{i}\,\beta\,S_{y})$. It is convenient to use

[TABLE]

where

[TABLE]

We have used the fact that the RHF wave function is an eigenfunction of $S_{y}$ with eigenvalue 0. We thus need transformed versions of $U_{0}$, $V_{+}$, and $W_{-}$, which we will respectively call $\tilde{U}$, $\tilde{V}$, and $\tilde{W}$. Writing $S_{y}$ in terms of the fermionic operators $c$ and $c^{\dagger}$ and using the commutator expansion yields

[TABLE]

We have defined additional operators in analogy with $U_{0}$, $V_{+}$, and $W_{-}$. Generically, the operators $\Theta_{k}$ are given by

[TABLE]

where $k$ is [math], $+$, or $-$ and where $\theta$ = $u$, $v$, or $w$ when $\Theta$ is $U$, $V$, or $W$, respectively.

At this point, the SGHF wave function could be written as

[TABLE]

where $N$ indexes the total excitation level. The $N$-fold excitation operator $C_{N}$ is given by

[TABLE]

where the action of the projection operator is understood to mean taking the $\Delta s_{z}=0$ part of the operator string.

One can introduce the definitions of $\tilde{U}$, $\tilde{V}$, and $\tilde{W}$ and use multinomial coefficients to expand the expressions for the rotated operators:

[TABLE]

This means that

[TABLE]

where the Kronecker delta arises from the $s_{z}=0$ projection operator and where the integral $\mathcal{I}(p,q,r)$ is given by

[TABLE]

Here, we have written $x=\cos(\beta)$ and simplified the exponent on $\sin(\beta)$; recalling that one factor of $\sin(\beta)$ is absorbed by the change of variables, for the rest we have

[TABLE]

where we have used the fact that the Kronecker delta enforces that $k+m=j+l+n-i$. Because the exponent on $\sin(\beta)$ is even, we have written

[TABLE]

At this point we resorted to computer algebra Wolfram Research, Inc. , and obtained a rather surprising result. Although we have been unable to rigorously prove it, we find that

[TABLE]

hence

[TABLE]

where

[TABLE]

and

[TABLE]

Equations 17 - 20 constitute the main result of this paper. We have checked that they yield precisely the same coefficients $C_{N}$ as are computed by the explicit form of Eqn. 13 through 30-fold excitations (i.e. $N=30$), though we cannot guarantee that they are correct for still higher excitation levels. It is important to note that where SUHF contains only single and double excitations and powers of themPiecuch et al. (1996); Qiu et al. (2016, 2017), SGHF additionally has triple and quadruple excitations which arise from the additional $S_{z}$ projections and which may have significant energetic consequences (see below).

Table 1 provides numerical values for the coefficients $\lambda_{ijk}$ through octuple excitations, to give an idea of how rapidly the coefficients decay, and compares to the analogous coefficients in coupled cluster theory. In coupled cluster, the wave function is a simple exponential:

[TABLE]

where different levels of theory correspond to different choices of the operator $T$. Since SGHF includes up through quadruple excitations, the natural comparison is coupled cluster with singles, doubles, triples, and quadruples (CCSDTQ). Omitting the singles for brevity, we would have

[TABLE]

where $T_{2}$, $T_{3}$, and $T_{4}$ are double, triple, and quadruple excitation operators. The coefficients in the SGHF wave function clearly decay much more rapidly with increasing excitation level than do the corresponding coefficients in coupled cluster. The same is true for spin projected UHF and for number projected BCS, which have respectively

[TABLE]

where both contain only double excitation operators which, like the SGHF excitation operators, are parameterized in terms of the Thouless transformation taking RHF to the broken-symmetry determinant. This observation for number and $S^{2}$ spin projection forms the basis for our work on what we call attenuated coupled cluster,Gomez et al. (2017) in which we add the SUHF or PBCS polynomial to the coupled cluster exponential, thereby substantially mitigating the failures of traditional coupled cluster theory in the stongly correlated limit. Both the SUHF and the PBCS polynomials are examples of (univariate) hypergeometric functions, as indeed is the exponential of coupled cluster theory; the SGHF polynomial is a multivariate generalization.

It is not terribly obvious that $C_{2}$, $C_{3}$, and $K_{4}$ really are singlet operators (that is, operators expressible in terms of the unitary group generators $E_{ai}$ of Eqn. 3). The form of $\kappa(X,Y)$ is the same as appears in the SUHF polynomial, however, and with a little algebra one finds that

[TABLE]

where $i$ and $j$ index occupied spatial orbitals and $a$ and $b$ index virtual spatial orbitals. Since $\kappa(X,Y)$ is obviously a singlet, it follows that both $C_{2}$ and $K_{4}$ are indeed singlet operators. We do not at present have a clean demonstration that $C_{3}$ actually is a singlet operator, though of course it must be.

A few brief comments are in order. First, we use $K_{4}$ instead of $C_{4}$ in Eqn. 18 simply because in SUHF, $C_{3}=0$ and $C_{4}=3/10\,C_{2}^{2}$; thus, $K_{4}$ is the new quadruply-excited part of the SGHF wave function. Second, both $C_{3}$ and $K_{4}$ are important parts of the SGHF wave function; if they are neglected, one essentially falls back to SUHF. Finally, if one sets $\bm{v}=\bm{w}=\bm{0}$, so that the Thouless transformation forms only a UHF wave function, our result properly reduces to SUHF.

To validate our SGHF expression, we have implemented SGHF in a full configuration interaction code, using the conjugate gradients method to optimize the parameters $t_{i}^{a}$, $u_{i}^{a}$, $v_{i}^{a}$, and $w_{i}^{a}$. Note that $\bm{u}=\bm{v}=\bm{w}=\bm{0}$ is a solution of the SGHF equations which simply corresponds to RHF, and $\bm{v}=\bm{w}=\bm{0}$ is also a solution and corresponds to SUHF. We must therefore provide nonzero initial guesses for $\bm{u}$, $\bm{v}$, and $\bm{w}$. Our results are unfortunately somewhat sensitive to these initial guesses, and this is complicated further by the difficulty in finding GHF solutions for use as an initial guess, since many systems do not spontaneously produce GHF determinants. We have found many GHF solutions in the periodic Hubbard Hamiltonian, but most are orthogonal to the RHF ground state and hence cannot be reached by a non-unitary Thouless transformation of the sort we haved used in our derivation.

Nevertheless, we have tested our polynomial for the 6-site Hubbard model at half-filling by carrying out SGHF calculations with real molecular orbitals (hence real $\bm{t}$, $\bm{u}$, $\bm{v}$ and $\bm{w}$) using the standard technique of integrating over the symmetry coherent statesJiménez-Hoyos et al. (2012) and compared to results using our SGHF polynomial. The energies agree to the convergence tolerance of the respective calculations (in this case, to better than six decimal places). Though the half-filled six-site Hubbard model does not test polynomial coefficients for excitations higher than $C_{6}$, it serves to validate both the coefficients to that order and the procedure we used to derive these coefficients. Note that our polynomial does not assume real wave function parameters and is equally valid for complex Thouless transformations; our code, however, is limited to a real Thouless transformation acting on a real RHF wave function. This means that the underlying broken-symmetry GHF determinant is limited to coplanar spins. Similarly, while we have tested it for the Hubbard Hamiltonian, our derivation is universal and applies to any generalized Hartree-Fock wave function, provided only that the broken-symmetry GHF determinant is not orthogonal to the RHF determinant used as a reference for the Thouless transformation.

II.1 Discussion

It may be illuminating to explain how we obtained the formula for $\lambda_{ijk}$ given in Eqn. 20.

First, then, we obtained the values of the coefficients by finding the explicit form of $C_{N}$ using Eqn. 13 and assuming it had the factorized form of Eqn. 17. We compared unique operator strings on each side of the foregoing equation (e.g. $U_{0}^{4}\,V_{+}\,W_{-}$) and solved for the coefficients $\lambda_{ijk}$ needed to satisfy the equality. Although we have not been able to automate this process, we have had little difficulty in carrying it out.

Inspired by the form of the SUHF wave function, we grouped all powers of $C_{2}$ for a given combination $C_{3}^{j}\,K_{4}^{k}$ and observed that

[TABLE]

accurately predicted all $\lambda_{ijk}$ that we had found. The special functions $\bar{i}_{n}$ are related to modified spherical Bessel functions of the first kind:

[TABLE]

The coefficients $b_{j,k}$ turned out to be

[TABLE]

though this was not initially obvious to us.

What was clear was that for a given power of $C_{3}$, increasing the power of $K_{4}$ by one is always accompanied by increasing the index or indices on the $\bar{i}_{n}$ by two. By the recursion relation of Eqn. 26b, this is tantamount to taking two derivatives, which permits us to rewrite the SGHF wave function as

[TABLE]

where we have introduced

[TABLE]

and where

[TABLE]

While we had not recognized the pattern in the $b_{j,k}$, to our gratified surprise we recognized the coefficients in the power series in $\mathcal{X}$ for a given power of $C_{3}$. Defining a second set of special functions $\bar{I}_{n}$ related to the modified (but not spherical) Bessel functions and satisfying

[TABLE]

we saw that we were apparently able to write the SGHF wave function in the compact and rather suggestive form

[TABLE]

where we have checked the sum on $j$ through $j=4$ and where the polynomials $P_{j}$ are

[TABLE]

Assuming the foregoing result to be true for all $ik$, we derived explicit expressions for the coefficients $\lambda_{ijk}$ for all $ik$ for a given $j$. Having done so for several values of $j$, the final result of Eqn. 20 was obvious.

Let us, then, be clear about what we do and do not mean to say. We do not claim to have proven that the SGHF wave function factorizes to all excitation levels, but it does factorize at all excitation levels we have been able to check. Likewise, though the coefficients of Eqn. 20 are correct in every case we have been able to verify, they were obtained by recognizing patterns in the coefficients for low excitation levels and assuming those patterns hold ad infinitum. This they need not do, and though we think it unlikely, a counterexample to either Eqn. 18 or Eqn. 20 might eventually be found.

III SGHF in Terms of Tensor Invariants

Writing the SGHF Thouless transformation in terms of the operators $U_{0}$, $V_{+}$, and $W_{-}$ simplifies the derivation, but comes at a cost: the three operators are treated differently in the final expression. A more symmetric expression can be derived by considering an alternative parameterization.

Let us first write down an alternative set of triplet single excitation operators to replace $U$, $V$, and $W$:

[TABLE]

where $\eta$ and $\zeta$ run over $\uparrow$ and $\downarrow$ and the matrices $\sigma$ are Pauli matrices. Insisting that $X+Y+Z=U_{0}+V_{+}+W_{-}$ reveals that

[TABLE]

What does this mean for the SGHF operators $C_{2}$, $C_{3}$, and $K_{4}$?

We begin by noting that $\kappa(A,B)$ is bilinear:

[TABLE]

where $\alpha$ and $\beta$ are numbers and $A$, $B$, $C$, and $D$ are operators. It is also symmetric: $\kappa(A,B)=\kappa(B,A)$. Together, these imply

[TABLE]

and hence

[TABLE]

Let us define a rank two tensor operator $\kappa$ with elements

[TABLE]

where $e_{1}=X$, $e_{2}=Y$, and $e_{3}=Z$. Then we have simply

[TABLE]

which is of course an invariant under rotations which mix $X$, $Y$, and $Z$. Similarly, one finds

[TABLE]

which is also an invariant under these rotations.

Writing $C_{3}$ in terms of $X$, $Y$, and $Z$ gives us

[TABLE]

where we have defined $X_{0}$, $X_{+}$, $X_{-}$, and so on in analogy with our other operators. Writing a rank three tensor with elements

[TABLE]

we see that

[TABLE]

which is again rotationally invariant.

Our real point here is that the SGHF wave function can be expressed in a manifestly symmetric form in terms of simple invariants of a rank two Cartesian singlet tensor (with elements $\kappa_{ij}$) and a rank three Cartesian singlet tensor (with elements $\kappa_{ijk}$).

IV Combining SGHF and Variational Coupled Cluster

Our main goal in this work is to present the particle-hole form of the SGHF wave function. Since, however, the reason we are interested in this form is so that we can combine SGHF and coupled cluster theory, we would be remiss if we did not try at least a preliminary investigation of this combination. As in Ref. Qiu et al., 2017 we use a full configuration interaction code to test the combination of SGHF with variational coupled cluster with singles and doublesBartlett and Noga (1988); Van Voorhis and Head-Gordon (2000); Cooper and Knowles (2010); Evangelista (2011) (VCCSD) and with variational coupled cluster with singles, doubles, and triples (VCCSDT). In analogy with Ref. Qiu et al., 2017, we will denote the combination of SGHF and VCCSD as SGVCCSD, and the combination of SGHF with VCCSDT as SGVCCSDT; similarly, the combination of SUHF with VCCSD or with VCCSDT will be denoted respectively by SUVCCSD and SUVCCSDT.

We should emphasize that it is not our intention to use methods like SGVCCSD for production-level calculations. Rather, we just want to see to what extent the method has promise; practical calculations would have to combine SUHF or SGHF with a more traditional similarity-transformed coupled cluster approach, which would lead to polynomial computational scaling rather than the combinatorial scaling forced by the use of variational coupled cluster.

We carry out these tests on the periodic Hubbard Hamiltonian in one dimension. The Hamiltonian models electrons on a lattice, and is given by

[TABLE]

Here, $i$ and $j$ index lattice sites and the notation $\langle ij\rangle$ means that sites $i$ and $j$ are connected in the lattice. We consider only nearest-neighbor hopping ($i=j\pm 1$) and periodic boundary conditions (sites $1$ and $N$ in an $N$-site lattice are also connected). The parameter $U/t$ controls the strength of correlation, with small $U/t$ corresponding to a weakly-correlated system and large $U/t$ corresponding to strong correlation.

We begin with the half-filled 6-site lattice, shown in Fig. 1. The left panel shows the total energy per electron for RHF, VCCSD, and VCCSDT in comparison to the exact result; we omit results from spin-projected methods as they are all essentially superimposable with the exact data. To see how well the spin-projected methods perform, particularly in combination with VCCSD or VCCSDT, we show the error per electron with respect to the exact result in the right panel. It is readily apparent that in this half-filled case, neither VCCSD nor even VCCSDT is of particularly high quality except at small $U/t$. For sufficiently large $U/t$, SUHF is very good. Adding VCCSD to SUHF improves the results significantly in the intermediate range of $U/t$ but has little effect for large $U/t$; this deficiency is remedied by adding triple excitations to yield SUVCCSDT. Remarkably, SGHF is already superior to SUVCCSD for most $U/t$, and SGVCCSD is better even than SUVCCSDT for the strongly correlated limit. Unlike when we add coupled-cluster correlations to SUHF, triple excitations appear to add no meaningful corrections to SGVCCSD for large $U/t$ though they are naturally important in the less strongly correlated limit. This suggests that the dominant effect of going from SUVCCSD to SUVCCSDT is that the triple excitations we add mimic the effect of the spin-projected GHF (which, of course, contains a triply-excited component). Note that lattice momentum symmetry eliminates single excitations in this problem, so SUVCCSD contains only even excitation levels. We should also point out that in the $U/t\to\infty$ limit, UHF becomes energetically exact; consequently so too do SUHF, SGHF, and their combinations with variational coupled cluster theory.

While adding SUHF or particularly SGHF to variational coupled cluster provides significant improvements over variational coupled cluster theory alone in the half-filled lattice, this is not always true. The half-filled Hubbard lattice is something of a best case scenario for methods like SGVCC, as variational coupled cluster breaks down here while SUHF or SGHF are exceptionally accurate in the strongly correlated limit. On the other hand, consider the 8-site lattice with six electrons, as shown in Fig. 2. Here, SUHF and SGHF, while significantly improving upon RHF, are no longer energetically exact in the large $U/t$ limit (i.e. they do not capture all strong correlations). In contrast, variational coupled cluster theory is already quite good, and is in fact better than SUHF or SGHF for all $U/t$. The result is that while adding SGHF to VCCSD provides a modest amount of improvement, essentially nothing is gained by adding SUHF, and VCCSDT is so close to the exact result that even adding SGHF to it has minimal effect. We should also point out that generally, our errors with respect to the exact result are significantly larger in this doped system.

V Symmetry Projection of More Elaborate Wave Functions

Traditionally, efficient formulations of symmetry projection have started by integrating the broken-symmetry wave function over a manifold of coherent states.Ring and Schuck (1980); Blaizot and Ripka (1985); Schmid (2004); Jiménez-Hoyos et al. (2012) While conceptually and computationally simple, this approach does not readily lend itself to combination with symmetry-preserving particle-hole excitation operators simply because the traditional formalism is based on combining broken-symmetry states. For example, consider adding coupled cluster to spin-projected UHF. If the cluster operator $T_{2}$ were expressed in terms of broken-symmetry orbitals, it would itself break symmetry, so one would need to spin project the unrestricted coupled cluster wave function.Duguet (2014); Duguet and Signoracci (2016); Qiu et al. If, on the other hand, the cluster operator were expressed as symmetry-preserving excitations, it would not have a particle-hole form in the basis of the broken-symmetry orbitals.

It is therefore important to ascertain how one could write the symmetry-projected wave function in terms of the symmetry-adapted orbitals. Unfortunately, for each symmetry (and possibly for each symmetry eigenvalue!) the particle-hole representation of the symmetry-projected state is different. This certainly presents a practical challenge in actually working out the appropriate projected state. While to date we have found that symmetry-projected wave functions can be written as polynomials of a small number of elementary excitation operators, we do not guarantee that this is always the case; if it is not, symmetry projection is further complicated by the need to identify increasingly cumbersome operators to account for increasingly higher excitation levels.

Nevertheless, the general program seems clear: one should find a Thouless transformation which transforms the symmetry-adapted reference determinant into the broken-symmetry determinant, and by working in the symmetry-adapted basis, one can with sufficient exertion carry out the symmetry projection analytically. We have derived but not implemented the first few terms in the number- and spin-projected generalized HFB wave function. Simpler is the form for a combination of GHF spin symmetry breaking and singlet-paired number symmetry breaking. That is, if the Thouless transformation creating the broken-symmetry mean-field is

[TABLE]

where

[TABLE]

then the number- and spin-projected wave function is just

[TABLE]

That is, because $\Sigma_{+}$ and $\Sigma_{-}$ preserve spin symmetry while $U_{0}$, $V_{+}$, and $W_{-}$ preserve number symmety, the wave function is just the product of the SGHF spin projection polynomial and the number-projected HFB polynomial. In obtaining this result we have made use of the fact that the various Thouless transformation operators of Eqn. 46 are all mutually commuting.

Perhaps more interesting yet, one can write a more general broken-symmetry coupled-cluster-style wave function in terms of symmetry breaking cluster operators acting on a symmetry-adapted determinant, and then project the result. For example, consider a number-projected case for simplicity, in which we write

[TABLE]

where $T_{\pm k}$ adds or removes $2k$ electrons from the system. The wave function would start as

[TABLE]

and would clearly be more complete than the simple number-projected HFB in which $T_{\pm 2}$ are neglected.

VI Conclusions

Incorporating symmetry projection and methods to describe dynamic correlation has been a long-standing goal of many researchers. This combination is greatly simplified by writing the symmetry-projected mean-field determinant in terms of particle-hole excitations out of some symmetry-adapted reference determinant. The particle-hole representation of number projected BCS is textbook material (see, e.g., Eqn. 8.196 of Ref. Ring and Schuck, 1980). The particle-hole form of singlet spin-projected UHF has also been recognized for quite some time,Piecuch et al. (1996) though only recently has it been realized that the SUHF wave function can be written as a polynomial of single and double excitation operators.Qiu et al. (2016, 2017) This work shows that the singlet SGHF also has such a polynomial form, though now we need not only single excitation operators ($T_{1}$) and double excitation operators ($C_{2}$) to define the polynomial, we must also include triple $(C_{3})$ and quadruple excitation operators ($K_{4}$), but apparently no further. These higher excitations appear only once the Thouless transformation from RHF to symmetry-broken Hartree-Fock has more than one mode (i.e. we require $U_{0}$, $V_{+}$, and $W_{-}$).

Having a representation of a projected Hartree-Fock wave function in terms of symmetry-adapted wave functions opens the way to several possible combinations of projected Hartree-Fock and coupled cluster theory. Here and in Ref. Qiu et al., 2017 we have explored the simple wave function ansatz $|\Psi\rangle=\exp(T)|\mathrm{PHF}\rangle$ in a variational framework as a precursor to more computationally useful approaches which use the correlator $\exp(T)$ to define a similarity-transformed Hamiltonian as in traditional coupled cluster theory. Our work on attenuated coupled clusterGomez et al. (2017) follows similar ideas. Much earlier work explored writing the PHF wave function as an exponential of higher-order cluster operators (e.g. $|\mathrm{SUHF}\rangle=\exp(T_{2}+T_{4}+T_{6}+\ldots)|\mathrm{RHF}\rangle$) and taking advantage of this structure, perhaps in combination with the factorized form of the operators $T_{n}$ obtained from PHF wave functions, to obtain accurate results for strongly correlated systems at the cost of a well-defined wave function.Piecuch et al. (1990); Piecuch and Paldus (1991); Paldus and Piecuch (1992); Piecuch et al. (1996) Most of this work to date has focused on merging SUHF and coupled cluster in some way, but it could be generalized to the SGHF case.

Clearly there is quite a lot of work to be done before we can optimally combine projected mean-field states and coupled cluster theory in practical calculations. Early returns, however, seem promising, and we have high hopes that in the fullness of time these methods will permit accurate and essentially black-box calculations for both strongly and weakly correlated systems.

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. is a Welch Foundation Chair (C-0036). TMH would like to thank Yiheng Qiu for helpful discussions and John Gomez and Carlos Jiménez-Hoyos for providing benchmark SGHF data for us to compare to.