Algebraic Varieties in Quantum Chemistry

We develop algebraic geometry for coupled cluster (CC) theory of quantum many-body systems. The high-dimensional eigenvalue problems that encode the electronic Schr¨odinger equation are approximated by a hierarchy of polynomial systems at various levels of truncation. The exponential parametrization of the eigenstates gives rise to truncation varieties. These generalize Grassmannians in their Pl¨ucker embedding. We explain how to derive Hamiltonians, we offer a detailed study of truncation varieties and their CC degrees, and we present the state of the art in solving the CC equations.


Introduction
Electronic structure theory is a powerful quantum mechanical framework for investigating the intricate behavior of electrons within molecules and crystals.At the core lies the interaction between particles, specifically the electron-electron and electron-nuclei interactions.Embracing the essential quantum physical effects, this theory is the foundation for ab initio electronic structure calculations performed by many researchers in chemistry and related fields, complementing and supplementing painstaking laboratory work.With its diverse applications in chemistry and materials science, electronic structure theory holds vast implications for the mathematical sciences.Integrating methods from algebra and geometry into this field leads to the development of precise and scalable numerical methods, enabling extensive in silico studies of chemistry for e.g.sustainable energy, green catalysis, and nanomaterials.The synergy between fundamental mathematics and electronic structure theory offers the potential for groundbreaking advancements in addressing these global challenges.
The electronic structure problem is the innocent-looking eigenvalue problem in (18).This has been under intense investigation since the 1920s, but it is still a formidable problem of contemporary science.The governing partial differential equation, i.e. the electronic Schrödinger equation (19), has 3d degrees of freedom, where d is the number of electrons [29].As d increases, the problem size of (18) grows exponentially.Efficient and tractable numerical schemes are essential for approximating the behavior of complex atoms and molecules [21].One class of widely employed high-accuracy methods rests on coupled cluster theory [15], which is considered the gold standard of quantum chemistry.Despite considerable successes, this approach has its limitations, leaving many opportunities for further developments.
In this article, we develop algebraic geometry for coupled cluster theory, following [10].Central to our investigations are the coupled cluster (CC) equations.This system of polynomial equations had already been studied by the quantum chemistry community in the 1990s.Researchers used homotopy continuation for numerical solutions, and they focused on enhancing computational capabilities.For details on this history we refer the interested reader to [17, and the references therein, or [10, Sections 1.1 and 1.2].In contrast, our present study goes beyond advancements in computations, delving deeply into the mathematical structures that underlie CC theory.It rests on nonlinear algebra [1,2,22,27].
On the geometry side, our point of departure is the Grassmannian Gr(d, n) with its n d Plücker coordinates.These represent quantum states for d electrons in n spin-orbitals.We introduce a family of projective varieties V σ in the same ambient space P ( n d )−1 , one for each subset σ of [d] = {1, 2, . . ., d}.The singleton σ = {1} yields the Grassmannian Gr(d, n).
We formulate the coupled cluster equations as a truncated eigenvalue problem on V σ .The truncation varieties V σ for σ = {1}, {2}, {1, 2}, {1, 2, 3} correspond to the CC variants CCS, CCD, CCSD, CCSDT, explained in e.g.[9,10].The number of complex solutions to the CC equations for a general Hamiltonian H is the CC degree, denoted CCdeg d,n (σ).This invariant reveals the number of paths that need to be tracked for finding all solutions.
Experts in nonlinear algebra will find these concepts to be consistent with notions they are familiar with.The truncated eigenvalue problem on V σ is reminiscent of the theory of eigenvectors for tensors [22,Section 9.1].The CC degree is an analog to the Euclidean distance (ED) degree [27, Section 2] and to the maximum likelihood (ML) degree [27,Section 3].
The present article launches an entirely new line of research.By contrast, previous mathematical investigations of CC theory were performed within a functional analytic framework.In that framework, Schneider performed the first local analysis of CC theory in 2009 which was based on Zarantonello's lemma [26].This approach was subsequently extended [24,25] and applied to different CC variants [18,19,11].Later Csirik and Laestadius established a more versatile framework for analyzing general CC variants which uses topological degree theory [6,7].Recent numerical analysis results regarding single reference CC were established by Hassan, Maday, and Wang using the invertibility of the CC Fréchet derivative [14].
We now summarize the organization and contributions of this paper.In Section 2 we present the exponential parametrization which expresses the quantum states in terms of the cluster amplitudes.This map is invertible.Theorem 2.5 gives formulas for all coordinates of the forward map and the backward map.This involves a master polynomial of degree d whenever n = 2d.For instance, for d = 2, this polynomial is the Plücker quadric ψ 12 ψ 34 − ψ 13 ψ 24 + ψ 14 ψ 23 .In general, its monomials correspond to uniform block permutations [23].
Section 3 introduces the truncation variety V σ which lives in the projective space P ( n d )−1 .Its restriction to an affine chart is a complete intersection, defined by some coordinates of the backward map.The homogeneous prime ideal of V σ is found by saturation (Theorem 3.2).Theorem 3.5 reveals the Grassmannian for σ = {1}.Proposition 3.7 features particle-hole symmetry (d, n) ↔ (n−d, n), and Theorem 3.10 identifies all σ for which V σ is a linear space.
Section 4 starts from high school chemistry.It explains the discretization process in coupled cluster theory.From the electronic Schrödinger equation we derive the Hamiltonian H, a symmetric matrix of size n d × n d , which serves as the parameter in the CC equations.The molecule lithium hydride (LiH), with d=4 electrons in n=8 orbitals, is our running example.The Hamitonians derived in Example 4.3 furnish the input data for Examples 6.4 and 6.5.
Our theme in Section 5 is the CC equations.We define them in ( 26)-( 27) via a rank constraint on V σ , and we give a reformulation in terms of cluster amplitudes in (30)-(31).For most CC variants used in computational chemistry, this agrees with the traditional formulation [15].But our equations differ for some others (Theorem 5.11).Starting from the general bound in Theorem 5.2, we offer a detailed study of the CC degrees of truncation varieties.
In Section 6 we turn to numerical solutions of the CC equations, both for generic Hamiltonians and for systems derived from chemistry.We present computations with the software HomotopyContinuation.jl [4], together with its certification feature [5].The use of monodromy loops is essential.Our findings show that the new theory leads to considerable practical advances.This is documented in Examples 6.2 and 6.3.Examples 6.4 and 6.5 offer case studies for lithium hydride (LiH) where d = 4, n = 8, and for lithium (Li) where d = 3, n = 8.For a comparison with previous work, [10,Section 6] reports that CCSD with three electrons in six spin-orbitals "supersedes the abilities of state-of-the-art algebraic geometry software", and [10, Theorem 4.10] offers the upper bound 2 27 = 134217728 for CCdeg 3,6 ({1, 2}).Yet, the true CC degree is 55, by Proposition 5.8.It is now instantaneous to solve these CC equations.In short, new algebraic geometry leads to progress in practice.

Exponential Parametrization
We work in the vector space H = ∧ d R n with its standard basis vectors e Motivated by nonlinear algebra [22,Chapter 5], we call ψ I the Plücker coordinates.Sometimes it is preferable to write the Plücker coordinates as c α,β , where α is a subset of [d] and β is a subset of [n]\[d] = {d + 1, . . ., n} of the same cardinality |α| = |β|.The c α,β are known as configuration interaction coefficients in quantum chemistry.The identification between these two systems of coordinates on the space of quantum states H = ∧ d R n is as follows: We think of the ψ I as the d × d minors of a d × n matrix, and we think of the c α,β as the minors of all sizes in a d × (n − d) matrix.These two sets have the same cardinality The title of this section refers to a birational map, defined shortly, between two copies of R ( n d ) .It restricts to a polynomial map with polynomial inverse on the affine hyperplane Later on, when we come to algebraic varieties, we shall pass from the vector space H to the projective space P ( n d )−1 = P(H).This is the projective closure of the affine space H ′ .Thus the above coordinates ψ I and c α,β also serve as homogeneous coordinates on P ( n d )−1 .See [8,Chapter 8] for basics on projective algebraic geometry with a view toward computation.
To define the exponential parametrization, we introduce our second vector space V.This is isomorphic to H, with coordinates indexed by [n]  d .The elements of V are called cluster amplitudes and we denote them by x = (x I ) I∈( [n] d ) .The cluster amplitudes x also have alternate coordinates that are indexed by minors of a d × (n − d)-matrix.As in (1), we set The level of a coordinate ψ I or x I is defined as the cardinality of I\ [d].Equivalently, the level of c α,β or t α,β equals |α| = |β|.For example, for d = 3 and n = 6, each of the spaces H and V has 20 coordinates: one of level 0, nine of level 1, nine of level 2, and one of level 3: The term level refers to the excitation level of the electrons in a chemical system.
Our workhorse is the nonlinear coordinate transformation between quantum states and cluster amplitudes.The basic ingredient is a lower-triangular matrix T (x) of square format The entry of T (x) in row J and column I is zero unless I\J ⊆ [d] and (J\I)∩[d] = ∅.If this holds then the matrix entry is ± t I\J, J\I , where the sign is defined as follows.Similarly to I, the sets J, I\J, J\I and I ∩ J are subsets whose elements are written in (increasing) order.To be precise, if , then the sequence (I ∩ J, I\J) is a permutation of I.The sign of t I\J,J \I in T (x) J,I is the sign of the permutation I → (I∩J, I\J) times the sign of the permutation J → (I ∩ J, J\I).
Using x-coordinates, the entry of T (x) in row J and column I equals the above sign times In conclusion, T (x) is a well-defined lower-triangular matrix of size n d × n d that depends linearly on the cluster amplitudes x.This matrix represents the cluster operator in [9,10].
The lower-triangular matrix T (x) is nilpotent of order d.This is shown in [10, Section 3.2] and also follows from equation (13).Hence, the matrix exponential is the finite sum In particular, the entries of the matrix exp(T (x)) are polynomials in x of degree at most d.
Thus the Grassmannian Gr(2, 5) ⊂ P 9 makes an appearance in the first column of exp(T (x)).
Returning to general d and n, the exponential parametrization is the map Here e [d] is the reference state in H ≃ R ( n d ) , i.e. the first standard basis vector.The transformation (5) gives a formula for the quantum states ψ in terms of the cluster amplitudes x.To be precise, each of the n d coordinates ψ I is a polynomial ψ I (x) in the n d unknowns x J .In the definition (5), we had assumed that x [d] = 1 and ψ [d] = 1.Geometrically, this means that we work in the affine spaces V ′ and H ′ , both of which are identified with R ( n d )−1 .Later on, we shall extend (5) to a birational automorphism of the projective space P ( n d )−1 .The reference coordinates x [d] = t ∅,∅ and ψ [d] = c ∅,∅ will then serve as homogenizing variables.
The formula ψ = exp(T (x))e [d] simply says that ψ equals the leftmost column vector of the matrix exp(T (x)).For instance, in Example 2.2, the formula for the Plücker coordinates of the quantum state ψ = (ψ 12 , ψ 13 , . . ., ψ 45 ) in terms of cluster amplitudes is given by the first column of (4).Here is a slightly larger example, where the matrices have size 20 × 20: ).The 20 coordinates in the formula (5) are as follows: In both examples, we can easily solve the equation ψ = exp(T (x)) e [d] for x.This is done inductively by level.For levels zero and one, we simply have x I = ±ψ I .At each larger level, we use the formulas for lower level coordinates x I in terms of the ψ J , and we substitute these into the equation.For instance, in Example 2.3, this yields the inversion formulas We next show that such inversion formulas can always be found.
Proposition 2.4.The map (5) from x-coordinates to ψ-coordinates has a polynomial inverse.Namely, x I is equal to ±ψ I plus a polynomial in ψ-coordinates of strictly lower level.
Proof.We proceed by induction on the level.At level zero, we have For level one, we already saw that x I = ±ψ I .If the index I has level r then the formula in (5) writes ψ I as ±x I plus a polynomial in variables x J of level < r.Each of these lower level x J can now be replaced with a polynomial in ψ, by the induction hypothesis.This yields the promised representation for x I as ±ψ I plus a polynomial in lower level ψ-coordinates.
Let x I (ψ) denote the polynomial that expresses the cluster amplitudes in terms of the Plücker coordinates.Conversely, ψ I (x) is a polynomial in cluster amplitudes.As is apparent from the proof of Proposition 2.4, the degree of each polynomial is the level of I and the variable x I occurs linearly in ψ I (x).Similarly, the variable ψ I occurs linearly in x I (ψ).
The monomials occurring in ψ I (x) are in natural bijection with those occurring in x I (ψ), and we shall give an explicit formula for these monomials and their coefficients.For this, it helps to note that, for each degree d, there is really only one ψ-polynomial and x-polynomial.

Here we refer to the symmetric group actions that permute the indices in [d] and in [n]\[d].
For fixed d and fixed level |I\[d]|, all ψ I (x) are in the same orbit, and ditto for x I (ψ).Each coordinate in the exponential parametrization is a replicate of a certain master polynomial.
It would be desirable to better understand our equations from the perspective of representation theory.In this setting, the master polynomials should be highest weight vectors.
The master polynomials of degree d are ψ I (x) and x I (ψ) where n = 2d and (5) and its inverse are obtained from these two by changing indices.For instance, all quadratic entries in the matrix (4) are replicates of ψ 34 (x) = x 14 x 23 −x 13 x 24 +x 34 .
By inverting the map (5), we find the following master polynomials of degree d ≤ 5:  d )−1 will be important in Section 3. We close this section by giving explicit combinatorial formulas for the master polynomials.Formulas for all other coordinates in (5) and their inverse are obtained by adjusting the indices.The number of monomials is found in The On-Line Encyclopedia of Integer Sequences, which is published electronically at http://oeis.org.Namely, it is the sequence We recall (e.g. from [23]) that a uniform block permutation is a partition which satisfies We denote by U d the set of all uniform block permutations of [2d].For details on the algebraic and combinatorial structures of U d we refer to [23, Section 2.3] and the references therein.
The connection to coupled cluster theory arises from identifying U d with the set of monomials that appear in the above polynomials ψ I (x) and x I (ψ).To see this, we pass to the coordinates in ( 1) and ( 2).The polynomial ψ I (x) is now written as c α,β (t), and x I (ψ) is now written as t α,β (c).With this notation, the master polynomials of degree [d] (c), and we write the monomial corresponding to a given uniform block permutation as The master polynomials are Z-linear combinations of these monomials for k = 1, 2, . . ., d.
Theorem 2.5.The coordinates in the exponential parametrization are In these formulas, the sign of an element π ∈ U d is the product of the signs of the permutations where each α i and β i is an increasing sequence.We also have ν Its term is the product of matrix entries of T (x) whose respective rows and columns are where A r = ∪ r i=1 α i and B r = ∪ r i=1 β i .One checks from the definition of T (x) that the sign of such an entry comes from the number of inversions of (B r−1 , β r ) and (α r , [d]\A r ) and from the sign (−1) |αr| (d−|αr|) .We take the product of the signs of these entries for all r ≤ k.Since is an even integer, the sign of π equals the product of the signs of the permutations in (8).
The set U d of uniform block permutations has a natural partial order, induced by the partial orders on the set partitions of [d] and [d].The Möbius function of the poset U d is given by µ(π) = (−1) k−1 (k − 1)!.For any uniform block partition ρ ∈ U d we can write Using Möbius inversion, we obtain the asserted formula for t [d], [d] in terms of the c π .

Truncation Varieties
In this section, we study the algebraic varieties that are promised in the title of this article.They are found by truncating the exponential parametrization (5) to a certain coordinate subspace.We consider the image of this truncation in H. Its closure in P ( n d )−1 is our variety.More precisely, let σ be a non-empty proper subset of [d] = {1, 2, . . ., d}, and define This is a linear subspace of the vector space V spanned by all basis vectors e J of level in σ.
The subspace V σ is the variety of the ideal The restriction of the exponential parametrization to the subspace V σ is injective.It maps V σ into the full space of quantum states H, and it maps further to the projective space We define the truncation variety V σ as the closure of the image of V σ under this map to P ( n d )−1 .Since the exponential parametrization is invertible, the dimension of the projective variety V σ is one less than the dimension of its linear space of parameters V σ .
The varieties V σ correspond to the various models in CC theory.For instance, in the notation of [10, Section 1.2], the index set σ = {1, 2} corresponds to CCSD, the index set σ = {1, 2, 3} corresponds to CCSDT, etc.But here we allow arbitrary truncation sets.For instance, taking σ = {2, 3} means that doubles and triples are included but singles are not.
Our next result characterizes the homogeneous prime ideals of the truncation varieties.It shows how to derive these ideals from the master polynomials that are given in Theorem 2.5.
Theorem 3.2.The homogeneous prime ideal of the truncation variety V σ is the saturation In particular, the restriction of V σ to the affine chart C ( The provided explicit description of the ideal of the truncation variety allows the computation of deg(V σ ) -and hence the bound in Theorem 5.2 -via the degree of the ideal.For a textbook introduction to the saturation in (11), see Definition 8 in [8, Section 4.4].Its role in the context of projective geometry is explained in [8,Section 8.5].These references and the following example are meant to help our readers in understanding Theorem 3.2.
Proof of Theorem 3.2.We write C[ψ] and C[x] for the rings of polynomial functions on H ′ and V ′ respectively.These are polynomial rings in n d − 1 variables, where The linear subspace V σ in (9) corresponds to the ideal P σ in (10).We write γ σ : C[x] → C[x]/P σ for the associated quotient map.By definition, the prime ideal of V σ ∩ C ( n d )−1 is the kernel of γ σ • ι.This equals ι −1 (P σ ), and it is a prime ideal because P σ is prime.Hence, This is the inhomogeneous prime ideal defining the irreducible affine variety V σ ∩ C ( d n )−1 .We pass to the prime ideal of the projective closure V σ ⊂ P ( n d )−1 by the saturation in (11).
Example 3.4 (d=3, n=6).There are six distinct truncation varieties V σ in P 19 .We compute their prime ideals by the formula in (11).Each item is indexed by the corresponding set σ: {2} The linear space V {2} ≃ P 9 is the zero set of the ten coordinates ψ I of level 1 or 3.

{3}
Here we obtain the line V {3} ≃ P 1 that is spanned by the two points e 123 and e 456 .
{2, 3} The linear space V {2,3} ≃ P 10 is the zero set of the nine coordinates ψ I of level 1.
These computations were carried out with the computer algebra system Macaulay2 [13].
We have already seen Grassmannians a few times for σ = {1}.This is a general result: Theorem 3.5.The truncation variety V {1} equals the Grassmannian Gr(d, n) in its Plücker embedding in P ( n d )−1 .The truncation varieties V σ are thus generalizations of Grassmannians.
Proof.We start with an alternative characterization of the matrix T (x).Recall from [9,10] that the excitation operators are the following graded endomorphisms of the exterior algebra: Here Here X α,β is the matrix for χ α,β restricted to H. Note that T (t) is our earlier matrix T (x).
We now consider the truncation to σ = {1}.Let L σ denote the subspace of L spanned by χ i,j : Operators in L σ are derivations, i.e. for τ ∈ L σ we have Let T k denote the matrix for the linear map τ ∈ L σ restricted to ∧ k R n .From ( 14) we infer We further note that T k has nilpotency k, that is T k+1 k = 0. Since the summands on the right hand side of (15) commute, we conclude that exp (T k+1 ) = exp(T k ) ∧ exp(T 1 ).The analogous property for Kronecker sums appears in [16,Theorem 10.9].By induction on k, Formula (16) shows that the first column of the matrix exp (T (t)) consists of the d × d minors of the first d columns of the n×n matrix exp (T 1 (t)) = Id n +T 1 (t).These columns are This holds because the operator τ ∈ L σ acts on the basis vectors e i of R n = ∧ 1 R n as follows: t i,j e j for i ≤ d and τ e i = 0 for i > d.
We conclude that the first column of exp (T (t)) gives the Plücker coordinates for Gr(d, n).
Remark 3.6.All Grassmannians are obtained from the polynomials xI (ψ) with |I\[d]| ∈ {2, 3, . . ., d}, by the saturation in (11).Starting with these polynomials for other level sets, we obtain all truncation varieties.This is the sense in which the V σ generalize Gr(d, n).
There is a natural isomorphism between the vector spaces ∧ d R n and ∧ n−d R n , and hence between corresponding projective spaces.This swaps the Grassmannians Gr(d, n) and Gr(n− d, n).The duality extends to all truncation varieties.This is the content of the next result, which is our algebraic interpretation of particle-hole symmetry in electronic structure theory.
for the second copy of P ( n d )−1 .The natural isomorphism in the statement of the proposition is given by relabelling as follows: One checks that our construction of the matrix T (x) and the exponential parametrization in Section 2 are invariant under this relabeling.It hence induces the desired isomorphism.
In light of this proposition, we shall assume n ≥ 2d in everything that follows.In Example 3.4 we saw a first census of truncation varieties.We next present two further cases.The last row lists the CC degrees, to be introduced in Section 5.The fourth row gives the number of minimal generators in degrees 1, 2, 3 of the ideal I(V σ ).All varieties live in P 34 .The first column is the Grassmannian Gr(3, 7).Among the other five, three are linear spaces.
Example 3.9 (d = 4, n = 8).The 14 varieties live in P 69 .Five of them are linear spaces: Our methodology for computing these degrees and CC degrees will be explained in Section 6.
We saw in our examples that the truncation variety V σ is a linear space for various subsets σ of [d].The final theorem in this section identifies those subsets for which this happens.
Theorem 3.10.The truncation variety V σ is a linear subspace of P ( n d )−1 if and only if the index set σ is closed under addition, i.e. if i, j ∈ σ with i + j ∈ [d] then i + j ∈ σ.
Proof.We identify V σ with its restriction to the affine chart Here a j ∈ Z\{0} and Next suppose that σ is not closed under addition.Fix the smallest k ∈ [d]\σ such that k = i + j for some i, j ∈ σ.By the same argument as above, ψ I = ±x I (ψ) = 0 for all ψ I of level i < k where i ∈ [d]\σ.Consider any level k equation ( 17) that holds on the variety V σ .Fix any degree-compatible monomial order.The initial monomial of x K (ψ) has degree > 1: where r j ≥ 2.
Since k is minimal, no monomial appearing in x I (ψ) ∈ P σ = I(V σ ) divides the above initial monomial in(x K (ψ)).Hence in(x K (ψ)) is a minimal generator for the initial ideal of P σ .This cannot be an initial ideal for a linear variety, and therefore V σ itself is not linear.

Discretization and Hamiltonians
The electronic structure Hamiltonian is a symmetric matrix H of size n d × n d , describing d electrons relative to a discretization with n spin-orbitals.The primary objective is to compute the quantum states that are eigenvectors of H, i.e. chemists wish to solve the equation In CC theory, this eigenvalue problem is replaced by a polynomial system known as the CC equations.These will be formulated in Section 5.In this section we explain where the Hamiltonian H comes from.The material that follows serves as an introductory guide for algebraists.No background in chemistry is assumed, beyond what is taught in high school.
The starting point of our analysis is the electronic Schrödinger equation From this differential equation, we shall derive a finite-dimensional eigenvalue problem (18).The unknown in (19) is the wave function Ψ(x 1 , x 2 , . . ., x d ).The arguments in this function are pairs x i = (r i , s i ) where i , r i , r i ) are points in R 3 that represent the positions of d electrons, and s i ∈ {±1/2} describes the electronic spin (vide infra).We assume that Ψ is sufficiently differentiable [29].By Pauli's exclusion principle, the wave function Ψ must be antisymmetric in its d arguments, i.e. if x i and x j are switched then Ψ is replaced by −Ψ.The Hamiltonian H in ( 19) is a second order differential operator which describes the behavior of d interacting electrons in the vicinity of d nuc stationary nuclei.The formula is The symbol ∆ r i in the leftmost sum denotes the Laplacian All other summands in (20) act on Ψ by multiplication.They contain constants which we now explain.The atoms and their nuclei are indexed by j = 1, 2, . . ., d nuc .The constant Z j is the jth nuclear charge.This is the atomic number listed in the periodic table, i.e.Z j is a positive integer.The position of the jth nucleus is the point R j ∈ R 3 which is also constant.We here consider only charge-neutral molecules.This means that d = dnuc j=1 Z j and d nuc ≤ d.The following molecule will serve as our running example, both here and in Section 6.
The next step is to construct a finite-dimensional space of functions, along with a suitable basis, which contains an approximate solution to the electronic Schrödinger equation (20).The restriction of H to that finite-dimensional space will be our symmetric matrix H.
There are many ways to select a suitable basis, and hence a discretization of H . See [20] for an overview.We apply the method called linear combination of atomic orbitals (LCAO).This is used widely in quantum chemistry.The LCAO method starts with atomic orbitals.
We select a set χ 1 , χ 2 , . . ., χ k } of atomic orbitals.These are sufficiently smooth functions χ i : R 3 → R which are linearly independent.Atomic orbitals encompass substantial physical principles.Readers can refer to [15, Chapters 5, 6 and 8] for a comprehensive explanation and motivation.Notably, atomic orbital basis sets for different atoms are well-documented and available through online data resources like www.basissetexchange.org.The number k of atomic orbitals is greater than or equal to the number d of electrons, i.e. d ≤ k.The equality d = k can hold.The number k of atomic orbitals determines the total count n of one-particle basis functions used in the discretization.We shall be led shortly to n = 2k.Example 4.2 (d = k = 4).Revisiting the lithium hydride (LiH) molecule, we select k = 4 atomic orbitals, three for lithium and one for hydrogen.A graphical representation of these four atomic orbitals is shown in Figure 1.The pictures are iso-surfaces for χ 1 , . . ., χ 4 .An iso-surface has the form {r ∈ R 3 : χ i (r) = c}, for a constant c that can be positive or negative.Before proceeding with the LCAO approach, we need to account for the electronic spin, a crucial degree of freedom in electronic structure theory that distinguishes physical states.The inclusion of spin doubles the size of our basis: each atomic orbital χ i will be replaced by two functions.Even though the electronic Schrödinger equation (20) does not explicitly include the electronic spin, it remains a significant factor in electronic structure calculations [29].This explains the equation n = 2k we stated before Example 4.2.Note that n ≥ 2d.
The electronic spin can take two possible values: "spin up" (+1/2) or "spin down" (-1/2).To incorporate this, we introduce a spin variable s ∈ {±1/2} and two binary functions The atomic orbitals can be separated into spatial and spin components, namely This factorization simplifies the treatment of electronic spin, making it possible to handle the spatial and spin degrees of freedom independently in calculations.In order to simplify notation, we replace R 3 by X := R 3 × {±1/2}, and we introduce compound coordinates x = (r, s) on X.Using these, we equip the atomic orbital space on X with the inner product where the indices of χ i and χ j in the right integral are understood as i and j modulo k + 1.
In principle, we find a Galerkin basis for H by passing from the ϕ i to d-particle functions.However, the LCAO method introduces an additional set of orthonormal functions known as molecular spin orbitals, which describe the behavior of individual electrons within the molecule.The molecular orbitals resemble the atomic orbitals.They are linear combinations: The functions ξ j and ϕ i span the same n-dimensional vector space.The determination of the expansion coefficients C j,i typically involves employing Hartree-Fock theory, as explained in [15,Chapter 10] or [21, Chapter 2.1].For our specific example, lithium hydride, the molecular orbitals obtained from (spin-restricted) Hartree-Fock theory are depicted in Figure 2. We shall now employ the molecular orbitals to construct d-particle functions.Since we are investigating electrons (i.e.fermions, which must adhere to Pauli's exclusion principle), the d-particle functions must be anti-symmetric [29,Chapter 4.2].Starting from the molecular orbitals, we form a Galerkin basis of n d anti-symmetric d-particle functions as follows: In order to calculate matrix elements, we equip the Galerkin space with the inner product The evaluation of this inner product requires the numerical evaluation of non-trivial integrals.
A common way to circumvent this intricate numerical integration is to use Gaussian-type atomic orbitals [15,Chapter 8].With this, we can proceed to discretize the Hamiltonian H .
The following real numbers are the entries in our symmetric matrix H of format n d × n d : This matrix H is known as the electronic structure Hamiltonian in first quantization.Its entries are parameters in the CC equations which will be introduced in the next section.Using PySCF [28], we evaluate the following integrals for all indices p, q, r, s ∈ [8]: This uses the molecular orbital basis {ξ 1 (x), . . ., ξ 8 (x)} obtained from Hartree-Fock theory.
For computing the matrix elements in (25), we expand the d-particle functions Φ I and Φ J using antisymmetry: The fact that each molecular orbital is a function in one variable allows us to factorize the integral expression.In this factorized expression h p,q and v p,q,r,s will make their appearance.Finally, we use the orthogonality of the molecular orbitals with respect to the inner product defined in (22).This yields the following expression for the entries of the desired matrix: In conclusion, by means of PySCF, we compute a 70×70 matrix H, to be used in Example 6.4.

The CC Equations
We now present the coupled cluster (CC) equations.These approximate the eigenvalue problem (18).The ambient space will be one of the truncation varieties V σ .The equations are determined by a Hamiltonian H as in Section 4. A first systematic study, with a focus on Newton polytopes, was undertaken in [10].In what follows we derive the CC equations from the perspective of algebraic geometry, leading to an alternative description.In Theorem 5.11 we examine to what extent our formulation here is equivalent to the one seen in [9,10].
Let σ be a non-empty proper subset of [d].The set of indices with level in σ is denoted We already saw that the dimension of the truncation variety V σ equals the cardinality | σ|.This translates into a system of quadratic equations on the projective variety V σ .This is the system of CC equations.The linear dependency condition in ( 26) can be expressed by the 2 × 2 minors of the (| σ|+1) × 2 matrix (Hψ) σ , ψ σ .This represents the truncation of ( 18): The number | σ| of constraints imposed by this equation equals the dimension of V σ , so we expect the CC equations to have a finite number of solutions for generic H.This is indeed the case.We call this number the CC degree of the variety V σ , denoted by CCdeg d,n (σ).
The CC degree is the number of solutions to the CC equations.Some non-trivial values of CCdeg d,n (σ) were shown in the tables of Examples 3.8 and 3.9.We have the following general upper bound for the CC degree in terms of the degree of the truncation variety V σ .Theorem 5.2.For any Hamiltonian H, the number of isolated solutions to (26) satisfies Proof.We set N = | σ| = dim(V σ ).For generic Hamiltonians H, the variety in P ( n d )−1 defined by the 2 × 2 minors of the (N +1) × 2 matrix (Hψ) σ , ψ σ has codimension N and degree N +1.Geometrically, this variety is a cone over a generic section of the Segre variety P 1 × P N ⊂ P 2N +1 .The first factor on the right-hand side in (28) is N + 1 = deg(P 1 × P N ).The intersection with V σ has expected dimension zero, but it can have higher-dimensional components.We are interested in the number of isolated solutions.By Bézout's Theorem, this number is at most the product of the degrees of the two varieties, seen on the right in (28).That bound on the number of isolated solutions also holds for special matrices H.
We note that the equality holds in (28) for the linear cases described in Theorem 3.10.
We expect (28) to be strict whenever deg(V σ ) ≥ 2. This holds whenever CCdeg d,n (σ) is known.There is a geometric explanation: the intersection of V σ with the variety of 2 × 2 minors of (Hψ) σ , ψ σ is not transverse on the hyperplane V (ψ [d] ).Finding the true CC degree is a problem in intersection theory, just like computing the degree of V σ itself.The two numbers are important for quantum chemistry because they govern the complexity of solving the CC equations.In particular, CCdeg d,n (σ) is the number of paths to be tracked when solving numerically with HomotopyContinuation.jl.This is discussed in Section 6.
Theorem 5.5.The CC degree associated with the Grassmannian Gr(2, n) equals This result was stated as a conjecture in the first version of this article from August 2023.It was proved in October 2024 in collaboration with Viktoriia Borovik.It is published in [3].
The inequality ( 28) is strict for Gr(2, n) because the CC equations admit extraneous solutions that lie on the hyperplane {ψ 12 = 0}.These form a higher-dimensional component which can be removed by saturation as in (11).We discuss this in detail for a small instance.
Indeed, the Grassmannian V {1} = Gr(2, 5) is cut out in P 9 by the 4 × 4 Pfaffians in a skewsymmetric 5 × 5 matrix (see [22,Example 4.9]).Let I be the ideal generated by these 5 Pfaffians plus the 7 2 = 21 maximal minors of the 7 × 2 matrix on the right.This gives the upper bound 35 = 5 × 7 in Theorem 5.2.The ideal I is radical, and it is the intersection of the desired ideal of codimension 9 and a linear ideal of codimension 7, namely ⟨ψ 12 , ψ 13 , ψ 14 , ψ 15 , ψ 23 , ψ 24 , ψ 25 ⟩.This extraneous component is responsible for the differerence 8 between the upper bound and the true CC degree, which is 27, as in Theorem 5.5.In three cases, the variety V σ is a linear space in P 19 , and the CC degree is dim(V σ ) + 1.In the other three cases, the CC degree is a bit below the bound given by Theorem 5.2.See Examples 3.8 and 3.9 for more comparisons between our upper bound and the CC degree.
Of special interest is the case when the truncation variety V σ is a hypersurface.This is the hypersurface defined by the master polynomial t Proof.We wish to count all scalars λ ∈ C such that the truncated eigenvalue equation (H − λ Id)ψ σ = 0 has a solution ψ in V σ .For this, we delete the last row of the matrix H − λ Id to get a matrix with one more column than rows.Using Cramer's Rule, we write the entries of ψ as signed maximal minors of that matrix.The last entry of ψ is a polynomial in λ of degree 2d d − 1, which is the size of these minors.All other entries of ψ = ψ(λ) are polynomials of degree 2d d − 2, because λ does not occur in the last column of our matrix.We substitute the vector ψ = ψ(λ) into the equation f We know that f has degree d and it is linear in the last variable Since H is generic, the polynomial is square free, and its number of complex zeros is given in (29).
We next express the CC equations in terms of the cluster amplitudes x I .Let z denote the restriction of the vector x to the coordinates in σ.To be precise, z is the vector of length n d which is obtained from x by setting x [d] = 1 and x J = 0 for all J ̸ ∈ σ.We identify C | σ| with the space of such vectors z.The truncation variety V σ has the parametric representation .
We substitute this parametrization into the (| σ|+1) × 2 matrix (Hψ) σ , ψ σ seen in (26).Therefore the CC equations ( 27) are equivalent to the following equations in the unknowns z: Here we require column vectors of length | σ| + 1 to be linearly dependent.For computations, it is advantageous to rewrite (30) as a square system of | σ| + 1 equations in | σ| + 1 unknowns: We now compare the system (31) with the traditional formulation of the CC equations, which was used in [10]  It turns out that -sometimes -the traditional formulation (34) yields a polynomial system that is fundamentally different from the system we derived in (31).The reason for this discrepancy is that the left multiplication with (32) need not commute with truncation.
The following result says that the discrepancy we discovered is actually not so bad.It characterizes all CC variants where the old and new formulation of the CC equations agree.
Proof.We consider a matrix A whose rows and columns are indexed by [n]  d .The identity Theorem 5.11 shows that the traditional formulation (34) coincides with our formulation (26) in all cases that have appeared in the computational chemistry literature [9,10,17] including CCS, CCD, CCSD, CCSDT.In particular, for electronic structure Hamiltonians, our formulation (31) contains only expressions where CC amplitudes appear at most to the fourth power, because of the special structure of these Hamiltonians.Scenarios where the two formulations differ, like σ = {2, 3}, are less relevant for coupled cluster theory.For us, ( 26) is more elegant than (34), and its algebraic degree is lower.This is why we refer to (26) as the CC equations.

Numerical Solutions
This section covers the state-of-the-art for computing all solutions to the CC equations.We use the formulation (31) as a square system with | σ| + 1 unknowns.Readers familiar with earlier formulations of the CC equations may consult Theorem 5.11 for the precise relation.
Our new approach allows for the solution of systems much larger than those in [10].This is accomplished by leveraging monodromy techniques.Throughout this section, we use Julia version 1.9.1 and HomotopyContinuation.jl version 2.9.2 [4].The python part is performed using Python 3.8.8 and PySCF 2.0.1 [28].Computations were done on the MPI-MiS computer server using four 18-Core Intel Xeon E7-8867 v4 at 2.4 GHz (3072 GB RAM).
The beginning of our experiments, for given d, n, is the choice of a symmetric matrix H.For a specific truncation level σ, we pick a pair (λ, z) ∈ C × C | σ| at random.We then construct a generic complex matrix H for which (λ, z) is a solution to (31).This is done easily by solving a linear system of equations.The reason is that (31) depends linearly on H.
Having the starting data H and (λ, z) we use the monodromy solver to find all solutions for H.This system reveals the CC degree, and is later used to compute quantum chemistry systems.The degree of V σ is found in a similar manner, by slicing V σ with an appropriate generic linear space.Whenever feasible, we use Macaulay2 [13] to validate the degrees and CC degrees we found numerically.In this manner we found the table entries in Examples 3.8, 3.9 and 5.7.Here is one more case:   low numbers on the diagonal.This is because CC theory is exact in its untruncated limit, i.e. the number of solutions is the matrix size.For CC at level [d − 1], the entries come from Proposition 5.8.The middle entries increase rapidly.The left column [1] concerns the Grassmannian Gr(d, 2d).
Even for larger cases, Theorem 5.2 gives a good upper bound on the number of solutions.The key ingredient is the degree of the truncation variety V σ .This is often easier to compute than CCdeg d,n (σ).For computing deg(V σ ), we used a range of techniques.First of all, the degree is one in the linear cases of Theorem 3.10.Second, for the CCS truncation (σ = {1}), the degree of the Grassmannian has an explicit description (cf.[22,Theorem 5.13]).Third, sometimes we can compute deg(V σ ) symbolically with Macaulay2 [13]; this requires either an explicit description of the ideal of the truncation variety, e.g., as provided in Theorem 3.2, or can sometimes be done by implicitization [22,Section 4.2] from the parametrization (5).The degree of the considered variety can be computed by the degree of the ideal.Finally, if this all fails, we use numerics, taking advantage of the fact that V σ is a complete intersection in affine space C ( n d )−1 , cut out by the polynomials x I (ψ) in (12).Namely, we intersect V σ with a generic affine-linear space of codimension n d − | σ| = codim(V σ ).The number of points in the intersection is the degree of V σ .Using an appropriate parametrization of the affine linear-space we may use the monodromy solver in order to compute the degree.
In conclusion, the inequality in Theorem 5.2 leads to upper bounds for the number of roots of the CC equations, even when the equations are too large to be solved completely.It is instructive to compare previously known bounds to those found by our new approach.Example 6.3 (Scaling of the number of roots).For d = 2 we consider the CC equations for singles (σ = {1}) and doubles (σ = {2}) investigating the scaling of the number of roots with respect to n. Figure 4 shows different bounds in a log-lin plot.The blue curve is the previous bound from [10,Theorem 4.10] and the green curve is our new bound from Theorem 5.2.We moreover show the exact number of roots; for CCD (right panel) this is given in Corollary 5.3 and for CCS (left panel), this is given in Theorem 5.5, here confirmed for n ≤ 10.The graphs show that the algebraic geometry in this paper leads to much improved bounds.Turning to solutions of the CC equations for real data, we use the parameter homotopy method to track paths from the generic starting solutions to solutions for a given quantum chemical Hamiltonian.The number of these paths is the CC degree.If this tracking converges, we find either a non-singular solution or a singular solution.A solution is singular if the Jacobian matrix of the polynomial system is singular (i.e.non-invertible), or if the winding number of the solution path is greater than one.A singular solution could indicate that the solution variety has an extraneous component of dimension ≥ 1.We observe that this is common in applications from quantum chemistry, such as those in Examples 6.4 and 6.5.However, inspecting the eigenvalues only, it appears that singular solutions still yield a good approximation see Example 6.4.For a general Hamiltonian H, all paths converge to non-singular solutions, and the number of solutions to (31) is exactly the CC degree.Note that the Hamiltonian H arising in quantum chemistry (cf.Section 4) is not generic, but has special structure.Therefore, the obtained number of solutions for the target system can be much smaller than the CC degree.The computation of the generic start system for σ = {1} takes 82 minutes and yields CCdeg 4,8 ({1}) = 154441, as predicted by Example 3.9.Tracking all paths yields 3 nonsingular solutions all of which are real.We also find 104641 singular solutions.Only 399 of them yield real energies.We use these for the comparison to the exact eigenvalues.These calculations take 11 minutes and 32 seconds.For σ = {2}, the computation of the generic start system takes 13 seconds and yields CCdeg 4,8 ({2}) = 73.Tracking all paths yields 36 non-singular solutions of which 24 are real and zero singular solutions.This takes less than one second.In Figure 5 we compare the exact eigenvalue spectrum with the energies obtained from CCS and CCD.An interesting observation is that CCS and CCD appear to approximate different subsets of eigenvalues that together cover the whole spectrum.
Since the CC degree for d = 4, n = 8 and σ = {1, 2} is currently unknown, we cannot apply our method to CCSD for lithium hydride yet.In order to investigate the approximation quality for CCSD, we lower the number of electrons by looking at the lithium atom.Example 6.5 (Lithium (d = 3, n = 8, σ = {1, 2})).We find one non-singular solution, which is real, and 2931 singular solutions.Compare this to CCdeg 3,8 ({1, 2}) = 145608.The runtime for the parameter homotopy was about one hour.We compare the eigenvalue spectrum of the 56 × 56 matrix H with the energies obtained from the CCSD computations in Figure 6.
For special Hamiltonians H, the number of isolated solutions to the CC equations can be much lower than the CC degree.This happens in applications, as seen in Examples 6.4  The CC degree is 9 for the Plücker quadric Gr (2,4).Let H be a general symmetric 6 × 6 matrix of rank r.For r = 1, 2, 3, 4, 5, the number of solutions to (31) is 1, 3, 5, 7, 9. We see this with Cramer's Rule as in proof of Proposition 5.8.
In conclusion, this article introduced a new formulation of the coupled cluster (CC) equations in electronic structure theory.This rests on a novel class of algebraic varieties, called truncation varieties.They live in the same projective space as the Grassmannian, which they generalize.Section 6 has demonstrated that our new approach leads to significant advances in numerically computing all roots of the CC equations.Current off-the-shelf software can now reliably solve instances of actual interest in quantum chemistry.These practical advances rest on the theorems in Sections 2, 3 and 5.We believe that those are of interest in their own right.Many open problems and new avenues of inquiry were presented.One concrete task is to find a formula for CCdeg d,n ({1}), which is the CC degree of the Grassmannian Gr(d, n).We hope that this topic will interest experts in intersection theory.

d
is a subset of [n] of cardinality d whose elements are always written in (increasing) order.Here d ≤ n are positive integers.The reference state is the first basis vector e [d] for [d] = {1, 2, . . ., d}.Vectors in H are called quantum states and they are written uniquely as linear combinations of the basis vectors:

Proposition 3 . 7 .
Fix a subset σ of [d] and let n ≥ 2d.There is a linear isomorphism between two copies of P ( n d )−1 which switches the truncation varieties V σ for (d, n) and for (n − d, n).Proof.The Plücker coordinates on the two spaces are ψ I with |I| = d and ψ I ′ with |I ′ | = n−d.Similarly, the coordinates in (1) are c α,β with α ⊂ [d] and β ⊂ [n]\[d] for the first copy of P ( n d )−1 , and

Example 3 . 8
(d = 3, n = 7).The six varieties correspond to the six columns in this table: d]| = k for each j, i.e. the levels of the variables in each monomial sum up to k.Since k ∈ [d]\σ and σ is closed under addition, each monomial contains some ψ I of level i < k and i ∈ [d]\σ.By induction, ψ I = 0.This now implies ψ K = 0.The restriction V σ ∩ H ′ is thus linear and therefore its projective closure V σ as well.

Example 4 . 1 (
Lithium hydride).This molecule has the formula LiH.It involves d nuc = 2 atoms, namely lithium Li and hydrogen H. Their atomic numbers are Z 1 = 3 and Z 2 = 1, so the number of electrons is d = Z 1 + Z 2 = 4.The two nuclei are fixed at locations R 1 and R 2 , whereas the four electrons have variables x 1 , x 2 , x 3 , x 4 .Here, Ψ is an unknown function of 16 scalars.It satisfies

Figure 2 :
Figure 2: Iso-surfaces of four molecular orbitals of lithium hydrate for s = +1/2.The isosurfaces correspond to the value ±0.04, with blue for +0.04 and red for −0.04.

Example 4 . 3 (
Lithium hydride).In our running example, we have d = 4 and n = 2d = 8.Our Hamiltonian H for LiH is a symmetric 70 × 70 matrix which is computed as follows.
holds for a general matrix B if and only if [ A ] σ,σ c = 0, where σ c = [d]\σ.Suppose that [ exp (T (z)) ] σ,σ c = 0. (35)Then the same block is zero in the inverse matrix exp (−T (z)), and therefore[ (H − λ Id) exp(T (z))e [d] ] σ = [ exp (−T (z)) ] σ,σ [ (H − λ Id) exp(T (z))e [d] ] σ = [ exp (−T (z))H exp(T (z))e [d] − λe [d] ] σ .This means that (30)-(31) is equivalent to (33)-(34).Conversely, if this holds then (35) must hold because the Hamiltonian H can be an arbitrary symmetric n d × n d matrix.We shall now prove that (35) holds if and only if σ has the form m[k].The (I, J) entry in exp (T (z)) is non-zero if and only if we can map from state e J to state e I via a composition of linear maps X α,β , where |α| = |β| ∈ σ.This is possible if and only if J\[d] ⊂ I\[d] and| I\[d] | − | J\[d] | = k∈σ p k k for some p k ∈ N.Given any 1 ≤ j < i ≤ d, we can always find I, J ⊆ [n] of levels i and j respectively such that J\[d] ⊂ I\[d].Hence the block [ exp (T (z)) ] σ,σ c is non-zero if and and only if∃ j ∈ [d]\σ ∃ i ∈ σ : i > j and j = i − k∈σ p k k for some p k ∈ N.(36)Let m and M be the minimal and maximal elements of σ.If M − p m m ∈ [d]\σ for some p m ∈ N + then (36) holds and [ exp (T (z)) ] σ,σ c is non-zero.Suppose now that M − p m m ∈ σ for all positive integers p m < M/m.Then M = km by the choice of m, and we have m[k] ⊆ σ.Next suppose m[k] = σ.No j ∈ [d]\σ with j < M is a multiple of m, and thus (36) fails.Finally suppose m[k] ⊊ σ.We take the smallest element i ∈ σ that is not a multiple of m.Then i − m ∈ [d]\σ and (36) holds for p m = 1.Therefore [ exp (T (z)) ] σ,σ c is non-zero.

Figure 4 :
Figure 4: Bounds to the number of roots of CCS (left panel) and CCD (right panel).

Figure 5 :
Figure 5: The real eigenvalues from exact diagonalization (FCI), shown as red bars, are compared to the energy spectra, shown in black, from CCS and CCD (in the STO-6G basis).

Figure 6 :
Figure 6: Energy spectra from exact diagonalization (FCI) and from CCSD (in the STO-6G basis).The left panel shows the full spectrum.The right panel shows the spectral values between −13.5 and −8.5.and 6.5.It would be interesting to understand this phenomenon for H on special loci in the space of symmetric matrices.The next example suggests such a study for low rank matrices.
3,4} ≃ P 53 .The other nine are listed here: ).The CC systems for the six varieties for three electrons in The number of real solutions (listed in row "#real") varies for different choices of real-valued Hamiltonians, unless V σ is a linear space.The counts 430, 1376, 658 (in row "#real") are from one representative sample run for a random real-valued H.The degree of V {1,2} is computed numerically.The runtimes in seconds are for solving and certifying for generic H. Example 6.2.The CC degrees found for various d, n and σ indicate the complexity of fully solving the CC equations.