Combining symmetry collective states with coupled cluster theory: Lessons from the Agassi model Hamiltonian

TL;DR

This paper develops a combined symmetry-projected coupled cluster approach to accurately model strongly correlated systems, demonstrated on the Agassi Hamiltonian with two broken symmetries, improving correlation energy calculations.

Contribution

It introduces a novel merger of symmetry projection with coupled cluster theory, extending previous simple model work to the complex Agassi Hamiltonian.

Findings

Effective in capturing correlation energies across the phase diagram.

Successfully models order parameters in strongly correlated regimes.

Addresses symmetry breaking in coupled cluster methods.

Abstract

The failures of single-reference coupled cluster for strongly correlated many-body systems is flagged at the mean-field level by the spontaneous breaking of one or more physical symmetries of the Hamiltonian. Restoring the symmetry of the mean-field determinant by projection reveals that coupled cluster fails because it factorizes high-order excitation amplitudes incorrectly. However, symmetry-projected mean-field wave functions do not account sufficiently for dynamic (or weak) correlation. Here we pursue a merger of symmetry projection and coupled cluster theory, following previous work along these lines that utilized the simple Lipkin model system as a testbed [J. Chem. Phys. 146, 054110 (2017)]. We generalize the concept of a symmetry-projected mean-field wave function to the concept of a symmetry projected state, in which the factorization of high-order excitation amplitudes in terms of low-order ones is guided by symmetry projection and is not exponential, and combine them with coupled cluster theory in order to model the ground state of the Agassi Hamiltonian. This model has two separate channels of correlation and two separate physical symmetries which are broken under strong correlation. We show how the combination of symmetry collective states and coupled cluster is effective in obtaining correlation energies and order parameters of the Agassi model throughout its phase diagram.

Figures

Tables

Table 1. Table 1 : Commutators of the 10 generators of the O(5) algebra of the Agassi model, in the format [row, column]. The mappings of all operators to fermions are given in Eqs. ( 20 ), ( 22 ), ( 47 ), and ( 48 ), and j 𝑗 j is 1/2 times the number of particles in the system.

	$J_{+}$	$J_{-}$	$J_0$	$A_{+1}^{†}$	$A_{+1}$	$A_{-1}^{†}$	$A_{-1}$	$A_0^{†}$	$A_0$	$N$
$J_{+}$	0	$2J_0$	$-J_{+}$	0	$-A_0$	$A_0^{†}$	0	$2A_{+1}^{†}$	$-2A_{-1}$	0
$J_{-}$	$-2J_0$	0	$J_{-}$	$A_0^{†}$	0	0	$-A_0$	$2A_{-1}^{†}$	$-2A_{+1}$	0
$J_0$	$J_{+}$	$-J_{-}$	0	$A_{+1}^{†}$	$-A_{+1}$	$-A_{-1}^{†}$	$A_{-1}$	0	0	0
$A_{+1}^{†}$	0	$-A_0^{†}$	$-A_{+1}^{†}$	0	$\left(N_{+1}-j\right)$	0	0	0	$J_{+}$	$-2A_{+1}^{†}$
$A_{+1}$	$A_0$	0	$A_{+1}$	$\left(j-N_{+1}\right)$	0	0	0	$-J_{-}$	0	$2A_{+1}$
$A_{-1}^{†}$	$-A_0^{†}$	0	$A_{-1}^{†}$	0	0	0	$\left(N_{-1}-j\right)$	0	$J_{-}$	$-2A_{-1}^{†}$
$A_{-1}$	0	$A_0$	$-A_1$	0	0	$\left(j-N_{-1}\right)$	0	$-J_{+}$	0	$2A_{-1}$
$A_0^{†}$	$-2A_{+1}^{†}$	$-2A_{-1}^{†}$	0	0	$J_{-}$	0	$J_{+}$	0	$\left(N-2j\right)$	$-2A_0^{†}$
$A_0$	$2A_{-1}$	$2A_{+1}$	0	$-J_{+}$	0	$-J_{-}$	0	$\left(2j-N\right)$	0	$2A_0$
$N$	0	0	0	$2A_{+1}^{†}$	$-2A_{+1}$	$2A_{-1}^{†}$	$-2A_{-1}$	$2A_0^{†}$	$-2A_0$	0

Equations

$$|\mathrm{RCC}⟩$$

$$T$$

$$|\Phi ⟩$$

$$T_1$$

$$Q_1$$

$$|\Phi ⟩$$

$$|\mathrm{UCC}⟩$$

$$U$$

$$|\mathrm{PUCC}⟩$$

$$|\mathrm{PRCC}⟩$$

$$|\mathrm{PHF}⟩$$

$$|\mathrm{PQCC}⟩$$

$$Q$$

$$|\mathrm{SC}⟩$$

$$|\mathrm{PQCC}⟩$$

$$|O_{\mathrm{A}}{O}_{\mathrm{B}}{O}_{\mathrm{C}}\dots ⟩$$

$$|\text{PHF}⟩$$

$$|\text{RCC}⟩$$

$$|\text{SC}⟩$$

$$|\text{PRCC}⟩$$

$$|\text{PQCC}⟩$$

$$H$$

$$J_0$$

$$J_{+}$$

$$J_{-}$$

$$A_{\sigma }^{†}$$

$$A_{\sigma }$$

$$\Pi$$

$$N$$

$$Q_1$$

$$|0⟩$$

$$\chi$$

$${\Sigma }_0$$

$$T$$

$$|\text{PHF}⟩$$

Full text

I Introduction

In single-reference coupled cluster (CC) theory, high-order particle-hole excitation amplitudes out of a reference determinant are factorized into products of lower-order particle-hole excitations via the exponential ansatz Coester (1958); Čížek (1966); Kummel et al. (1978); Shavitt and Bartlett (2009). Single-reference CC methods are at the center of modern quantum chemistry calculations Bartlett and Musiał (2007) because of their optimal combination of computational affordability and quantitative accuracy. Similarly, there has been a renaissance of the method in nuclear physics where high precision studies of medium-mass nuclei close to magic numbers were successfully performed in the last decade Dean and Hjorth-Jensen (2004); Gour et al. (2006); Hagen et al. (2014). However, single-reference CC is only accurate when applied to systems characterized by weak correlations, a category that excludes systems exhibiting such important properties as superconductivity or superfluidity Henderson et al. (2014); Duguet and Signoracci (2017), or nuclear deformation Duguet (2015), and such ubiquitous phenomena as bond breaking in quantum chemistry Bulik et al. (2015). Under weak correlation, the underlying reference wave function that coupled cluster usually takes as a given, and which in single-reference CC must be a single determinant, is a qualitatively good approximation of the ground state eigenfunction. Only if this condition is satisfied can the coupled cluster expansion of the wave function be truncated at a low order Bulik et al. (2015).

A single determinant labors under significant constraints in modeling real many-body wave functions. In particular, when the correlations are strong, it is impossible for a single-determinantal wave function to simultaneously exhibit accurate total energy and to respect the physical symmetries of the Hamiltonian - that is, in order to get a qualitatively good total energy, a mean-field treatment of the problem must break the real physical symmetries that the true many-body wave function would exhibit; this is called the “symmetry dilemma” Löwdin (1963). Coupled cluster methods can be built off of symmetry-broken determinants, such as unrestricted, quasiparticle, or Bogoliubov coupled cluster methods Henderson et al. (2014); Signoracci et al. (2015). These usually produce good total energies, but for finite systems, this comes at the cost of other properties of the wave function, including quantum numbers defined by physical symmetries such as angular momentum and particle number Sheikh and Ring (2000); Cui et al. (2013) which are lost in the broken-symmetry treatment. In the thermodynamic limit, this is not an issue, but to obtain good properties for finite systems, it is usually necessary to obtain explicitly symmetry-adapted wave functions.

It is possible to restore good quantum numbers to a broken-symmetry mean-field wave function by projecting out its symmetry-adapted component Ring and Schuck (2004); Blaizot and Ripka (1986); Schmid (2004); Scuseria et al. (2011); Bertsch and Robledo (2012); Jiménez-Hoyos et al. (2012). The resulting wave function is a linear combination of non-orthogonal determinants with good symmetries and qualitatively good energies under strong correlation. However, despite this qualitative superiority to single-reference CC, projected mean-field wave functions are still not exact and a good deal of “leftover” dynamic (i.e. weak) correlation remains to be accounted for. Furthermore, the static/strong correlation included in the symmetry projection of a mean-field reduces to the broken-symmetry mean-field in the thermodynamic limit Jiménez-Hoyos et al. (2012). A formalism for restoring particle number in quasiparticle coupled cluster has been presented Duguet and Signoracci (2017) but, to the best of our knowledge, not yet implemented.

We have recently explored the possibility of interpolating between projected mean-field and CCSD for the attractive pairing (“reduced BCS”) model exhibiting mean-field number symmetry breaking Degroote et al. (2016); Dukelsky et al. (2016), as well as explicitly combining the two approximations for molecules exhibiting mean-field spin symmetry breaking Qiu et al. (2017). In a previous work Wahlen-Strothman et al. (2017), we used the Lipkin model Lipkin et al. (1965); Meshkov et al. (1965); Glick et al. (1965) as a testbed to develop projected restricted coupled cluster (PRCC), projected quasirestricted coupled cluster (PQCC), and projected unrestricted coupled cluster (PUCC). All three methods combine mean-field symmetry breaking, symmetry-projection operators, and coupled cluster theory in various manners. They are described in detail in Sec. II.1 below.

Here we expand on this development using the Agassi model Hamiltonian Agassi (1968); Davis and Heiss (1986) as a second testbed which has more elaborate physics. The Agassi model is a schematic version of the paring-plus-quadrupole model Bes and Sorensen (1969) extensively used in nuclear physics to describe deformed superconducting nuclei. It combines the Lipkin model with the two-level pairing model Högaasen-Feldman (1961). Although its Hamiltonian is not integrable, the number of states in its collective Hilbert space grows slowly with the size of the system, as $O(n^{2})$, and therefore the Hamiltonian can be diagonalized for systems of up to hundreds of particles. However, despite its simplicity, the Agassi model is rich in non-trivial phenomena; it has two separate physical symmetries that are broken at the mean-field level in various regions of its phase diagram Davis and Heiss (1986). We test two of our new methods (PRCC and PQCC) on this model and find that they produce good agreement with exact results over the entire phase diagram. The successful tests of these methods on the Lipkin and now Agassi model Hamiltonian lay the foundation for their future application to realistic finite many-body systems. An Appendix details some algebraic properties and basis functions that are useful in performing calculations on the Agassi model.

II Theory

II.1 Symmetry-adapted formalisms

The restricted coupled cluster (RCC) wave function is

[TABLE]

in which $|0\rangle$ is a symmetry-adapted determinant, usually obtained from a mean-field approximation such as the restricted Hartree–Fock (RHF) method. The cluster operator, $T$, consists of a sum of particle-hole excitations out the reference determinant, weighted by excitation amplitudes (“$T$ amplitudes”), that do not change any of the quantum numbers associated with symmetries of the Hamiltonian. It can be expanded, and the expansion can be truncated, in terms of the number of particle-hole excitations out of the reference determinant,

[TABLE]

Because of the exponential ansatz [Eq. (1)], high-order particle hole excitations arise both from high-order terms in the expansion of the cluster operator and from products of lower-order terms. For instance, triple excitations originate both from $T_{3}$ and from $T_{1}\times T_{2}$, in terms of Eq. (2).

In the standard RCC algorithm, the Hamiltonian is subject to a similarity transformation which renders it non-Hermitian, $\bar{H}=e^{-T}He^{T}$, and the $T$ amplitudes are obtained by stipulating that the symmetry-adapted reference is a right-hand eigenstate, $\bar{H}|0\rangle=E|0\rangle$. The similarity-transformed Schrödinger equation is then left-projected against excited states, leading to a set of nonlinear equations defining residuals, $\{R_{1},R_{2},R_{3},\ldots\}$, which must be made to vanish by varying $T$ and which are coupled to each other Fan and Piecuch (2006). That is, $R_{n}$ depends on at least $T_{n}$, $T_{n+1}$, and $T_{n+2}$, so that the equations for the residuals of each order must be solved simultaneously. In practical applications this coupling chain is broken by setting amplitudes of an immediately higher order than a certain cutoff to zero; for example, in RCC with double excitations (RCCD), $T_{3}$ and $T_{4}$ are assumed to vanish. In the regime of dynamic (weak) correlation, this is a reasonable assumption and such methods give accurate energies. However, under the regime of static (strong) correlation, $T_{3}$ and $T_{4}$ become non-negligible and RCCD breaks down; the nonlinear equations for $T$ frequently do not have solutions or lead to complex energies. Alternatively, the cluster operator $T$ can be calculated using a variational method; this approach always yields real energies, but under strong correlation variational RCC (vRCC) undercorrelates substantially Qiu et al. (2017). Either way, RCC is not an appropriate ansatz for systems characterized by strong correlation.

Strong correlation also corresponds to spontaneous symmetry breaking at the mean-field level. In other words, a symmetry-adapted state such as $|0\rangle$ ceases to be the lowest-energy solution to the mean-field equations. A lower-energy broken-symmetry determinant appears, which we label $|\Phi\rangle$ and which is related to the symmetry-adapted state via a Thouless transformation,

[TABLE]

omitting a normalization factor. Here, $Q_{1}$ consists of weighted (by “$Q$ amplitudes”) single particle-hole excitations out of the symmetry-adapted reference which, unlike $T_{1}$, change the quantum numbers associated with unitary symmetries of the Hamiltonian. Note that particle-hole excitations in the same basis, such as those in $T_{1}$ and $Q_{1}$, commute with one another. An example of a symmetry-broken mean-field method is the spin-unrestricted Hartree–Fock (UHF) method used in quantum chemistry. The molecular Hamiltonian used in quantum chemistry commutes with the magnitude of the total electron spin angular momentum, $S^{2}$ and its $z$-axis component, $S_{z}$, and the two corresponding quantum numbers are labeled $S$ and $M_{S}$, respectively. The UHF method allows for mean-field symmetry breaking of the spin magnitude, $S$, but not $z$-axis component, $M_{S}$. The corresponding $T_{1}$ and $Q_{1}$ operators are

[TABLE]

where $t_{i}^{a}$ and $q_{i}^{a}$ are the $T$ and $Q$ amplitudes and $E_{a}^{i}$ and $S_{a}^{i}$ are particle-hole excitations from the $i$th occupied spatial orbital to the $a$th unoccupied spatial orbital. $E_{a}^{i}$ is an excitation which leaves both $S$ and $M_{S}$ unchanged, while $S_{a}^{i}$ is an excitation that leaves $M_{S}$ unchanged but increases $S$ by 1 Qiu et al. (2016, 2017). Combining these two excitations in the exponential ansatz [Eq. (3)] results in a UHF state, $|\Phi\rangle$, which (in general) contains multiple components with different $S$, so that $|\Phi\rangle$ itself cannot be described by any single spin quantum number.

In the Agassi model we work with below, there are no single particle-hole excitations that preserve the symmetries of the Hamiltonian, so $T_{1}$ is nonexistent. Therefore, for the rest of this work, we drop $T_{1}$. The definition of the broken-symmetry determinant is restated as

[TABLE]

Since symmetry-broken determinants appear when RCC breaks down, a simple solution is to use the symmetry-broken determinant itself as the reference for a coupled cluster expansion, resulting in the unrestricted coupled cluster (UCC) wave function,

[TABLE]

where $U$ consists of particle-hole excitations out of $|\Phi\rangle$ and, similar to $T$, is expanded and truncated in terms of particle-hole excitation order,

[TABLE]

Neither $|\Phi\rangle$ nor the UCC wave function has good quantum numbers associated with symmetries, and the individual terms in $U$ cannot easily be matched to specific transitions between symmetry quantum numbers. Although total energies calculated with the UCC wave function are usually very accurate, the results for wave function properties other than the energy are less reliable for finite systems because of the loss of symmetry.

In a previous work Wahlen-Strothman et al. (2017), we experimented with a method that applies a symmetry projection operator to a UCC wave function, which we named projected-unrestricted CC (PUCC). The PUCC wave function is

[TABLE]

where $P$ selects those components of the wave function with the desired symmetry quantum numbers. PUCC was applied to the Lipkin model system Wahlen-Strothman et al. (2017) and was found to improve on the accuracy of the total energies calculated with UCC, especially in the limit of very strong correlation. The Agassi Hamiltonian we will investigate below breaks number symmetry at the mean-field level under strong correlation, and the corresponding “unrestricted” CC method is quasiparticle or Bogoliubov CC Henderson et al. (2014); Signoracci et al. (2015). Here, we do not further discuss UCC or PUCC, which will be presented in due time, and focus on the symmetry adapted variants discussed in detail below.

An alternative to PUCC explored for the Lipkin model Wahlen-Strothman et al. (2017) is to substitute the symmetry-adapted RCC cluster operator, $T$, for $U$ in Eq. (9). The whole expression can then be written in the symmetry-adapted basis by using Eq. (6) for the broken-symmetry determinant, $|\Phi\rangle$. This is the projected-restricted CC (PRCC) wave function,

[TABLE]

Note that the projection operator commutes with a symmetry-adapted excitation operator such as $T$, so that we can move it to the right,

[TABLE]

thus expressing PRCC as a coupled cluster ansatz with a symmetry-adapted cluster operator and a projected Hartree–Fock (PHF) wave function as the reference “determinant.” However, a PHF state is actually a linear combination of several non-orthogonal determinants Jiménez-Hoyos et al. (2012) which arise from applying the projection operator to the broken-symmetry mean-field determinant. In the language of particle-hole excitations, PHF wave functions contain excitations to all orders. To see how this comes about, we can closely examine the PHF wave function

[TABLE]

where $K_{2}$ consists of symmetry-adapted double excitations which are products of single excitations in $Q_{1}$, and $F$ is some polynomial of $K_{2}$ to all orders. As discussed above, the powers of $Q_{1}$ arising from the exponential creates components of the wave function that mixes the symmetry quantum numbers albeit in a “controlled” manner. In selecting particular terms corresponding to the target quantum numbers from the exponential, the projection operator creates a new polynomial, $F$, in which high-order particle-hole excitations are expressed in terms of products of lower-order ones by some non-exponential formula. We have found in previous works that projecting out spin symmetry in this way leads to a hyperbolic sine function Qiu et al. (2016), and that projecting out number symmetry for the reduced BCS Hamiltonian results in a modified Bessel function of the first kind Degroote et al. (2016). Attempts to model these functions with an exponential ansatz will fail and this is at the root of why single-reference restricted coupled cluster is unstable under strong correlation.

Another alternative to PUCC involves a straightforward generalization of Eq. (10),

[TABLE]

where $Q$ now contains $Q_{1}$, as in PRCC, but also can contain higher-order symmetry-quantum-number-switching particle-hole excitations,

[TABLE]

This was named the projected quasirestricted coupled cluster (PQCC) method in Ref. Wahlen-Strothman et al., 2017. PRCC can be characterized as a form of PQCC in which $Q$ is restricted to $Q_{1}$. As in PRCC, each excitation in $Q$ can unambiguously be assigned a specific transition between symmetry quantum numbers, because the particle-hole basis is symmetry-adapted. However, unlike PRCC, PQCC cannot generally be interpreted as coupled cluster with a PHF reference state, because $Q$ is not limited to one-body excitations and $e^{Q}|0\rangle$ is therefore not a Thouless transformation of the symmetry-adapted reference determinant. We use the name “symmetry collective state,” or “SC,” to refer to this sort of generalization of PHF,

[TABLE]

A PHF wave function is the simplest example of a symmetry collective state. All symmetry collective states are symmetry-adapted and contain particle-hole excitations to all orders out of a symmetry-adapted reference determinant. PQCC can be described as the combination of a coupled cluster ansatz with a symmetry collective state reference,

[TABLE]

In the test on the Lipkin model Wahlen-Strothman et al. (2017), PUCC was found to give superior accuracy in the strong-correlation limit, whereas PQCC was found to give better accuracy in the recoupling region, i.e., at correlation strengths near where the lowest-energy mean-field wave function transitions between symmetry-adapted and symmetry-broken. Also, PQCC was found to have slightly superior accuracy in the energies compared to PRCC at higher orders of truncation (i.e., triples or higher) because of the additional excitations included in $Q$.

Throughout the remainder of this work, we introduce a symbolic nomenclature in which terms all of the wave functions hitherto described (other than UCC and PUCC) can be expressed:

[TABLE]

where $O_{\mathrm{A}}$, $O_{\mathrm{B}}$, and $O_{\mathrm{C}}$ are particle-hole excitation operators in the basis of the symmetry-adapted reference state, $|0\rangle$. In this nomenclature, PHF, RCC, a general symmetry collective state, PRCC, and PQCC are

[TABLE]

Note that in the case of methods such as RCC, which include no symmetry-breaking operators, the projection operator has no effect.

II.2 Agassi Model

In this work, we apply the PRCC and PQCC methods to a more sophisticated testbed than the Lipkin model on which they were originally developed: the Agassi model. According to the labeling conventions of Davis and Heiss Davis and Heiss (1986), the Agassi model describes a set of $2j$ spinless fermions occupying two sets of single-particle states, labeled $\sigma=-1$ and $\sigma=+1$, each of which is $2j$-fold degenerate. The individual states in each level are furthermore labeled with an index $m$, with $-j\leq m\leq j$ and $m\neq 0$. The Hamiltonian operator for this model is

[TABLE]

where $\epsilon$, $g$, and $V$ are adjustable parameters and where the one-body operators, $J$ and $A$, are

[TABLE]

That the size of the Hilbert space for this model grows quadratically with the number of particles, as mentioned in Sec. I, is demonstrated in the Appendix.

However, despite this simplicity, the physics of the Agassi model is rich in non-trivial phenomena. The Agassi model is, in essence, a union of the Lipkin model Lipkin et al. (1965); Meshkov et al. (1965); Glick et al. (1965) and the two-level pairing model Högaasen-Feldman (1961). The Hamiltonian, Eq. (19), combines two qualitatively distinct channels of many-body correlation. The second term in Eq. (19), involving $J_{-}$ and $J_{+}$ and with strength parameter $V$, is equal to the correlation term of the Lipkin model Hamiltonian. It moves pairs of particles in either the lower or upper level “vertically” up or down, between $\sigma=-1$ and $\sigma=+1$, but cannot change the “horizontal” indices $m$. On the other hand, the third term in Eq. (19), involving $A_{\sigma}$ and $A_{\sigma}^{\dagger}$ and with strength parameter $g$, is equivalent to the two-body term in the pairing Hamiltonian. It moves pairs of particles with opposite $m$ indices (e.g. $\sigma,m$ and $\sigma,-m$) to another pair of states at $\sigma^{\prime},m^{\prime}$ and $\sigma^{\prime},-m^{\prime}$. A diagrammatic representation of all this is depicted in Fig. 1.

The two distinct channels of correlation correspond to the breaking of two separate symmetries of the Agassi Hamiltonian at the mean-field level when the respective correlation strength parameters are large enough. The two relevant symmetries here are parity, which is inherited from the Lipkin model and which breaks at the mean-field level when $V$ is large, and particle number, which breaks at the mean-field level when $g$ is large, as in the two-level pairing model Högaasen-Feldman (1961). The parity symmetry operator, $\Pi$, and particle number operator, $N$, are

[TABLE]

and the Hamiltonian has $[H,\Pi]=[H,N]=0$. A state has $\Pi$-eigenvalue $+1$ or $-1$ depending on whether it has an even or odd number of particle-hole excitations from the lower level to the upper level and whether $j$ is even or odd.

The general mean-field wave function of the Agassi model which can break these two symmetries is given by Eq. (6) with

[TABLE]

where the symmetry-adapted reference determinant is

[TABLE]

and where $|-\rangle$ is the physical vacuum. Eq. (24) describes the full occupation of the lower level of states ($\sigma=-1$) created by the maximum possible number of $A_{-1}^{\dagger}$ pair creation operators, and the full vacancy of the upper level ($\sigma=+1$), indicated by the absence of any $A_{+1}^{\dagger}$ pair creation operators or any $J_{+}$ particle-hole excitations. In Eq. (23), the two terms correspond to excitations along the two Agassi correlation channels and to breaking of the two Agassi symmetries. Nonzero $q_{\mathrm{L}}$ corresponds to a parity-broken mean-field and nonzero $q_{\mathrm{P}}$ corresponds to a number-broken mean-field. Throughout the rest of the paper, we refer to the symmetry-adapted reference state as the “RHF” state and the general mean-field determinant $|\Phi\rangle$ as the “UHF” state which is equivalent of the Hartree-Fock-Bogoliubov vacuum of Ref. Davis and Heiss, 1986. Figure 2 depicts the mean-field phase diagram; here and throughout the following we have made the change to dimensionless variables of Davis and Heiss Davis and Heiss (1986) from $\epsilon$, $V$, and $g$ to

[TABLE]

which are defined so that mean-field parity symmetry breaking occurs for $\chi>1$ and mean-field number symmetry breaking occurs for $\Sigma_{0}>1$.

The breaking of these symmetries at the mean-field level signals the impending failure of typical single-reference methods to model the structure of the Agassi model (or of any model Hamiltonian, for that matter). Figure 3 tracks this failure for the RCCD method applied to an Agassi model with $j=20$. The RCCD wave function [$|T_{2}\rangle$ in the nomenclature of Eq. (17)] is given by Eq. (1) with

[TABLE]

In Fig. 3, we obtain the amplitudes $t_{\mathrm{LL}}$ and $t_{\mathrm{PP}}$ through the standard CC algorithm based on a similarity transformation; the existence of real solutions to the corresponding equations is not guaranteed. Under weak correlation; i.e., $\chi<1$ and $\Sigma_{0}<1$, the RCCD method is nearly exact. However, near $\chi=1$ or $\Sigma_{0}=1$, RCCD begins to overcorrelate, and starting shortly thereafter, we were unable to find solutions to the nonlinear equations. Under strong correlation, $T_{3}$ and $T_{4}$ becomes non-negligible; RCCD calculations assume that it is zero and thus fail.

The failure of restricted CC is due to its inability to model the strong correlation captured by the symmetry breaking at $\chi=1$ and $\Sigma_{0}=1$ in the mean-field wave function. Starting at $x=1$ in Fig. 3, the UHF solution to the Agassi model begins to diverge from the RHF state. The energy errors of the UHF solution are certainly much lower than the errors of the RHF energy. Going from UHF to PHF [$|Q_{1}\rangle$ in terms of Eq. (17)] further improves the energies by a significant margin. (A variation-after-projection method is used to obtain the PHF results in Fig. 3 and everywhere else in this work.) For the Agassi model the PHF wave function is explicitly

[TABLE]

where $Q_{1}$ is given by Eq. (23) and where $I_{0}$ is a modified Bessel function of the first kind. The projection operator selects those terms with even powers of $J_{+}$, which preserve parity symmetry, and those terms with equal powers of $A_{+1}^{\dagger}$ and $A_{-1}$, which preserve particle number. The resulting excitations are separable into a product of a hyperbolic cosine of $J_{+}^{2}$ (as reported in the work on Lipkin parity projection in Ref. Wahlen-Strothman et al., 2017) and a Bessel function of $A_{+1}^{\dagger}A_{-1}$ (as reported in the work on number projection in Ref. Degroote et al., 2016).

The PHF wave function is exact for the Agassi model in the strong attractive pairing limit ($\Sigma_{0}\to\infty$ and $\Sigma_{0}\gg\chi$) and in the simultaneous thermodynamic and strong Lipkin correlation limit ($\chi\to\infty$, $\chi\gg\Sigma_{0}$, and $j\to\infty$). The behaviors of the red ($\Sigma_{0}>\chi$) and green ($\chi>\Sigma_{0}$) curves in Fig. 3 are suggestive of these limits. There is still visible error in the PHF energies, however, especially near the recoupling region. We attribute this to weak (dynamical) correlation and therefore we expect that this error can be substantially cured with our PRCC approach. Augmenting Eq. (28) with the doubles cluster operator of RCCD results in a PRCC wave function which, in the nomenclature of Eq. (17), is labeled $|T_{2}Q_{1}\rangle$:

[TABLE]

Note that, although the $|T_{2}Q_{1}\rangle$ wave function has more parameters than the PHF ($|Q_{1}\rangle$) wave function, the excitations are still separable into a product of two polynomials, one for the $J_{+}^{2}$ excitations related to the Lipkin correlation channel and one for the $A_{+1}^{\dagger}A_{-1}$ excitations related to the pairing correlation channel.

There is one limit of strong correlation in which PHF is not exact regardless of system size: the simultaneous strong correlation limit along both channels, $\Sigma_{0}=\chi\to\infty$. The PHF energy error along that line in Fig. 3 (blue) grows monotonically with correlation strength instead of converging, suggesting a non-dynamical correlation effect not captured by PHF in this region of the phase diagram. We conjecture, and confirm in Sec. III below, that the relatively poorer performance of PHF in this region is related to the lack of any terms coupling correlation along the two channels. An analogy can be made to the difference between a wave function obtained in an adiabatic approximation separating two different interaction channels and one which includes non-adiabatic coupling. We anticipate that the latter wave function will outperform the former in modeling a system in which the characteristic excitations of the two channels have similar energy scales.

The lowest-order way to couple $J_{+}^{2}$ excitations to $A_{+1}^{\dagger}A_{-1}$ excitations is to include a higher-order $Q$ term:

[TABLE]

Although we call this $Q_{2}$, we note that a number-breaking, parity-preserving term with $A_{+1}^{\dagger}A_{+1}^{\dagger}+A_{-1}A_{-1}$ as the operator part should formally be added to Eq. (30). We don’t expect this latter excitation to be important, so we omit it here. $Q_{2}$ as given by Eq. (30) changes the quantum numbers of parity and number symmetry simultaneously. Adding $Q_{2}$ to the excitations included in PHF generates a more general symmetry collective state labeled $|Q_{1}Q_{2}\rangle$:

[TABLE]

Unlike PHF ($|Q_{1}\rangle$), $|Q_{1}Q_{2}\rangle$ is not a symmetry projection of a mean-field wave function, because the exponential of $Q_{2}$ cannot be considered a Thouless transformation. Making the same modification to the $|T_{2}Q_{1}\rangle$ wave function generates a PQCC method labeled $|T_{2}Q_{1}Q_{2}\rangle$,

[TABLE]

We expect $|Q_{1}Q_{2}\rangle$ and $|T_{2}Q_{1}Q_{2}\rangle$ to outperform PHF ($|Q_{1}\rangle$) and $|T_{2}Q_{1}\rangle$, respectively, in the region of the phase diagram where both correlation strength parameters are simultaneously high.

Another potential way to couple the two correlation channels is in the cluster operator; i.e., by including terms such as

[TABLE]

which is similar to $Q_{2}$ in that it excites along both correlation channels simultaneously. However, in this work we seek a method that can handle strong correlation at the doubles level or lower, and $T_{4}$ contains quadruple excitations. Furthermore, CC methods excel at describing weak, dynamic correlation, and the coupling between pairing and Lipkin channels is expected when correlation along both channels is strong. For these reasons we do not here pursue methods using $T_{4}$.

III Results

We test a PRCC wave function ($|T_{2}Q_{1}\rangle$) and a PQCC wave function ($|T_{2}Q_{1}Q_{2}\rangle$) on the Agassi model and examine both the accuracy of the computed energies and the quality of the wave function, as measured by several order parameters. We compare $|T_{2}Q_{1}\rangle$ and $|T_{2}Q_{1}Q_{2}\rangle$ to the results from exact diagonalization of the Agassi model (“FCI”) as well as a variational implementation of RCCD ($|T_{2}\rangle$) and the PHF ($|Q_{1}\rangle$) and $|Q_{1}Q_{2}\rangle$ wave functions. We take two separate approaches to determining the $T$ and $Q$ amplitudes: first we test a variational formalism in which all amplitudes are determined by minimizing the energy expectation value over a given wave function, and then we test a similarity-transformed basis approach more common to coupled cluster theory. The former is robust and has guaranteed solutions for any point on the phase diagram (although obtaining those solutions may not always be computationally easy), whereas the latter is more practical for large systems where FCI codes are not available. We continue to use the $j=20$ system which is large enough to significantly reduce finite-size effects.

III.1 Variational Energies

Our variational treatments of CC theories are here denoted by the prefix “v”. We make a Hermitian energy expectation value stationary with respect to all wave function parameters simultaneously. For instance, the vRCCD energy is obtained by solving the set of equations,

[TABLE]

In this section we evaluate the accuracy of various wave functions in recovering the correlation energy. “Correlation energy” in this case is defined with relation to the RHF energy:

[TABLE]

RHF is used as a reference point instead of the symmetry-broken mean field because the latter happens to give exact energies in various limits of the Agassi model Henderson et al. (2014); Wahlen-Strothman et al. (2017), limiting its usefulness as a reference in measuring the accuracy of approximate methods.

The error from FCI of the vRCCD ($|T_{2}\rangle$ wave function) correlation energy across the phase diagram is plotted in Fig. 4(a). The vRCCD method fails quantitatively outside of the weak-correlation region at low $\chi$ and low $\Sigma_{0}$; at higher correlation strengths vRCCD systematically undercorrelates by 5-10%. There is no failure to converge the equations, as depicted in Fig. 3, but that is because of the inherent robustness of the variational approach, as contrasted to the projective method more commonly used in coupled cluster calculations Shavitt and Bartlett (2009) and which RCCD (without the “v”) utilizes. The vRCCD results still confirm that coupled cluster is poorly suited to modeling strong correlation.

In contrast, the correlation energy error for PHF ($|Q_{1}\rangle$ wave function) is plotted in Fig. 4(b). The PHF energy has notably poorer accuracy than the vRCCD energy in the weakly-correlated region of the phase diagram. However, at either strong correlation limit, PHF outperforms vRCCD, undercorrelating by only up to 1% in the strong Lipkin correlation region ($\chi>\Sigma_{0}$) and achieving accuracy to within 0.01% of the exact correlation energy in the strong pairing region ($\Sigma_{0}>\chi$). The better performance in the strong pairing limit is attributable to the fact that PHF reduces to the PBCS wave function in that limit, and in the strong-correlation limit of the attractive pairing model the PBCS wave function is exact Richardson (1977); Román et al. (2002); Degroote et al. (2016). The limit of the PHF wave function under strong Lipkin correlation, meanwhile, is also exact for the Lipkin model, but only in the simultaneous thermodynamic and strong-correlation limits. The worse performance of PHF in the strong Lipkin region compared to the strong pairing region is thus attributable to finite-size error.

Surprisingly, the $|Q_{1}Q_{2}\rangle$ wave function does not appear to improve significantly on the PHF results; indeed, the difference between the two in Figs. 4(b) and 4(c) is barely visible. Because these are variational methods and because $|Q_{1}Q_{2}\rangle$ has one more parameter than does $|Q_{1}\rangle$, the energy of the $|Q_{1}Q_{2}\rangle$ wave function is indeed lower than that of PHF, but not by a significant amount. This seems to falsify our hypothesis discussed in Sec. II.2 that $|Q_{1}Q_{2}\rangle$ would make a critical improvement in the region where both correlation strengths are simultaneously large. The only region of the phase diagram where any improvement over PHF is visible is the strong-Lipkin region that, as discussed above, is subject to finite-size error.

Regardless, none of these three methods alone provide good-quality correlation energies across the entire phase diagram. However, they complement one another when combined. The energy error of vPRCC ($|T_{2}Q_{1}\rangle$), which is a combination of PHF and vRCCD, is plotted in Fig. 5(a). The vPRCC energy is uniformly superior to both vRCCD and PHF, which is perhaps unsurprising since these are all variational methods and the $|T_{2}Q_{1}\rangle$ wave function has the greatest number of parameters of these three. However, the details are significant: the vPRCC energy is nearly exact everywhere except in the regions near phase transitions. The worst result is along the diagonal ray $\chi=\Sigma_{0}>1$, where vPRCC undercorrelates 1-2%. In this region of the phase diagram, both number and parity symmetry are broken at the mean-field level, neither correlation channel dominates the other, and an adiabatic-like ansatz in which the excitations are factorized into products of excitations along each channel is inappropriate.

Going from vPRCC to vPQCC (that is, going from $|T_{2}Q_{1}\rangle$ to $|T_{2}Q_{1}Q_{2}\rangle$) improves the quality of the correlation energy in this region by an order of magnitude, as also shown in Fig. 5(a). Now the hypothesis we discussed in Sec. II.2, that $Q_{2}$ is critical to include when both correlation strengths are large, is confirmed. Clearly, the basic reason for this improvement is that $Q_{2}$, which is included in the PQCC ($|T_{2}Q_{1}Q_{2}\rangle$) and $|Q_{1}Q_{2}\rangle$ wave functions but not the PRCC ($|T_{2}Q_{1}\rangle$) or RCCD ($|T_{2}\rangle$) wave functions, couples Lipkin-like excitations to pairing-like excitations at the level of the symmetry collective state, which is the part of our wave function that is designed to account for strong-correlation effects. However, it does seem that this effect only appears when combined with other double excitations, as in $|T_{2}Q_{1}Q_{2}\rangle$, and not when it is the sole double term, as in $|Q_{1}Q_{2}\rangle$.

III.2 Quality of the variational wave function

While PRCC ($|T_{2}Q_{1}\rangle$) and PQCC ($|T_{2}Q_{1}Q_{2}\rangle$) yield highly accurate energies, it remains to see whether those wave functions are accurate for properties other than energies. We seek to test the accuracy of the Agassi wave function by probing several order parameters derived from the wave function. The simplest is a measure of the average level of excitation; that is, a measure of the number of particles in the upper level ($\sigma=+1$) on average:

[TABLE]

This is a straightforward measure of the strength of correlation in the wave function and does not distinguish between Lipkin-channel and pairing-channel correlation: the closer our wave function is to RHF, the closer $n$ will be to zero.

We additionally probe order parameters which track the effects of the two correlation channels independently. For mean-field wave functions that are allowed to break symmetry, the obvious candidates involve expectation values of the excitation operators that produce the broken-symmetry mean field:

[TABLE]

Note that $\Delta$ is the gap in the BCS sense Ring and Schuck (2004). However, these parameters cannot be used as written because the approximate wave functions we are using (not to mention the exact wave function) are symmetry-adapted, so $J$ and $\Delta$ as defined in Eqs. (38) and (39) are zero by construction. Therefore, for the purpose of tracking correlation strength in these wave functions, we modify Eqs. (38) and (39) to read

[TABLE]

which, if evaluated with a UHF wave function, give the same results as Eqs. (38) and (39) in the thermodynamic limit.

Figure 6 plots the values of $n$, $J$, and $\Delta$ calculated using Eqs. (37), (40), and (41) from the exact (FCI) wave function. It can be seen that $n$ is large whenever the system is strongly correlated along either channel, whereas $J$ and $\Delta$ are large only when Lipkin correlation or pairing correlation, respectively, is strong.

Holding the above in mind, we examine the errors from FCI of the three above-described order parameters for the vPRCC ($|T_{2}Q_{1}\rangle$) and vPQCC ($|T_{2}Q_{1}Q_{2}\rangle$) methods, which are plotted for two lines through the phase diagram in Fig. 7. These order parameters are estimated to within 0.5% in most regions of the phase diagram by both vPRCC and vPQCC. In the weakly correlated region, the values of all three order parameters are nearly zero; in the strongly Lipkin-correlated regions and strongly pairing-correlated regions, $|T_{2}Q_{1}\rangle$ and $|T_{2}Q_{1}Q_{2}\rangle$ display minor errors. However, in the simultaneous strongly Lipkin-correlated and pairing-correlated region (the center of Fig. 7(b)), the errors of $|T_{2}Q_{1}\rangle$ estimation of all three parameters jump to 1-10%. Nevertheless, the $|T_{2}Q_{1}Q_{2}\rangle$ wave function substantially cures this error. These results further support the inference above that $|T_{2}Q_{1}Q_{2}\rangle$ accounts for coupling between the two correlation channels which is important when the two are of similar strength and both symmetries are broken - the region where $|T_{2}Q_{1}Q_{2}\rangle$ is a significant improvement over $|T_{2}Q_{1}\rangle$ in the order parameters is the same region where it is an improvement in the energies.

III.3 The Similarity-Transformation Formalism

Although the variational method is robust, it is not practical for realistic systems because variational coupled cluster equations do not truncate until the order of excitations reaches the number of particles Voorhis and Head-Gordon (2000). Here we experiment with replacing the variational formalism with a more standard projective formalism Shavitt and Bartlett (2009), in which the Hamiltonian is similarity-transformed with the coupled-cluster operator and then left-projected against selected excitations.

The Schrödinger equation, when approximated with a PRCC or PQCC ansatz, is

[TABLE]

where $|Q\rangle$ is the relevant symmetry collective state [Eqs. (17) and (18)], i.e. PHF ($|Q_{1}\rangle$) for $|T_{2}Q_{1}\rangle$ or $|Q_{1}Q_{2}\rangle$ for $|T_{2}Q_{1}Q_{2}\rangle$. We multiply by $e^{-T}$ on the right to obtain our similarity transformation,

[TABLE]

which renders the transformed Hamiltonian, $\bar{H}$, non-Hermitian. We project the Schrödinger equation against the left-hand eigenstate of $\bar{H}$, which we parameterize as

[TABLE]

where the prime on $\langle Q^{\prime}|$ indicates that the $Q$ amplitudes defining the symmetry collective state are, in general, allowed to differ from those in the right-hand eigenstate, $|Q\rangle$, and where $Z$ contains de-excitation operators conjugate to the excitations in $T$. For instance, $Z_{2}$ is

[TABLE]

The overall energy expression is

[TABLE]

and we make the energy stationary with respect to all $T$, $Z$, and $Q$ amplitudes.

Note that we have only similarity-transformed with the CC part of our method, and left the symmetry collective state essentially variational, or rather bivariational, as the bra and ket have different $Q$ amplitudes. That is, we have not explicitly performed a polynomial similarity transformation Degroote et al. (2016); Qiu et al. (2016), $F^{-1}(Q)HF(Q)$ where $F(Q)$ is defined by $F(Q)|0\rangle=Pe^{Q}|0\rangle$. If we had, it would be necessary to determine the form of the left-hand eigenstate of the doubly-similarity-transformed Hamiltonian, which is not trivial. This problem is discussed for the case of spin-projected unrestricted Hartree–Fock (SUHF) in Ref. Qiu et al., 2016.

Figure 8 shows the error from FCI of the PRCC ($|T_{2}Q_{1}\rangle$) and PQCC ($|T_{2}Q_{1}Q_{2}\rangle$) methods, respectively. The accuracy is slightly reduced in both cases from the variational results, but not qualitatively different. $|T_{2}Q_{1}\rangle$ is extremely accurate in all except recoupling regions and diagonal strong correlation, and $|T_{2}Q_{1}Q_{2}\rangle$ remedies its defect in the latter, with correlation energies accurate to within 1%. These results confirm that the combination of coupled cluster with symmetry collective states improves dramatically on the results demonstrated in Fig. 3 in Sec. II.2; at every point on the phase diagram, real solutions to the equations for the $T$ amplitudes are obtained, despite the fact that we have discarded the variational method.

IV Conclusion

The PRCC and PQCC methods, initially tested on the Lipkin model Hamiltonian Wahlen-Strothman et al. (2017), are shown to very effectively estimate the wave function of the Agassi model, recovering more than 99% of correlation energy and predicting order parameters of the wave function with very high accuracy. While the previous work Wahlen-Strothman et al. (2017) demonstrated the principle by which such methods could model systems under strong correlation, here we have shown how they can be used for systems with multiple distinct correlation channels, including multiple simultaneous motifs of strong correlation. Certainly, more work is necessary in generalizing this formalism to realistic nuclear and electronic structure problems of nuclear physics, quantum chemistry, and condensed matter physics. However, the fact that methods with between 5 and 8 parameters including no excitation amplitudes beyond the level of doubles are so powerful for this simple model system is very promising about the potential of symmetry-collective-state-plus-CC approaches to strong correlation.

Acknowledgements.

We thank Jacob M. Wahlen-Strothman and Dr. Matthias Degroote for providing reference data for testing purposes, and Dr. Thomas M. Henderson for assistance in overcoming computational problems. 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 (No. C-0036). Computational resources for this work were supported in part by the Big-Data Private-Cloud Research Cyberinfrastructure MRI-award funded by NSF under grant CNS-1338099 and by Rice University. J.D. acknowledges support from the Spanish Ministry of Economy and Competitiveness and FEDER through Grant No. FIS2015-63770-P.

Appendix A Algebra of the Agassi model

The Hamiltonian of the Agassi model [Eq. (19)] is given in Sec. II.2 terms of seven one-body operators, whose mapping to fermions is given in Eqs. (20). If we add to these seven the number operator [$N$, Eq. (22)] along with $A_{0}^{\dagger}$ ($A_{0}$), which creates (destroys) a particle $\sigma=+1$ and another one in $\sigma=-1$,

[TABLE]

then the combined set of ten one-body operators closes an O(5) algebra. All commutators among the ten $A$, $J$, and $N$ operators are summarized in Table 1. Note that the SU(2) algebra of the Lipkin model and the three SU(2) algebras of the pairing model are subalgebras of this O(5) Klein et al. (1982). Other useful one-body operators can be defined as linear combinations of these ten generators; for instance, the number of particles in the upper and lower levels are given by

[TABLE]

Following Ref. Klein et al. (1982), a complete (although non-orthogonal) set of states for the Agassi model can likewise be defined in terms of pair creation ($A^{\dagger}$) operators acting on the physical vacuum:

[TABLE]

which is complete because all other generators in the $O(5)$ algebra annihilate the physical vacuum to the right. Furthermore, if we seek only the states with $2j$ particles, then we have $n_{-}+n_{+}+n_{0}=j$ and the size of the Hilbert space is therefore quadratic in $j$, as asserted in Sec. II.2. As for parity symmetry, even or odd sectors correspond simply to even or odd $(n_{0}-j)$. The symmetry-adapted non-interacting ground state ($|0\rangle$ in the main text) in these terms is $|j,0,0\rangle$.

By repeatedly applying the commutation relationships summarized in Table 1, the effects of all ten generators on a state $|n_{-},n_{+},n_{0}\rangle$ can be evaluated:

[TABLE]

Equations (50e), (50g), and (50i) can be utilized to construct recursive equations that are used to evaluate the overlap matrix of the states $|n_{-},n_{+},n_{0}\rangle$. Keeping in mind that states with differing $N$ and $J_{0}$ eigenvalues are orthogonal, and that $\langle 0,0,0|0,0,0\rangle=\langle-|-\rangle=1$, we have

[TABLE]

Given the overlap matrix [Eqs. (51)] and the actions of all ten generators on the non-orthogonal basis [Eqs. (50)], an orthonormal basis and matrix elements of all ten operators are easily obtained. Since, as demonstrated, the size of the basis (at least in the symmetry-adapted space) grows quadratically with $j$, the overlap and (symmetry-adapted products of) operator matrices in the symmetry-adapted space require only $O(j^{4})$ storage space. The largest computational difficulty is numerical instability in systems of $10^{2}$ particles or more, arising from overlap matrix elements with values that differ by more than 16 orders of magnitude, which is partially resolved by using quadruple-precision storage of floating-point numbers. This means that FCI and all of the symmetry-adapted methods described in this work are easily implemented using complete-basis representations of all operators. The symmetry-broken mean-field is another matter, however this can be calculated very easily using simple closed-form equations as explained by Davis and Heiss Davis and Heiss (1986).