Exploring Ground and Excited States via Single Reference Coupled-Cluster Theory and Algebraic Geometry

Abstract

The exploration of the root structure of coupled cluster equations holds both foundational and practical significance for computational quantum chemistry. This study provides insight into the intricate root structures of these non-linear equations at both the CCD and CCSD level of theory. We utilize computational techniques from algebraic geometry, specifically the monodromy and parametric homotopy continuation methods, to calculate the full solution set. We compare the computed CC roots against various established theoretical upper bounds, shedding light on the accuracy and efficiency of these bounds. We hereby focus on the dissociation processes of four-electron systems such as (H₂)₂ in both D_2h and D_∞h configurations, H₄ symmetrically distorted on a circle, and lithium hydride. We moreover investigate the ability of single-reference coupled cluster solutions to approximate excited state energies. We find that multiple CC roots describe energies of excited states with high accuracy. Remarkably, our investigations reveal that for systems like lithium hydride, CC not only provides high-accuracy approximations to several excited state energies but also to the states themselves.

1 Introduction

Coupled cluster (CC) theory is a high-precision wave function method in computational quantum chemistry ¹. Its origins date back to 1958 when Coester proposed to use an exponential parametrization of the wave function ², i.e.

$$|{\Psi}\rangle=e^{\hat{T}}|{\Phi_{0}}\rangle,$$ (1)

where $\hat{T}$ is the new unknown, and $|{\Phi_{0}}\rangle$ is a chosen reference Slater determinant (commonly the Hartree-Fock state). This parametrization was independently derived by Hubbard ³ and Hugenholtz ⁴ in 1957 as an alternative to summing many-body perturbation theory contributions order by order. A pivotal moment in the development of CC theory was Čížek’s seminal work in 1966 ⁵. Čížek presented the first derivation of the CC working equations together with corresponding simulations. This work established CC theory within the theoretical framework that we recognize today, including second quantization applied to many-fermion systems, normal ordering, contractions, Wick’s theorem, normal-ordered Hamiltonians and Goldstone-style diagrammatic techniques. There exist several reviews of CC theory, summarizing different aspects of the method’s success over the past decades, see e.g. Kümmel ⁶, Čížek ⁷, Bartlett ⁸, Paldus ⁹, Arponen ¹⁰ and Bishop ¹¹.

Employing the exponential ansatz in Eq. (1) yields the following reformulation of the stationary Schrödinger equation:

$$\left\{\begin{aligned} \langle{\Phi_{0}}|e^{-\hat{T}}\hat{H}e^{\hat{T}}|{\Phi_{0}}\rangle&=E,\\ \langle{\Phi}|e^{-\hat{T}}\hat{H}e^{\hat{T}}|{\Phi_{0}}\rangle&=0&\forall|{\Phi}\rangle\perp|{\Phi_{0}}\rangle.\\ \end{aligned}\right.$$ (2)

Employing the standard Galerkin projection, the latter condition yields the projected CC equations (vide infra), i.e. a set of nonlinear algebraic equations. Typically, a single solution to this set of equations is computed using (quasi) Newton-type methods ¹². However, as a set of nonlinear algebraic equations, the CC equations possess multiple roots, raising the following natural questions:

• How many solutions to the CC equations exist?

• Which solution describes the ground state best?

• Do some of the solutions describe excited states?

Especially the last question is controversial, to say the least. On one side, it can be proven that as long a state is not strictly orthogonal to the chosen reference state, CC theory can – in its untruncated form – describe that state exactly. On the other side, excited states commonly have poor overlap with the Hartree-Fock ground state reference – hindering the convergence of (quasi) Newton-type methods. Besides such theoretical arguments, it has been numerically shown that higher-order CC solutions can approximate excited state energies well while the CC solutions provide poor approximations to the actual states ¹³. Conventionally, excited states are computed in a post-single-reference CC manner. This is done by either constructing and diagonalizing the similarity-transformed Hamiltonian in equation-of-motion CC ^{14, 15, 16, 17, 18, 19, 20} or by examining the poles of the linear-response function in linear-response CC ^{16, 21, 22, 23}. In either case, the excited states provided by these methods are inherently biased toward the lowest eigenstate upon which they are built. As a direct consequence, the accuracy of such approaches is known to be quite poor for certain classes of eigenstates that are energetically higher than the ground state ^{24, 25, 26, 27, 28, 29}.

Besides the fundamental nature of questions i)–iii), the existence of multiple roots to the CC equations can present challenges for the practical implementations of the method. For example, the roots may be difficult to adequately converge to ^{5, 30, 31, 32}, and the convergence properties of the iterative solution methods can strongly depend on the employed initial guess ^{33, 34}. See also Ref. 35 for a mathematical exposition of this problem and Ref. 36 for a potential remedy.

In this work, we elaborate on recent advances coming from algebraic geometry that provide new ideas on how to bound the number of roots to the CC equations. These new and significantly improved bounds can then be used to compute all roots to the CC equations, answering the above-posed questions i) – iii). In particular, we demonstrate that higher-order solutions of the CC equations can effectively approximate not only the energies of excited states but also the eigenstates themselves.

1.1 Previous works and perspective

We highlight that the inaugural investigation of the root structure of CC methods dates back to 1978 when Zivkovic and Monkhorst scrutinized the singularities and multifarious solutions in single-reference CC equations ³⁷. Building upon this work, Paldus and colleagues subsequently conducted mathematical and numerical analyses during the early 1990s, focusing on the manifold of solutions of both single-reference and state-universal multi-reference CC equations, while also elucidating their singularities and analytical characteristics ^{38, 39}. In 1998, Kowalski and Jankowski revived homotopy methods in connection with the single-reference CC theory and used them to solve a CCD system ⁴⁰. Kowalski with Piecuch then extended the application of homotopy methods to the CCSD, CCSDT, and CCSDTQ ⁴¹ for a 4-electron system described in a minimal basis set. From these investigations, they derived the $\beta$-nested equations and proved the Fundamental Theorem of the $\beta$-NE Formalism ⁴², which enabled an explanation of curves connecting multiple solutions of the various CC polynomial systems, i.e. from CCSD to CCSDT, CCSDT to CCSDTQ, etc. This Kowalski-Piecuch homotopy was recently analyzed in-depth using topological degree theory ^{43, 44}. Moreover, Kowalski with Piecuch used homotopy methods to obtain all solutions of the generalized Bloch equation ⁴⁵ and explain symmetry breaking processes in the CC formulation ^{46, 47, 48, 49, 50}. It is noteworthy that these fundamental results of CC theory have spurred the development of the method of moment CC approach ^{51, 52, 53, 54}, which in turn has led to the highly efficient completely renormalized CC approach ^{55, 56, 57}. This emphasizes the pivotal role of this research trajectory in driving methodological innovations within CC theory. A perspective on the recent use of homotopy methods in the context of quantum chemistry can be found in Ref. 58.

A more algebraic and mathematically driven approach to this problem has recently gained momentum. Initialized by some of the authors, a first algebraic perspective was introduced in Ref. 59. This was followed by an in-depth algebraic investigation leading to the introduction of the CC truncation varieties ⁶⁰ – which are the subject of this work. This work moreover contained open mathematical questions and conjectures, one of which was proven shortly after ⁶¹. More generally, these works are part of a presently increasing interest in the mathematical understanding of CC theory, driven by the applied mathematics community. In this context, Schneider performed the first local analysis of CC theory in 2009 employing a functional analytic framework ⁶². Subsequently, this approach was refined ^{63, 64} and applied to different CC methods ^{65, 66, 67}. A more versatile framework for analyzing general CC variants using topological degree theory was later proposed by Csirik and Laestadius ^{43, 44}. The most recent numerical analysis results regarding single reference CC were established by Hassan, Maday, and Wang ^{68, 69}.

2 Coupled-cluster theory

The starting point of CC theory is the exponential ansatz in Eq. (1), leading to the CC energy expression

$$E=\langle{\Phi_{0}}|e^{-\hat{T}}\hat{H}e^{\hat{T}}|{\Phi_{0}}\rangle,$$ (3)

where the cluster operator $\hat{T}$ is determined through the CC equations

$$0=\langle{\Phi}|e^{-\hat{T}}\hat{H}e^{\hat{T}}|{\Phi_{0}}\rangle,\qquad\forall|{\Phi}\rangle\perp|{\Phi_{0}}\rangle.$$ (4)

The cluster operator is a linear combination of elementary particle-hole excitation operators, i.e.

$$\hat{T}({\bf t})=\sum_{\mu>0}t_{\mu}\hat{X}_{\mu},$$ (5)

where we highlight the dependence of $\hat{T}$ on its expansion coefficients by writing $\hat{T}(\bf t)$. Recall that elementary particle-hole excitation operators $\hat{X}_{\mu}$ are merely compositions of projections onto a subset of single-particle basis elements (the occupied orbitals) and wedge products with another subset of single-particle basis elements (the virtual orbital). It is therefore convenient to label the elementary particle-hole excitation operators using multi-indices $\mu>0$ that clarify the projections and wedge products involved ^{12, 62}. Therefore, acting on the reference state $|{\Phi_{0}}\rangle$, the elementary particle-hole excitation operators define the Hilbert space of functions that are $L^{2}$-orthogonal to $|{\Phi_{0}}\rangle$. We henceforth set $\mathcal{V}:=\{|{\Phi_{0}}\rangle\}^{\perp}$ and note that $\hat{T}({\bf t})|{\Phi_{0}}\rangle\in\mathcal{V}$ for all CC amplitudes $\mathbf{t}$. Hence, we can express Eq. (4) in the more common form as a set of projective equations, i.e.

$$0=\langle\Phi_{\mu}|e^{-\hat{T}({\bf t})}\hat{H}e^{\hat{T}({\bf t})}|\Phi_{0}\rangle,\quad\forall\mu>0,$$ (6)

where $|{\Phi_{\mu}}\rangle=\hat{X}_{\mu}|{\Phi_{0}}\rangle$ are the excited Slater determinants spanning $\mathcal{V}$. In practice, the expansion in Eq. (5) is commonly limited to e.g. only containing one and two particle-hole excitations leading to CCSD. In this case, the projective equations (6) are also restricted to the corresponding set of excited Slater determinants, yielding a square system of polynomial equations that we seek to solve. Hence, a key object in CC theory is the set

$$\mathcal{S}=\{{\bf t}\in\mathbb{F}^{K}~{}|~{}\langle\Phi_{\mu}|e^{-\hat{T}({\bf t})}\hat{H}e^{\hat{T}({\bf t})}|\Phi_{0}\rangle=0,\quad\forall\mu>0\},$$ (7)

where $\mathbb{F}$ is the considered number field (either $\mathbb{R}$ or $\mathbb{C}$) and $K$ is the “system size” (determined by the number of correlated electrons, the size of one-particle basis functions as well as further selection rules). Note that $\mathbb{F}$ can be $\mathbb{R}$ or $\mathbb{C}$ depending on whether we are seeking real- or complex-valued CC amplitudes. We emphasize that although the CC polynomial coefficients are real, the roots of the polynomial system may not be. Mathematically, the theory describing complex-valued solutions is simpler and more complete than the theory describing real-valued solutions. Therefore, we shall consider $\mathbb{F}=\mathbb{C}$ for mathematical considerations, however, solving the truncated CC equations over $\mathbb{F}=\mathbb{C}$ may yield unwanted complex valued energies.

Mathematically, the set $\mathcal{S}$ is an (algebraic) variety. Since the concept of varieties may be unconventional within the computational chemistry community, we provide an illustrative example. A variety is a mathematical object defined by the zeros of a set of polynomial equations. Consider the multivariate polynomial system

$$\left\{\begin{aligned} p_{1}(x,y,z)&=x^{2}+y^{2}-z^{2},\\ p_{2}(x,y,z)&=\frac{x+y}{\sqrt{2}}+z-1.\end{aligned}\right.$$ (8)

The zeros of $p_{1}$ form a cone, whereas the zeros of $p_{2}$ form a plane; the set of simultaneous roots defines a variety that corresponds to a parabolic curve, see Fig. 1. As a curve, this variety is one-dimensional, embedded in the ambient space which is three-dimensional in this case. A central quality for this work is the degree of a variety which is defined to be the maximal number of intersection points of the variety with a general linear space of complementary dimension. By requiring the linear space to be general we ensure the intersection is always finite. For the variety in Fig. 1 we see that intersecting the yellow curve with a plane (a linear space of complimentary dimension) yields at most two intersection points, hence, its degree is two. We add that over the complex field, the condition of “maximal number of intersection points” is redundant for the definition of degree since the number of intersection points is the same for generic complex spaces of complementary dimension.

Figure 1: Root structure corresponding to Eq. (8). The cone describes the roots of $p_{1}$, and the plane describes the roots of $p_{2}$. The simultaneous roots are described by the yellow curve that defines the solution variety of the multivariate polynomial system.

We conclude this section by noting that in the case of untruncated CC, amplitudes that solve the full CC equations describe eigenstates of the Hamiltonian. Therefore, the cardinality of $\mathcal{S}$ is exactly the number of eigenstates of $\hat{H}$ that are intermediately normalized and the number of roots to the CC equations is bounded from above by the number of Slater determinants. However, when truncations are imposed, the cardinality of $\mathcal{S}$ increases, i.e. truncated CC theory yields (some) unphysical solutions. Therefore, it becomes less clear what the different elements in $\mathcal{S}$ describe. Understanding the variety $\mathcal{S}$ and characterizing (some of) its elements are the subject of this manuscript.

3 Truncation varieties

In CC theory one encounters various algebraic equations that define different varieties. The truncation varieties ⁶⁰ arise from the exponential parametrization, i.e.

$${\rm exp}:\mathbb{V}\to\mathcal{V},~{}\mathbf{t}\mapsto e^{\hat{T}(\mathbf{t})}|\Phi_{0}\rangle,$$ (9)

where $\mathbb{V}$ is the CC amplitude space. In the untruncated CC method, the exponential map is bijective. However, when truncations are introduced, this property no longer holds. Truncations correspond to restrictions of the exponential map to a subspace, e.g. for CCSD the exponential map is restricted to $\mathbb{V}_{\{1,2\}}\subset\mathbb{V}$ consisting of cluster amplitudes with only one and two particle-hole excitations. The corresponding truncation variety, $V_{\{1,2\}}$ is then defined as the smallest variety containing ${\rm exp}(\mathbb{V}_{\{1,2\}})$. This construction of the truncation variety can be generalized to any CC truncation: Let $\sigma$ denote the CC truncation rank, i.e. $\sigma\subset[\![N]\!]=\{1,\dots,N\}$, where $N$ denotes the number of electrons in the system. Then, the smallest variety containing ${\rm exp}(\mathbb{V}_{\sigma})$ is the $\sigma$-truncation variety denoted $V_{\sigma}$.

Remarkably, the truncation varieties can be used to express the projective CC equations (6) as an eigenvalue problem, i.e.

$$[\hat{H}|\Phi\rangle]_{\sigma}=E[|\Phi\rangle]_{\sigma},\quad|\Phi\rangle\in V_{\sigma},$$ (10)

where $[\,\cdot\,]_{\sigma}$ is the projection onto coordinates with excitation level in $\sigma\cup\{0\}$, i.e.

$$[|\Psi\rangle]_{\sigma}=\langle\Phi_{\mu}|\Psi\rangle\quad{\rm for}~{}|\mu|\in\sigma\cup\{0\}.$$ (11)

This formulation is indeed equivalent to the projective CC equations in Eq. (6) for CCD, as well as for rank complete CC truncations, i.e. CCSD, CCSDT, CCSDTQ etc., see Theorem 5.11 in Ref. 60. The number of roots for these CC variants is bounded by the number of solutions to Eq. (10) for a generic Hamiltonian matrix $\hat{H}$. We denote this number the CC degree of the corresponding truncation variety. The CC degree serves as a complexity measure for solving the CC equations and Bézout’s theorem provides the upper bound

$$\mathrm{CCdeg(\sigma)}\leq(\dim(V_{\sigma})+1)\deg(V_{\sigma}),$$ (12)

where $\dim(V_{\sigma})$ is the dimension of the truncation variety and $\deg(V_{\sigma})$ is its degree. The dimension of the truncation variety is well understood. The reason is that cluster analysis allows us to bijectively map cluster amplitudes to CI-coefficients using polynomial equations, see e.g. Refs. 12, 70, 32 or Proposition 2.4 in Ref. 60 for a mathematical proof. Then, since $V_{\sigma}$ arises as a deformation of the linear space $\mathbb{V}_{\sigma}$ under a polynomial with polynomial inverse, the dimension of $V_{\sigma}$ is equal to the dimension of $\mathbb{V}_{\sigma}$. For example in the case of CCSD, we obtain

$$\dim(V_{\sigma})=\dim(\mathbb{V}_{\sigma})=n_{\rm virt}n_{\rm occ}+\binom{n_{\rm virt}}{2}\binom{n_{\rm occ}}{2}.$$ (13)

The degree of ${V}_{\sigma}$ on the other hand is not as easily computable and requires a more profound understanding of ${V}_{\sigma}$. For special cases, however, explicit formulas for the degree have been established; for example, for CCS, the truncation variety is the Grassmanian (see Theorem 3.5 in Ref. 60) which is mathematically well understood. In other cases we rely either on symbolic methods, provided by e.g. Macaulay2 ⁷¹, or on numerical methods, accessible in e.g. HomotopyContinuation.jl ⁷², to compute the degree.

Despite the complexity of calculating $\deg(V_{\sigma})$ the bound given by Eq. (12) provides significant insight into bounding the CC degree, since it is less complex to calculate $\deg(V_{\sigma})$ than to calculate the CC degree directly. Moreover, this upper bound is notably better than the previously known upper bounds obtained by applying Bézout’s theorem directly to the projective CC equations, see Section 6 in Ref. 60.

4 Computational methods

In this section, we outline the computational procedure employed to compute the CC degree and the CC roots for a given CC truncation rank $\sigma\subset[\![N]\!]$. We begin by introducing the CC family, a parametric polynomial family for Eq. (10), i.e.

$$\mathcal{F}_{\rm CC}(\sigma)=\{[(\hat{G}-E\cdot I_{K})e^{\hat{T}(\mathbf{t})}]_{\sigma}~{}\big{|}~{}\hat{G}\in\mathbb{F}^{K\times K}\},$$ (14)

where the energy $E\in\mathbb{C}$ and cluster amplitudes $\mathbf{t}\in\mathbb{V}_{\sigma}$ are the variables, and the entries of the matrix $\hat{G}$ are the parameters. We will ultimately set $\hat{G}=\hat{H}$ to be the Hamiltonian of interest, but $\mathcal{F}_{\rm CC}(\sigma)$ is a much broader family of CC-like systems with $\hat{G}\in\mathbb{F}^{K\times K}$ an arbitrary matrix. The parameter continuation theorem provides an upper bound to the number of roots of any polynomial system $F_{\hat{G}}\in\mathcal{F}_{\rm CC}$. Moreover, this bound is tight and equality is obtained for all polynomial systems with parameters in a certain dense subspace of $\mathbb{F}^{K\times K}$. We subsequently assume the validity of the parameter continuation theorem and refer the interested reader to Section 3 in Ref. 73. The upper bound provided by the parameter continuation theorem yields the CC degree mentioned in the previous section.

The numerical procedure used in this work consists of two steps: First, we solve a general system of equations from the CC family, $\mathcal{F}_{\rm CC}$, using monodromy methods. Second, we use parameter homotopy continuation to track the solutions from this solved but generic system of equations to solutions of the CC equations describing a specific chemical system of interest. Both steps are subsequently outlined in detail.

We emphasize that, at their core, both methods use a homotopy to track new solutions from known solutions ^{74, 75, 76, 77}. More precisely, we define a homotopy for a (piecewise smooth) path $\gamma:[0,1]\to\mathbb{F}^{K\times K}$ in the parameter space via

$$H(E,\mathbf{t},\lambda)=F_{\gamma(\lambda)}(E,\mathbf{t}),$$ (15)

where $H(E,\mathbf{t},\lambda)$ describes the continuous deformation of a start system $G(E,\mathbf{t})=H(E,\mathbf{t},1)$ to a target system $F(E,\mathbf{t})=H(E,\mathbf{t},0)$. The individual solution paths to the system $H(E,\mathbf{t},\lambda)$ are then tracked from $G(E,\mathbf{t})$ to $F(E,\mathbf{t})$ as $\lambda\to 0$. This is accomplished via an ordinary differential equation, known as the Davidenko differential equation ^{78, 79}, i.e.

$$\frac{\partial}{\partial\mathbf{x}}H(\mathbf{x},\lambda)\left(\frac{\mathrm{d}}{\mathrm{d}\lambda}\mathbf{x}(\lambda)\right)+\frac{\partial}{\partial\lambda}H(\mathbf{x},\lambda)=0,$$ (16)

initialized by a root $\mathbf{x}$ of the start system, i.e. $G(\mathbf{x})=0$ where we set $\mathbf{x}=(E,\mathbf{t})\in\mathbb{C}\times\mathbb{V}_{\sigma}$.

4.1 Monodromy

The monodromy solver is a numerical technique specifically designed for solving generic polynomial systems within a parameterized family. We emphasize that the monodromy solver requires the system to be generic, meaning it is not suitable for solving a specific set of CC equations directly. However, it can be employed to establish a start system $F_{\hat{G}}\in\mathcal{F}_{\rm CC}$ within the CC family. This start system is then used by the parameter homotopy continuation method to solve a (physical) target system of interest. Note that the monodromy method is required only once for a particular system configuration, i.e. the number of electrons, the number of orbitals, and the CC truncation level $\sigma\in[\![N]\!]$.

We start with a root $\mathbf{x}\in\mathbb{C}\times\mathbb{V}_{\sigma}$ of a generic system $F_{\hat{G}}\in\mathcal{F}_{\rm CC}$ in the CC family and a loop $\gamma:[0,1]\to\mathbb{F}^{K\times K}$ in the parameter space, based at the generic matrix $\hat{G}$, i.e. $\gamma(0)=\gamma(1)=\hat{G}$. Moreover, we let $\mathbf{x}(\lambda)$ be the solution path of $H(E,\mathbf{t},\lambda)=0$ along $\gamma$, starting at $\mathbf{x}(1)=\mathbf{x}$. Since $\gamma$ is a loop, the start system in this homotopy is equal to the end system. In particular, the solution path ends again at a root of $F_{\hat{G}}$. The monodromy solver exploits the fact that the root $\mathbf{x}(0)$ need not necessarily be equal to the root $\mathbf{x}(1)$. In fact, let $(\mathbf{x}_{1},\dots,\mathbf{x}_{M})$ be all roots to $F_{\hat{G}}$, then the roots $(\mathbf{x}_{1}(0),\dots,\mathbf{x}_{M}(0))$ define a mere permutation the system’s roots.

We exemplify this procedure with the simple parametric polynomial family:

$$\mathcal{F}=\{x^{3}-6\,x^{2}+11\,z\,x-6~{}\big{|}~{}z\in\mathbb{C}\},$$ (17)

where $x\in\mathbb{C}$ is the variable and $z\in\mathbb{C}$ is the parameter. For $z=1$, the corresponding polynomial

$$x^{3}-6\,x^{2}+11\,x-6=(x-1)(x-2)(x-3)$$ (18)

has the roots $x_{1}=1,x_{2}=2,x_{3}=3$. The monodromy solver then tracks these roots along the path $\gamma:[0,1]\to\mathbb{C},t\mapsto e^{i2\pi t}$. Note that $\gamma$ is a loop based at $1$. Indeed, this procedure permutes the roots, see Figure 2.

Figure 2: The solution paths of the roots of $x^{3}-6x^{2}+11zx-6$ for $z(\lambda)=\gamma(\lambda)$ along the path $\gamma:[0,1]\to\mathbb{C},t\mapsto e^{i2\pi t}$.

This approach can be effectively utilized to find all solutions to a generic system of equations. For instance, consider a scenario where the only known solution for the polynomial in Eq. (18) is $x_{2}=2$. By tracking this solution along the unit circle, we can determine another solution, $x_{1}=1$. When we extend this tracking across sufficiently many loops, the method reveals all possible roots of the system.

This process is central to the monodromy solver. The solver begins with a single known solution and numerically traces known solutions through multiple loops within the parameter space. This tracking continues until no additional solutions are found over a predetermined number of loops, leading to the conclusion that all possible solutions to the system have been numerically identified.

The parameters in the CC family $\mathcal{F}_{\rm CC}$ are linear, hence, we can find a generic start system along with a root by choosing a pair $(E,\mathbf{t})\in\mathbb{C}\times\mathbb{V}_{\sigma}$ at random, and then compute one generic solution $\hat{G}_{0}$ fulfilling

$$[(\hat{G}_{0}-E\cdot I_{K})e^{\hat{T}(\mathbf{t})}]_{\sigma}=0.$$

We then use the monodromy solver to compute all roots to $F_{\hat{G}_{0}}$.

4.2 Parametric homotopy continuation

Applying the parametric homotopy continuation method to find the roots of the CC equations corresponding to a system of interest is now straightforward. The goal is to solve

$$F_{\hat{H}}(E,{\mathbf{t}})=[(\hat{H}-E\cdot I_{K})e^{\hat{T}(\mathbf{t})}]_{\sigma}=0$$ (19)

for a given Hamiltonian $\hat{H}$, arising from a physical electronic system. By construction, both $\hat{H}$ and $\hat{G}_{0}$ are in the parameter space of $\mathcal{F}_{\rm CC}$. In particular, we can define a path $\gamma$ in the parameter space that continuously transforms $\hat{G}_{0}$ into $\hat{H}$. Since we have computed all roots to the system $F_{\hat{G}_{0}}$, the Davidenko differential equation (16) allows the tracking of solutions of $F_{\hat{G}_{0}}$ to $F_{\hat{H}}$ using e.g. Newton’s method.

As $\lambda\to 0$, we may encounter different scenarios when getting close to $\lambda=0$. Here, the endgame ⁸⁰ should be employed to determine the target solutions and their “type”. We illustrate such different scenarios in Figure 3.

Figure 3: Sketch of possible homotopy paths. The solid line shows a path with no finite limit as $\lambda\to 0$, the dashed lines have the same limit of a singular root, and the dotted-dashed line has a unique limit of a nonsingular root.

We conclude this section by highlighting the progress that has been made within homotopy continuation based software development exploiting parallel implementations, which can significantly extend their application areas in the coming years. In particular, we here want to emphasize PHCpack ⁸¹, Bertini, ⁸² HOM4PS, ^{83, 84} NAG4M2 ⁷⁴ and HomotopyContinuation.jl. ⁷²

5 Numerical results

In this section we investigate the root structure and potential energy curves (PECs) for ground and excited states at the level of CCD and CCSD. We report PECs corresponding to different CC solutions obtained by the above-outlined computational procedure. The obtained energies and corresponding eigenstate approximations are then compared to exact eigenpairs of the respective Hamiltonian. Besides the computation and analysis of the CC states, we compare different bounds to the number of CC roots with the computed number of CC solutions along different bond stretching processes. We restrict our investigations to four electron systems, such as (H₂)₂ in D_2h and D_∞h configurations, H₄ symmetrically disturbed on a circle ⁸⁵, and LiH. The reported computations are performed in first quantization and the considered systems are all at half-filling, i.e. four electrons in eight spin orbitals. The one-body and two-body integrals for the considered systems are obtained from the Python-based Simulations of Chemistry Framework (PySCF) ^{86, 87, 88}. For the algebraic simulations, we used the monodromy and parametric homotopy continuation methods in HomotopyContinuation.jl. ⁷²

5.1 CC bounds for generic four electron systems at half-filling

We emphasize that the theoretical bounds hold for generic Hamiltonians, and therefore only depend on the number of electrons, the number of spin orbitals, and the CC truncation level.

Coupled cluster doubles:

The roughest bound is obtained by applying Bézout’s theorem directly to the CCD equations. Note that for four electron systems, the CCD equations are quadratic. Moreover, in first quantization, there are ${4\choose 2}^{2}=36$ doubly excited Slater determinants. Hence, the Bézout bound directly applied to the CCD equations yields $2^{36}$ as an upper bound for the number of roots. The second bound we compare with is the bound in Eq. (12). Again, in first quantization, there are 36 CCD amplitudes, hence, ${\rm dim}(V_{\{2\}})=36$. Moreover we compute the degree of $V_{\{2\}}$ using Macaulay2 yielding ${\rm deg}(V_{\{2\}})=2$. Overall, this yields $74$ as an upper bound. For a generic Hamiltonian, we can moreover compute the exact number of roots to the CCD equations using the monodromy method which yields CCdeg${}_{4,8}(\{2\})=73$. This quantifies the improvement provided by the bound in Eq. (12) over the direct application of Bézout’s theorem.

Coupled cluster singles and doubles:

When considering CCSD instead of CCD these numbers dramatically increase. Frist, we note that the Bézout bound directly applied to the CCSD equation yields the upper bound $3^{n_{s}}4^{n_{d}}$, where for four electrons in eight spin orbitals $n_{s}=16$ and $n_{d}=36$. This leads to approximately $2\cdot 10^{29}$ solutions. Second, we note that ${\rm dim}(V_{\{1,2\}})=52$ and numerical calculations yield ${\rm deg}(V_{\{1,2\}})=442066$, hence, the bound provided by Eq. (12) yields approximately $2\cdot 10^{7}$ solutions. Numerical calculations show that CCdeg${}_{4,8}(\{1,2\})=16\,952\,996$. This clearly illustrates the improvement provided by the bound in Eq. (12) over the direct application of Bézout’s theorem.

To connect these theoretical bounds to practical systems, we subsequently calculate the number of CC solutions that result in real-valued energies, as well as those that accurately predict excited state energies at different bond lengths for the considered systems. Additionally, we report the number of FCI solutions. This comparison allows us to assess the accuracy of the bounds established for generic Hamiltonians.

5.2 Dissociation of lithium hydride – CCD

We initiate our energy and state approximation investigations by looking at the dissociation process of lithium hydride. We consider bond distances ranging from $1.375$ to $5.95$ bohr. Diagonalizing the Hamiltonian yields the full spectrum reported in Figure 16(a) (see supplementary material) together with all roots to the CCD equations. Due to the large number of FCI solutions, we minimize the distance between FCI solutions and CCD solutions showing a much clearer picture of potentially well-resolved PECs, see Figure 16(b) in the supplementary material. This graph allows us to extract the PECs of the different states that are well-approximated using CCD, see Figure 4(a). We find that 10 PECs are approximated up to $5\cdot 10^{-2}$ hartree for the entire dissociation process, however, we also note that for various bond distances, more CCD energies approximate FCI energies well, see Figure 16(b) and Table 1. For the 10 PECs, we moreover compute the overlap between the CCD states and the approximated eigenstates, see Figure 4(b). Remarkably, three of the four lowest energy states are well approximated both in terms of energy and states. The highest energy state stands out as well since it is well-approximated near the equilibrium, both in terms of energy as well as the eigenstate. However, as we transition away from the equilibrium the level of approximation deteriorates. Note that poor overlap may be detected when the eigenspace is higher dimensional, that is if the corresponding eigenvalue has multiplicity two or higher. Therefore some of the higher energy states reporting bad overlap might have CCD states that are close to the actual eigenspace – or even within the eigenspace – since we are only reporting overlap with one representative from the eigenspace. A detailed investigation of this is left for future works.

(a)

(b)

Figure 4: (a) PECs of LiH that can be accurately approximated by CCD energies. (b) Overlaps of the CCD states with the corresponding eigenstate.

Investigating the root structure, we observe that the total number of CCD roots fluctuates along the bond stretching process, see Table 1. We find that any bound for a generic Hamiltonian severely overestimates the true number of roots for this system.

Bond distances	1.375	1.6	1.9	2.2	2.5	3.1	4.0	4.9	5.8
$\#$ CCD real	33	40	33	40	34	28	30	33	28
$\#$ CCD approx.	22	25	20	23	18	19	17	18	14

Table 1: The number of CCD roots in LiH dissociation for selected bond distances. By "$\#$ CCD approx." we denote the energetically relevant CCD roots. For comparison, we recall that the Bézout bound led to $2^{36}$ solutions, the bound in Eq. (12) led to 74 solutions, and CCdeg_4,8({2}) is 73.

5.3 Lithium hydride – CCSD

We moreover perform CCSD computations on lithium hydride at the bond distance 2.875 bohr – close to the molecule’s equilibrium. Computing all roots to CCSD equations yields $15\,954$ roots, $2\,170$ of which yield real-values energies, and $1\,280$ are real-valued solutions. Visualizing all $15\,954$ CCSD energies shows that the energies tend to cluster around the FCI solutions, see Figure 5(a). Being able to compare with the true Hamiltonian spectrum reveals that only 26 of these CCSD solutions yield energies that are close to an FCI energy up to $10^{-3}$ hartree, see Figure 5(b).

(a)

(b)

Figure 5: (a) FCI solutions (red lines) together with all CCSD energies (b) CCSD energies that approximate FCI energies up to $10^{-3}$ hartree.

For all 70 eigenstates, we compute the energetically closest CCSD state together with the overlap with the corresponding states, see Figure 6. We find that three of the 26 energetically relevant CCSD states have good overlaps with the targeted states. Therefore there are at least three CCSD states that are well approximated both in terms of energy and state.

Figure 6: Overlap of the energetically closest CCSD state with the corresponding eigenstate. The blue dots show the 26 energetically relevant CCSD states.

5.4 Dissociating systems of hydrogen

We now investigate the dissociation of (H₂)₂ planar model systems in different geometries ³⁸.

5.4.1 (H₂)₂ dissociation (D_2h symmetry)

The bond stretching procedure for (H₂)₂ in D_2h symmetry is sketched in Figure 7. We consider an intra-molecular bond distance of $R=1.4$ bohr, which corresponds to the equilibrium geometry of H₂; this is kept fixed during the dissociation process. The inter-molecular distance is varied from 1.5 to 4.0 bohr.

Figure 7: Schematic depiction of the dissociation process of (H₂)₂ in D_2h configuration. The parameter $R=1.4$ bohr.

The full spectrum together with the CCD solutions can be found in Figure 13(a) in the Supplementary material, and the PECs together with the closest CCD solutions are reported in Figure 13(b). We extract the PECs that are well-approximated using CCD and report them in Figure 8(a). Here, we observe that eight CCD energies approximate FCI PECs well. For these eight CCD states, we compute the overlap with the corresponding eigenstate, see Figure 8(b). We see that the ground state is well approximated, both in terms of energy as well as the state. The following six states (i.e. states two to seven in Figure 8(b)) have intermediate overlaps with the targeted states; the last state has a poor overlap with the targeted state. This illustrates that CC amplitudes could yield good energy approximations while providing poor approximations to the actual eigenstates ¹³.

(a)

(b)

Figure 8: (a) PECs of (H₂)₂ in D_2h symmetry that are accurately described by CCD energies. (b) Overlap of the CCD states with the corresponding eigenstate.

In Table 2 we report the number of CCD roots yielding real-valued energies and energetic relevant CCD roots for selected bond distances. We observe that the number of CCD roots fluctuates along the bond stretching procedure.

Bond distances	1.6	1.9	2.1	2.5	3.0	3.9
$\#$ CCD real	66	45	60	64	55	64
$\#$ CCD approx.	32	26	33	32	33	34

Table 2: The number of roots along the (H₂)₂ dissociation in D_2h configuration for selected bond distances. By "$\#$ CCD approx." we denote the energetically relevant CCD roots. For comparison, we recall that the Bézout bound is $2^{36}$, the bound in Eq. (12) yields 74, and CCdeg_4,8({2}) is 73.

5.4.2 (H₂)₂ dissociation (D_∞h symmetry)

Next, we investigate the (H₂)₂ dissociation in D_∞h symmetry, see Figure 9 for a sketch of the dissociation process. The intra-molecular distance is again set to $R=1.4$ bohr and is kept fixed during the dissociation process. The inter-molecular distance is varied from 1.5 to 4.0 bohr.

Figure 9: Schematic depiction of the dissociation process of (H₂)₂ in D_∞h configuration. The parameter $R=1.4$ bohr.

The full spectrum and CCD solutions is presented in Figure 14(a) in the Supplementary material, with the PECs and their closest CCD solutions detailed in Figure 14(b). We have identified PECs that are accurately approximated by CCD, as shown in Figure 10(a), where four CCD energies closely match the FCI PECs. The overlap of these four CCD states with their respective eigenstates is analyzed in Figure 10(b). Consistent with observations from the (H₂)₂ in D_2h configuration, the ground state demonstrates a high degree of accuracy, while the other three states exhibit intermediate to poor overlap with their targeted eigenstates.

(a)

(b)

Figure 10: (a) PECs of (H₂)₂ in D_∞h symmetry that are accurately described by CCD energies. (b) Overlap of the CCD states with the corresponding eigenstate.

Investigating the root structure, we observe that – similar to (H₂)₂ in D_2h symmetry – the total number of CCD roots fluctuates along the bond stretching procedure, see Table 3. Note that for (H₂)₂ in D_∞h symmetry, the total number of CCD roots is lower than the number of CCD roots (H₂)₂ in D_2h symmetry. In particular, any bound for a generic Hamiltonian severely overestimates the true number of roots for this system.

Bond distances	1.6	1.9	2.1	2.5	3.0	3.9
$\#$ CCD real	46	42	40	38	37	35
$\#$ CCD approx.	13	19	12	14	16	21

Table 3: The number of roots along the (H₂)₂ dissociation in D_∞h configuration for selected bond distances. By "$\#$ CCD approx." we denote the energetically relevant CCD roots. For comparison, we recall that the Bézout bound led to $2^{36}$ solutions, the bound in Eq. (12) led to 74 solutions, and CCdeg_4,8({2}) is 73.

5.5 H₄ disturbed on a circle

We proceed by investigating a variant of the H₄ model consisting of four hydrogen atoms symmetrically distributed on a circle of radius $R=\sqrt{2}$ bohr ⁸⁹, see Fig. 11.

Figure 11: Schematic depiction of the H₄ model undergoing a symmetric disturbance on a circle modeled by the angle $\Theta$.

The full spectrum together with the CCD solutions can be found in Figure 15(a) in the Supplementary material, and the PECs together with the closest CCD solutions are reported in Figure 15(b). We extract the PECs that are well-approximated using CCD which are reported in Figure 12(a). Note that it is well-known that for this system CCD does not perform well in the region close to 90^∘, due to strong degeneracies ³⁸. This results in a dramatic reduction of the total number of roots, we find three roots at $90^{\circ}$ and six roots at $89^{\circ}$. For completion, these energies are included in the supplementary material, however, since the focus of this work is on the root structure and potentially physically relevant roots accessible by single-reference CC theory, we focus on cases where single-reference CC theory provides reasonable approximations and therefore excluded these points in Figure 12(a). Outside of this challenging region, we find five solutions that approximate PECs well. For these five solutions, we compute the overlap with the corresponding eigenstate, see Figure 12(b). Consistent with the previous H₄ systems, this model system shows the effect of CCD states accurately resolving excited energies but providing poor overlap to the exact eigenstates. The only state that shows a good approximation to its eigenstate is the CCD ground state.

(a)

(b)

Figure 12: (a) PECs of H₄ symmetrically distributed on a cycle that can be accurately reconstructed by CCD roots. (b) Overlap of the CCD states with the corresponding eigenstate.

Similar to (H₂)₂ in D_∞h symmetry, we observe that generic bounds severely overestimate the number of CCD roots, though this effect is further amplified for H₄ symmetrically distributed on a circle, see Table 4.

Angles in degrees	45	55	65	75	85
$\#$ CCD real	35	36	35	36	28
$\#$ CCD approx.	18	18	19	21	22

Table 4: The number of roots of H₄ symmetrically distributed on a circle for selected bond distances. By "$\#$ CCD approx." we denote the energetically relevant CCD roots. For comparison, we recall that the Bézout bound led to $2^{36}$ solutions, the bound in Eq. (12) led to 74 solutions, and CCdeg_4,8({2}) is 73.

6 Conclusion

The introduction highlighted critical unanswered questions in the study of the coupled cluster (CC) equations’ algebraic structures, including the total number of solutions, the best approximation element, and the feasibility of approximating excited states with higher-order roots of single reference CC equations. The ability to compute the full solution spectrum to the CC equations enabled us to tackle these questions for concrete physical systems.

Most remarkably, our results demonstrate that for lithium hydride in multiple geometries, higher-order CC solutions not only provide high-accuracy approximations to the energies but also to the states themselves. For a broader range of systems including multiple (H₂)₂ planar model systems in different geometries, we found that high-order CC solutions accurately predict several full potential energy curves, though they possibly provide a poor approximation to the eigenstates. While this effect has been previously reported, it requires further investigation to fully understand it and ultimately leverage it for practical applications.

Additionally, we examined the total number of roots and the effectiveness of various bounds developed for generic Hamiltonians. Our simulations indicate that while the most recent bounds perform exceptionally well for generic systems, they significantly overestimate the number of solutions in physical systems. More precisely, for all considered systems, we computed between $30$-$60$ real solutions of which only $10$-$30$ approximate physical states, compared to CC degree of $72$. This excess of solutions gets amplified for CCSD: Solving the CC equations of a generic Hamiltonian yields $16\,952\,996$ solutions, compared to $1280$ solutions that yield real-valued energies of the CC equations describing LiH. This overestimation by a factor of over ten thousand indicates that to further improve the bounds of the CC equations, specific structures of the Hamiltonian have to be taken into account.

The authors thank Piotr Piecuch for insightful discussions.

7 Dissociating systems of hydrogen

We begin by reporting the full energy spectrum of the Hamiltonian describing (H₂)₂ in D_2h symmetry together with all CCD roots, see Figure 13(a). Given the large number of solutions, we minimize over the difference between spectral values of the Hamiltonian and the found CCD energies. This provides further insight into energies that are well approximated by CCD, see Figure 13(b).

(a)

(b)

Figure 13: (a) Full spectrum (b) Minimal difference.

Similarly, we proceed for (H₂)₂ in D_∞h symmetry. We reporting the full energy spectrum of the Hamiltonian describing (H₂)₂ in D_∞h symmetry together with all CCD roots, see Figure 14(a) Given the large number of solutions, we minimize over the difference between spectral values of the Hamiltonian and the found CCD energies. This provides further insight into energies that are well approximated by CCD, see Figure 14(b).

(a)

(b)

Figure 14: (a) Full spectrum (b) Minimal difference.

8 H₄ disturbed on a circle

Similarly, we proceed for H₄ disturbed on a circle. We report the full energy spectrum of the Hamiltonian describing H₄ disturbed on a circle together with all CCD roots, see Figure 15(a). Given the large number of solutions, we minimize over the difference between spectral values of the Hamiltonian and the found CCD energies. This provides further insight into energies that are well approximated by CCD, see Figure 15(b).

(a)

(b)

Figure 15: (a) Full spectrum (b) Minimal difference.

9 Lithium hydride

We report the full energy spectrum of the Hamiltonian describing LiH together with all CCD roots, see Figure 16(a) Given the large number of solutions, we minimize over the difference between spectral values of the Hamiltonian and the found CCD energies. This provides further insight into energies that are well approximated by CCD, see Figure 14(b).

(a)

(b)

Figure 16: (a) Full spectrum (b) Minimal difference.