Diagrammatically size-consistent and basis-set-free vibrational coupled-cluster (XVCC) theory for both zero-point energies and transition frequencies of a molecule, the latter through the equation-of-motion (EOM) formalism, is defined for an *n*th-order Taylor-series potential energy surface (PES). Quantum-field-theoretical tools (the rules of normal-ordered second quantization and Feynman–Goldstone diagrams) for deriving their working equations are established. The equations of XVCC and EOM-XVCC including up to the *m*th-order excitation operators are derived and implemented with the aid of computer algebra in the range of 1 ≤ *m* ≤ 8. Algorithm optimizations known as strength reduction, intermediate reuse, and factorization are carried out before code generation, reducing the cost scaling of the *m*th-order XVCC and EOM-XVCC in an *n*th-order Taylor-series PES (*m* ≥ *n*) to the optimal value of *O*(*N*^{m+⌊n/2⌋}), where *N* is the number of modes. The calculated zero-point energies and frequencies of fundamentals, overtones, and combinations as well as Fermi-resonant modes display rapid and nearly monotonic convergence with *m* towards the exact values for the PES. The theory with the same excitation rank as the truncation order of the Taylor-series PES (*m* = *n*) seems to strike the best cost-accuracy balance, achieving the accuracy of a few tenths of cm^{−1} for transitions involving (*m* − 3) modes and of a few cm^{−1} for those involving (*m* − 2) modes. The relationships between XVCC and the vibrational coupled-cluster theories of Prasad and coworkers and of Christiansen and coworkers as well as the size-extensive vibrational self-consistent-field and many-body perturbation theories are also elucidated.

## I. INTRODUCTION

Coupled-cluster (CC) theory^{1,2} was introduced by Coester and Kümmel^{3,4} in nuclear physics and soon adopted by Čížek and coworkers^{5,6} in the field of quantum chemistry, where it was fully developed into the most accurate, versatile, and widely applicable electronic structure theory.^{7–9} Being rigorously (i.e., diagrammatically) size consistent,^{10,11} it was applied to problems in solid-state physics.^{12–18} It also witnessed a revival in nuclear physics,^{19,20} as the Hamiltonian for nucleons became increasingly more accurate. Another area where the applicability of this theory is being explored is quantum nuclear dynamics in a bound potential energy surface (PES), i.e., anharmonic vibrational structure theory.^{21,22}

Following its applications to model PES’s,^{23–25} vibrational CC (VCC) theory was pioneered and extensively developed by Prasad and coworkers for general molecules with truncated Taylor-series PES’s.^{26–34} They considered both zero-point energies and transition frequencies using the equation-of-motion (EOM)^{27,32} or multireference^{29,34} formalism. They also calculated spectra and band intensities in either the time-dependent^{26,28,30,31} or independent^{33} framework. These methods are based on a single Hartree product of harmonic-oscillator wave functions as a reference, of which the mode frequencies and coordinate center are varied to minimize its energy. Their working equations are derived using the second-quantized ladder operators, which are valid for harmonic oscillators. The reference wave function adopted can be identified as that of the size-extensive vibrational self-consistent-field (XVSCF) method with an anharmonic geometry correction, i.e., XVSCF[*n*] (*n* is the rank of the Taylor-series PES).^{35}

Christiansen and coworkers^{36–42} proposed alternative formulations of VCC theory and its excited-state analog, VCC linear response (VCCLR). They may be viewed as a direct extension of vibrational correlation theories pioneered by Bowman and coworkers and by Ratner and coworkers, such as vibrational configuration-interaction (VCI)^{43} and vibrational Møller–Plesset perturbation (VMP) theories,^{44} and share similar design philosophy. They all rely on a single-reference Hartree product of one-mode functions (“modals”) along normal coordinates furnished by vibrational self-consistent-field (VSCF) theory.^{45–47} These modals are determined so as to minimize the VSCF energy and are numerically defined on a grid or a basis set. These vibrational correlation theories tend to be implemented in symbolic (as opposed to matrix-algebraic) algorithms, in which excited Hartree products are manipulated during a calculation (not during formula derivation) with the Hamiltonian matrix elements in modals evaluated by quadrature.

In the meantime, it has been shown that VSCF theory is not diagrammatically size consistent.^{48} Unlike a truncated CI for electrons, which gives zero electron-correlation energy per electron in the thermodynamic limit, VSCF yields a finite and thus physically meaningful anharmonic correction in this limit. Nevertheless, this correction is due to a minute portion of a large number of calculated contributions, most of which are shown to vanish in the same limit as they correspond to disconnected diagrams. This renders VSCF an exceedingly inefficient method for large molecules and solids. When such vanishing contributions are erased from the VSCF equation, we arrive at a streamlined, rigorously size-consistent, but less accurate variant of VSCF, which we call XVSCF(*n*) or XVSCF[*n*],^{35,49,50} depending on whether an anharmonic geometry correction is suppressed (in the former) or not (in the latter). XVSCF is no longer exact for any one-dimensional potential, which VSCF is, in the infinite basis-set limit. Instead, each modal of XVSCF (or that of VSCF in the thermodynamic limit) is shown to be a harmonic-oscillator wave function with its frequency and center dressed with mean-field anharmonic effects,^{48,49} which is consistent with Makri’s theorem.^{51} Since the Hamiltonian matrix elements in harmonic-oscillator modals are expressed simply with force constants, XVSCF and vibrational correlation theories based on it can be defined compactly with force constants and excitation amplitudes. Grids, basis sets, and quadrature are no longer necessary.

Insofar as the defining feature of CC theory is its size consistency, which also ensures its uniformly high accuracy across system sizes, its formalism should be stipulated with many-body-theoretical techniques such as normal-ordered second quantization and Feynman–Goldstone diagrams and be based on a rigorously size-consistent reference wave function. A normal-ordered second-quantized form of the vibrational Hamiltonian has been introduced recently.^{52} Remarkably, its constant term and excitation amplitudes are identified as the energy, gradient, and frequency expressions of XVSCF and of the closely related self-consistent-phonon (SCP) method^{53–57} (but not of VSCF). This is reminiscent of the normal-ordered electronic Hamiltonian, in which the expressions of the Hartree–Fock (HF) energy and Fock operator naturally emerge.^{2} This observation highlights the central place XVSCF and SCP occupy in anharmonic vibrational structure theory for large molecules and solids and suggests that they provide an ideal reference wave function for vibrational correlation theories for such systems. The establishment of XVSCF,^{35,49,50} along with the elucidation of their relationship to the normal-ordered vibrational Hamiltonian,^{52} therefore, sets the stage for a more complete development of a VCC theory.

In this article, we present a rigorously (i.e., diagrammatically) size-consistent VCC theory for the ground and excited states, the latter through the EOM formalism. We call them XVCC*m* and EOM-XVCC*m* theories, respectively, where *m* stands for the highest connected excitation rank considered. We propose a rational and expedient method of deriving their algebraic definitions using the rules of normal-ordered second quantization and Feynman–Goldstone diagrams. We demonstrate its utility through the derivation of low-order XVCC and EOM-XVCC equations. The formalisms thus derived are valid for any single-reference Hartree product of harmonic-oscillator modals along delocalized, rectilinear coordinates, such as XVSCF(*n*),^{49} XVSCF[*n*],^{35} and SCP (Refs. 53–57) as well as the harmonic approximation. Some diagrammatic contributions are shown to vanish in the XVSCF and SCP reference, making them more favorable, although higher-order XVCC and EOM-XVCC are increasingly insensitive to the reference choice. XVCC and EOM-XVCC with the XVSCF[*n*] reference are found to be equivalent to the VCC formalism first developed by Prasad and coworkers,^{27,32} shedding light on the relationship between the two VCC theories and XVCC theory: the VCC theory of Prasad and coworkers and XVCC theory may be diagrammatically size-consistent and infinite-basis counterparts of the VCC theory of Christiansen and coworkers in a truncated Taylor-series PES, if no other approximations are made.

We furthermore derive and implement high-order XVCC*m* and EOM-XVCC*m* equations up to *m* = 8 or octuple excitations, with the aid of computerized derivation, optimization, and code synthesis. We apply these methods to the zero-point energies and frequencies of the fundamentals, overtones, and combinations as well as Fermi-resonant modes of the water and formaldehyde molecules, using quartic force fields (QFF’s).^{58} A comparative analysis is performed for the results of these methods and of XVSCF(*n*), XVSCF[*n*], size-extensive second-order vibrational many-body perturbation theory (XVMP2),^{59} and full VCI. It shows that XVCC*m* and EOM-XVCC*m* with *m* equal to the truncation rank of the Taylor-series PES are uniformly accurate and stable for both resonant and nonresonant modes and are likely the methods with the greatest practical utility of all considered.

## II. THEORY

### A. XVCC formalism

We express the *m*th-order XVCC (XVCC*m*) wave function for the zero-point (ground) state (labeled by subscript “0”) as

where |Φ_{0}〉 is a reference wave function. The latter, in turn, is a Hartree product of harmonic-oscillator modals (which include but are not limited to those from the harmonic approximation to the PES),

where |*n _{i}*〉 denotes the harmonic-oscillator wave function with quantum number

*n*and frequency

_{i}*ω*along the

_{i}*i*th normal, first-order Dyson,

^{50}or other delocalized, rectilinear coordinate. Operator $T\u02c6$ is the sum of the single, double, through

*m*-fold cluster excitation operators,

with

where *τ*_{i1i2⋯in} is an excitation amplitude, the value of which is to be determined by solving the amplitude equations (see below). Operator $a\u02c6i\u2020$ and its adjoint $a\u02c6i$ are the ladder operators for the harmonic-oscillator modals, i.e.,

The braces {⋯} in Eq. (5) bring the operators enclosed by them into a normal order.^{52} Two or more of the mode indices (*i*_{1}, *i*_{2}, …, *i _{n}*) in Eq. (4) are permitted to be coincident.

We require that the wave function of Eq. (1) satisfy the vibrational Schrödinger equation in the space of Hartree products reachable by $T\u02c6$ from |Φ_{0}〉,

and

for all *ν*_{1}, …, *ν _{n}* and each

*n*(1 ≤

*n*≤

*m*). Here, $E0(m)$ is the XVCC

*m*total energy and

The number of the amplitude equations expressed collectively by Eq. (9) is equal to the number of unknown excitation amplitudes {*τ*} in Eq. (4). Furthermore, it can be shown below that $E0(m)$, which is also unknown, is eliminated from Eq. (9). Hence, by solving Eq. (9), we obtain {*τ*}. Substituting the latter in the energy equation [Eq. (8)], we determine $E0(m)$.

The reference wave function |Φ_{0}〉 can be furnished by any of the following: XVSCF(*n*),^{49} XVSCF[*n*],^{35} SCP,^{53–57} or the harmonic approximation. They differ from one another in the coordinates, their center, and the modal frequencies,^{50} but they all share the same form: the Hartree product of harmonic-oscillator modals in delocalized, rectilinear coordinates. The XVCC*m* formalisms are essentially unchanged by the choice of the references, except that some terms may vanish in some references. The foregoing definition of XVCC*m* and the working equations derived in the following are valid for any of them.

The rather restricted form of reference wave functions is dictated by the fact that every modal of a Hartree product in delocalized, rectilinear coordinates (such as normal coordinates) reduces to a harmonic oscillator in the thermodynamic limit, as shown both analytically^{48} and numerically^{49} in the case of VSCF and also in a more general sense by Makri.^{51} This makes it convenient to express the formalisms of vibrational correlation methods in terms of the ladder operators and force constants (as opposed to molecular integrals) and also renders the methods free of numerically defined basis functions, as discussed below. Its drawback is a lack of choices of coordinates^{60–62} and the poor performance of the reference methods for strongly anharmonic potentials such as a Morse potential and a double-well potential, for overtones and combinations, and for anharmonic resonances. As shown below, some of these issues are ameliorated by XVCC, but others (such as a double-well problem) are unlikely to be resolved.

One may notice at this stage a striking difference between electronic CC and XVCC. Since the one-mode basis functions of XVCC are harmonic-oscillator wave functions, which need not be stored numerically, XVCC can take account of an infinite number of them implicitly. This means that, even when the Hamiltonian is truncated after the *n*-body term (*n* ≥ 3) for a molecule with *N* vibrational degrees of freedom, neither XVCC*n* nor XVCC*N* is exact mathematically, though XVCC*m* is shown below to be rapidly convergent with *m* in practice. In other words, XVCC is an infinite-basis or basis-set-free method despite the finiteness of its equations. This is where XVCC and Christiansen’s VCC (Ref. 36) also differ from each other fundamentally. The excitation space reached by $T\u02c6$ is limited by the number of basis functions included in VCC, whereas it is not in XVCC.

The VCC ansatz of Prasad and coworkers^{27,32} is the same as above, but it is predicated on a reference Hartree product of Gaussian modals along normal coordinates, whose centers and widths are varied to minimize its energy. This reference method is identified as the XVSCF[*n*] method.^{35} Therefore, XVCC*m* with the XVSCF[*n*] reference corresponds to the VCC theory of Prasad and coworkers,^{27,32} although XVCC*m* is defined for more general reference Hartree products, including XVSCF(*n*), SCP, and the harmonic approximation.

### B. Algebraic derivation of XVCC

Here, we illustrate an algebraic derivation of the working equations of XVCC*m* using the normal-ordered second quantization and Wick’s theorem.^{52} We use, as an example, XVCC1 or XVCC with single excitations (XVCCS) for a vibrational Hamiltonian with a QFF.

This Hamiltonian can be written in the normal-ordered second-quantized form^{52} as

with

where it should be understood that there is term-by-term correspondence between the left- and right-hand sides of each equation.

The constant term in Eq. (11), $E0(0)$, is the total ground-state energy expression of XVSCF,^{35,49,52} regardless of the actual reference method adopted in the calculation, and is given by

where *V*_{0} is the value of the PES at the minimum, *ω _{i}* is again the frequency of the

*i*th modal, and

*F*’s are the force constants. The Hamiltonian amplitudes, {

*W*}, are given by

in a QFF. These amplitudes are invariant to any permutation of indices. With the XVSCF(4) reference, $Wii+=\omega i$ holds. In the XVSCF[4] reference, additionally, *W _{i}* = 0 is true. Furthermore, $Wij+=0$ for

*i*≠

*j*with the SCP reference. The fact that these key quantities of the XVSCF methods naturally emerge in the expression of the normal-ordered Hamiltonian underscores their fundamental significance in many-body vibrational theory.

^{52}

From Eq. (8), we find the XVCCS energy equation to be

where $E0(1)$ is the XVCCS total energy. The Taylor series of the exponential operator is truncated after the $T\u02c614$ term because the higher-order terms cannot contribute in a QFF (see below). The XVCCS amplitude equation, Eq. (9), similarly reduces to

in a QFF. Here, the Taylor series can be terminated after the $T\u02c615$ term (see below).

Wick’s theorem states that the vacuum expectation value of a product of two or more normal-ordered series of operators is the sum of all full contractions with no internal contractions.^{2,52,63} It also states that the vacuum expectation value of just one normal-ordered series is zero.^{2,52,63}

Using this theorem, one can evaluate the first term of Eq. (21) as

The second term has only one nonzero full contraction,

The third term has two nonzero full contractions,

The fourth and fifth terms are evaluated similarly, yielding

and

Since the Hamiltonian with a QFF has only up to four-body operators, the $T\u02c615$ term and higher in the exponential operator cannot produce a full contraction.

Together, the XVCCS energy is given by

It may be noticed that each of the second and subsequent terms is connected in the sense that every factor has at least one common summation index with another factor. This ensures the extensivity of $E0(1)$.^{10,11}

The connectedness of the energy equation is a universal feature of XVCC*m* for all *m*. For future convenience, therefore, we introduce an XVCC*m* energy equation alternative to Eq. (8), which reads

where $H\u02c6N=H\u02c6\u2212E0(0)$, $\Delta E0(m)=E0(m)\u2212E0(0)$, and (⋯)_{C} deletes terms that are disconnected. We call this the energy equation of XVCC*m*.

The XVCCS amplitude equation can likewise be derived by Wick’s theorem. The first term of Eq. (22) has one nonzero full contraction,

The second term has two nonzero full contractions,

The third term has four such contractions,

The remainders are evaluated similarly, yielding

The higher-order terms in the exponential operator cannot produce a full contraction and are zero.

Combining these, we find the XVCCS amplitude equation to read

Note that the first five terms in the right-hand side are disconnected in the sense that there is a factor (*τ*_{ν}) that has no common summation index with any other factor; each of them is, therefore, a simple product of two factors.

Factoring these disconnected terms, we can rewrite the above equation as

where we use the fact that the sum in the parenthesis is equal to $E0(1)$ as per Eq. (28). We can, therefore, further simplify the amplitude equation into

Notice that this equation consists of connected terms only.

In fact, we can recast the amplitude equations of XVCC*m* in general into a form that is equivalent to, but more convenient than Eq. (9),

for all *ν*_{1}, …, *ν _{n}* and all

*n*(1 ≤

*n*≤

*m*). We call this the

*T*amplitude equation of XVCC

_{n}*m*. The diagrammatic derivations of XVCC and the EOM-XVCC formulations to be discussed below are based on Eqs. (29) and (40).

In the foregoing derivation, the connected structure of the XVCC energy and amplitude equations is only illustrated with the XVCCS example. In Appendix A, we document an alternative derivation that proves more rigorously the connectedness of these equations at any order.^{1,2}

### C. Diagrammatic derivation of XVCC

The equations of XVCC can be derived even more expediently with Feynman–Goldstone diagrams, which is equivalent to the normal-ordered second-quantized derivation; a diagram corresponds to all topologically related Wick’s contractions. Here, we introduce the rules of drawing closed diagrams for XVCC energy equations and open diagrams for XVCC amplitude equations as well as the rules of translating them into algebraic expressions. These rules differ from those for XVMP energies and self-energies.^{64} We illustrate their application to XVCC*m* with *m* = 1 (XVCCS) and *m* = 2 (XVCCSD) in this subsection.

We recast each term in the normal-ordered Hamiltonian, $H\u02c6N$, into diagrammatic forms, as shown in Fig. 1. Each diagram represents one of the normal-ordered operators defined by Eqs. (12)–(15). The cluster excitation operators, ${T\u02c6n}$, are expressed diagrammatically in Fig. 2.

Table I lists the rules to generate diagrams of the XVCC*m* energy and amplitude equations in an *n*th-order Taylor-series PES. Applying these rules, we find the energy diagrams of XVCCS in a QFF (*n* = 4) to be those in Fig. 3. The energy diagrams of XVCCSD in a QFF are the union of Figs. 3 and 4. Hence,

where the summands are the algebraic interpretation of the diagrams. Since each $T\u02c6k$ vertex in an energy diagram has only upgoing lines (signifying the fact that they describe an excitation process), all of its lines have no choice but to connect with the $H\u02c6N$ vertex, making the whole diagram always connected when closed, guaranteeing its thermodynamic extensivity as per the extensive diagram theorem.^{10,11,65}

(1) | Draw an $H\u02c6N$ vertex having up to n lines. |

(2) | Draw 0 to $nT\u02c6k$ vertices (1 ≤ k ≤ m) beneath the $H\u02c6N$ vertex. |

(3) | Connect the lines of the $T\u02c6k$ vertices with the lines of the $H\u02c6N$ |

vertex without altering (upgoing or downgoing) line directions, | |

leaving l external (i.e., dangling) lines to form a connected | |

diagram, which is either closed (l = 0) or open (l > 0). | |

(4) | Repeat steps (1)–(3) to enumerate all topologically distinct |

diagrams. If permutation of $T\u02c6$ vertices allows two diagrams | |

to overlap, they are deemed topologically equivalent. |

(1) | Draw an $H\u02c6N$ vertex having up to n lines. |

(2) | Draw 0 to $nT\u02c6k$ vertices (1 ≤ k ≤ m) beneath the $H\u02c6N$ vertex. |

(3) | Connect the lines of the $T\u02c6k$ vertices with the lines of the $H\u02c6N$ |

vertex without altering (upgoing or downgoing) line directions, | |

leaving l external (i.e., dangling) lines to form a connected | |

diagram, which is either closed (l = 0) or open (l > 0). | |

(4) | Repeat steps (1)–(3) to enumerate all topologically distinct |

diagrams. If permutation of $T\u02c6$ vertices allows two diagrams | |

to overlap, they are deemed topologically equivalent. |

Table II provides a set of rules to translate these diagrams into algebraic expressions. For example, diagram *E*_{4} of Fig. 3 can be labeled as shown in Fig. 5 using rules 1 through 5 of the table. Following the remaining rules, we find the algebraic expression of the whole diagram to be

where the 1/4! factor comes from rule 8 because there are four equivalent $T\u02c61$ vertices. It is straightforward to confirm that the algebraic interpretations of *E*_{1}, *E*_{2}, and *E*_{3} obtained with Table II agree exactly with the second, third, and fourth terms of the XVCCS energy expression [Eq. (28)], respectively.

(1) | Label each internal line with a mode index (i, j, k, …). |

An internal line is a line terminated by vertices at both ends. | |

(2) | Label each external line with a mode index (ν, μ, λ, …). |

An external line is a line terminated by a vertex at only one end. | |

(3) | Associate each $H\u02c6N$ vertex with a Hamiltonian amplitude W. |

The indices of W are given by the labels of the lines connected | |

to that vertex. A two-line $H\u02c6N$ vertex is associated with W^{−} | |

if both lines are upgoing or downgoing; with W^{+} otherwise. | |

(4) | Associate each $T\u02c6k$ vertex with an excitation amplitude τ. |

The indices of τ are given by the labels of the lines connected | |

to that vertex. | |

(5) | Associate each $R\u02c6k$ vertex with an excitation amplitude r. |

The indices of r are given by the labels of the lines connected | |

to that vertex. | |

(6) | Sum over all internal line indices (i, j, k, …). |

(7) | Sum over all permutations of external line indices (ν, μ, λ, …). |

(8) | For each set of n equivalent vertices, multiply by a factor of |

(1/n!). If permutation of two vertices leaves the diagram | |

unchanged topologically, the two vertices are equivalent. | |

(9) | For each set of m equivalent lines, multiply by a factor of |

(1/m!). If permutation of two lines leaves the diagram | |

unchanged topologically, the two lines are equivalent. |

(1) | Label each internal line with a mode index (i, j, k, …). |

An internal line is a line terminated by vertices at both ends. | |

(2) | Label each external line with a mode index (ν, μ, λ, …). |

An external line is a line terminated by a vertex at only one end. | |

(3) | Associate each $H\u02c6N$ vertex with a Hamiltonian amplitude W. |

The indices of W are given by the labels of the lines connected | |

to that vertex. A two-line $H\u02c6N$ vertex is associated with W^{−} | |

if both lines are upgoing or downgoing; with W^{+} otherwise. | |

(4) | Associate each $T\u02c6k$ vertex with an excitation amplitude τ. |

The indices of τ are given by the labels of the lines connected | |

to that vertex. | |

(5) | Associate each $R\u02c6k$ vertex with an excitation amplitude r. |

The indices of r are given by the labels of the lines connected | |

to that vertex. | |

(6) | Sum over all internal line indices (i, j, k, …). |

(7) | Sum over all permutations of external line indices (ν, μ, λ, …). |

(8) | For each set of n equivalent vertices, multiply by a factor of |

(1/n!). If permutation of two vertices leaves the diagram | |

unchanged topologically, the two vertices are equivalent. | |

(9) | For each set of m equivalent lines, multiply by a factor of |

(1/m!). If permutation of two lines leaves the diagram | |

unchanged topologically, the two lines are equivalent. |

Likewise, the application of the rules in Table II to the diagrams in Fig. 4 leads to the XVCCSD energy equation,

The factor of (1/2!)^{3} in the penultimate term (diagram *E*_{7}) arises from the fact that there are one pair of equivalent $T\u02c62$ vertices and two pairs of equivalent lines.

Table I also outlines rules to enumerate all diagrams in the XVCC*m* amplitude equations. Unlike energy diagrams, a diagram in the *T _{l}* amplitude equation is open with

*l*external lines. Such diagrams can be disconnected, in which case they must be deleted as per rule 3.

Figure 6 lists all four connected diagrams of the XVCCS *T*_{1} amplitude equation in a QFF, which can thus be written as

where *S _{i}* is the algebraic interpretation of the corresponding diagram. The

*T*

_{1}amplitude equation of XVCCSD in a QFF is diagrammatically expressed by Figs. 6 and 7. The

*T*

_{2}amplitude equation of the same method is defined with the diagrams in Fig. 8. The

*T*

_{1}and

*T*

_{2}amplitude equations of XVCCSD are thus formally written as

The interpretation rules in Table II apply to open diagrams in these amplitude equations also. For example, lines and vertices in diagram *S*_{4} can be labeled as shown in Fig. 9. Its algebraic expression then reads

Applying the same rules to diagrams *S*_{1}, *S*_{2}, and *S*_{3} and substituting the results in Eq. (46) reproduces the *T*_{1} amplitude equation [Eq. (39)] derived algebraically.

where $T1(1)$ is given algebraically in Eq. (39).

Similarly, the *T*_{2} amplitude equation of the same method can directly result from Fig. 8 and reads

where $P\u02c6(\nu 1\cdots \nu n)$ yields the sum of *n*! permutations of indices {*ν*_{1}, …, *ν _{n}*} in the term it acts on (cf. rule 7 of Table II). For instance,

In the same way, the diagrammatic and algebraic energy and amplitude equations of XVCC*m* with *m* = 3 (XVCCSDT) in a QFF are derived. They are documented in Appendix B. For 4 ≤ *m* ≤ 8, diagrams and equations were generated with the aid of computer algebra based on the rules provided in Tables I and II and are documented in the supplementary material.^{66}

### D. EOM-XVCC formalism

In a vibrational analysis, it is usually the frequencies (excitation energies) that are of greater interest than the zero-point energies (total energies in the ground state). Here, we stipulate and implement an excited-state XVCC theory using the EOM formalism,^{67–72} which is equivalent to CC linear response^{73–78} for excitation energies.

The EOM-XVCC*m* theory can be viewed as a CI for excited states using the XVCC effective Hamiltonian $H\u0304N$ defined as

where $T\u02c6$ is truncated after the *m*-fold cluster excitation operator [Eq. (3)] and all of its amplitudes should have been determined by a preceding XVCC*m* calculation.

Let $R\u02c6$ be a linear excitation operator truncated after the same rank,

with

where *r*_{i1i2⋯in} is an excitation amplitude. We require that the excited-state wave function $R\u02c6|\Phi 0\u3009$ satisfy the Schrödinger equation with the modified Hamiltonian $H\u0304N$ in the space of Hartree products reachable by $R\u02c6$ from |Φ_{0}〉. That is,

for all *ν*_{1}, …, *ν _{n}* and all

*n*(1 ≤

*n*≤

*m*), where $\Delta Ek(m)=Ek(m)\u2212E0(0)$ and $Ek(m)$ is the total EOM-XVCC

*m*energy of the

*k*th excited state. The constant term

*r*

_{0}is formally necessary in Eq. (56) for excited states of the same symmetry as the ground state. Although $R\u02c6$ and its amplitudes vary with

*k*and should thus carry a

*k*identifier, it is suppressed for brevity.

Equations (59) and (60) define EOM-XVCC*m* as an eigenvalue problem with a non-Hermitian matrix of $H\u0304N$, which thus has both left and right eigenvectors with $R\u02c6$ corresponding to the latter. They can, however, be further simplified by using the fact that *T* amplitudes in $H\u0304N$ satisfy the XVCC*m* equations. Equations (29) and (40) indicate that |Φ_{0}〉 satisfies the Schrödinger equation with $H\u0304N$ in the same Hartree-product space. This implies

because $R\u02c6\u2020|\Phi 0\u3009\u222a{R\u02c6\u2020|\Phi \nu 1\nu 2\cdots \nu n\u3009}$ spans the space no greater than |Φ_{0}〉∪{|Φ_{ν1ν2⋯νn}〉}. Subtracting these from Eqs. (59) and (60), we obtain

where $\omega k(m)=\Delta Ek(m)\u2212\Delta E0(m)$ and is the EOM-XVCC*m* transition energy (frequency) to the *k*th excited state.

Equation (64) is the ansatz for EOM-XVCC*m* for frequencies, which can also be written as

for all *ν*_{1}, …, *ν _{n}* and all

*n*(1 ≤

*n*≤

*m*). We call this the

*R*amplitude equation of EOM-XVCC

_{n}*m*. The nested connectedness requirement is always met when the term is simply connected because both $T\u02c6$ and $R\u02c6$ are excitation operators and can only directly connect with $H\u02c6N$, but not with each other. The connectedness requirement, however, deletes all terms involving

*r*

_{0}because the latter has no index (or a line in its diagrammatic representation) to be shared with other factors.

Equation (65), therefore, refers to a set of equations whose number is equal to the number of unknowns (excluding *r*_{0}) minus one (because $\omega k(m)$ is also unknown) and, with an appropriate normalization condition, can determine $\omega k(m)$ and all the amplitudes in $R\u02c6$ except *r*_{0}. The value of *r*_{0}, in turn, can be determined by Eq. (63), though it is not necessary if only frequencies are sought. Equation (65) satisfies the intensive diagram theorem,^{10,11} ensuring the size consistency of EOM-XVCC*m* for total excited-state energies and frequencies, despite its CI-like structure.

Just like the electronic counterpart, EOM-XVCC*m* can also be derived by applying the first-order time-dependent perturbation (linear-response) theory to XVCC*m*. This derivation (not repeated here because it is identical to the electronic case^{78}) elegantly explains why the truncation ranks of $T\u02c6$ and $R\u02c6$ should be the same. EOM-XVCC*m* with the XVSCF[*n*] reference is the same as excited-state VCC of Prasad and coworkers.^{27,32}

### E. Algebraic derivation of EOM-XVCC

Here, we illustrate the normal-ordered second-quantized derivation of the amplitude equations of EOM-XVCC*m* using *m* = 1 (EOM-XVCCS) in a QFF as an example. The working equations of EOM-XVCC*m* with *m* = 2 (EOM-XVCCSD) are documented in Appendix C.

The left-hand side of Eq. (65) in the case of *m* = 1 reads

where the nested connectedness means simple connectedness in this case and terms containing higher powers of $T\u02c61$ cannot yield any connected contributions.

The first term of Eq. (66) has only one full contraction,

which is connected. The next term has two connected full contractions,

The last term yields the connected full contraction,

Together, we obtain the *R*_{1} amplitude equation of EOM-XVCCS in a QFF as

### F. Diagrammatic derivation of EOM-XVCC

The diagrammatic forms of the $R\u02c6n$ operators are given in Fig. 10. They are isomorphic to the diagrams of the $T\u02c6n$ operators with the only difference being the colors of the vertices. The rules to generate all diagrams entering the amplitude equations of EOM-XVCC*m* are given in Table III. The rules for interpreting these diagrams algebraically are provided in Table II.

(1) | Draw an $H\u02c6N$ vertex having up to n lines. |

(2) | Draw a $R\u02c6h$ vertex (1 ≤ h ≤ m) and zero to $(n\u22121)T\u02c6k$ vertices |

(1 ≤ k ≤ m) beneath the $H\u02c6N$ vertex. | |

(3) | Connect the lines of the $R\u02c6h$ and $T\u02c6k$ vertices with the lines of |

the $H\u02c6N$ vertex without altering (upgoing or downgoing) line | |

directions, leaving l external (i.e., dangling) lines to form | |

a connected diagram. | |

(4) | Repeat steps (1)–(3) to enumerate all topologically distinct |

diagrams. If permutation of $T\u02c6$ vertices allows two | |

diagrams to overlap, they are deemed topologically equivalent. |

(1) | Draw an $H\u02c6N$ vertex having up to n lines. |

(2) | Draw a $R\u02c6h$ vertex (1 ≤ h ≤ m) and zero to $(n\u22121)T\u02c6k$ vertices |

(1 ≤ k ≤ m) beneath the $H\u02c6N$ vertex. | |

(3) | Connect the lines of the $R\u02c6h$ and $T\u02c6k$ vertices with the lines of |

the $H\u02c6N$ vertex without altering (upgoing or downgoing) line | |

directions, leaving l external (i.e., dangling) lines to form | |

a connected diagram. | |

(4) | Repeat steps (1)–(3) to enumerate all topologically distinct |

diagrams. If permutation of $T\u02c6$ vertices allows two | |

diagrams to overlap, they are deemed topologically equivalent. |

As an illustration, we derive the *R*_{1} amplitude equation of EOM-XVCCS in a QFF diagrammatically. Diagrams in this equation are drawn in Fig. 11. The *R*_{1} amplitude equation can thus be formally written as

Diagram $S3\u2032$ may be labeled as shown in Fig. 12 according to the rules in Table II. The corresponding algebraic translation is

where the 1/2! factor is due to a pair of equivalent $T\u02c61$ vertices. The other two diagrams can be interpreted similarly, yielding the same *R*_{1} amplitude equation as Eq. (70).

The diagrammatic and algebraic equations of EOM- XVCC*m* for *m* = 2 through 8 were generated with the aid of computer algebra based on the rules provided in Tables II and III. The EOM-XVCCSD equations are documented in Appendix C, while the EOM-XVCC*m* equations for 3 ≤ *m* ≤ 8 are documented in the supplementary material.^{66}

## III. AUTOMATED COMPUTER IMPLEMENTATION

The XVCC*m* and EOM-XVCC*m* methods with *m* = 1 through 8 in a QFF were implemented. The Hamiltonian employed was the pure vibrational one, neglecting any coupling between vibrations and rotations. The XVSCF(4), XVSCF[4], or the harmonic reference wave function was used (“4” standing for the QFF). The choice of the reference wave function is indicated by a suffix as in XVCC*m* (4), XVCC*m*[4], or XVCC*m*{4}, respectively; “{*n*}” suffix indicates the use of a harmonic reference with the *n*th-order PES.

The *T _{n}* amplitude equations of XVCC

*m*were solved with a standard Newton–Krylov algorithm,

^{79}implemented in igmres, available from scipy. The initial values of

*τ*were zero. The

*R*amplitude equations of EOM-XVCC

_{n}*m*with

*m*≥ 2 were solved with the iterative Arnoldi method, a non-Hermitian generalization of the Lanczos method,

^{80–82}available in arpack.

^{83}The eigenvalue convergence criterion was set to machine precision. The Arnoldi method does not require the matrix elements of $H\u0304N$ stored in memory; instead, it evaluates the left-hand sides of the

*R*amplitude equations [Eq. (65)] directly for several trial

_{n}*r*vectors. This allows a small subset of low-lying excited-state roots to be determined in a memory-efficient manner. On the other hand, the EOM-XVCCS problem was solved by diagonalization of the $H\u0304N$ matrix, solving Eqs. (59) and (60) for all excited-state energies instead.

The XVCC*m* and EOM-XVCC*m* diagrammatic equations were derived by a symbolic algebra program using the rules presented in Tables I and III. These diagrammatic equations were then translated by the same symbolic algebra program into the algebraic form according to the interpretation rules given in Table II. Terms that originate from a single diagram and are thus related by a permutation operator were identified and consolidated. When the resulting expressions were evaluated literally, the operation costs of XVCC*m* and EOM-XVCC*m* in an *n*th-order Taylor-series PES would scale as steeply as *O*(*N*^{m+n}), where *N* is the number of modes. To reduce the cost, three types of algorithm optimizations were performed computationally before the automatic synthesis of c++ codes was carried out.^{84,85} The optimizations considered in this work were the strength reduction, intermediate reuse, and factorization.^{84} The details of each optimization procedure are given below.

Every tensor has a high degree of index permutation symmetry such as *τ*_{i1i2i3} = *τ*_{i1i3i2} = *τ*_{i2i1i3} = *τ*_{i2i3i1} = *τ*_{i3i1i2} = *τ*_{i3i2i1}. Exploiting it can compress the memory spaces of the *m*-fold cluster and linear excitation amplitudes and of *m*th-order force constants by a factor of *m*! and the operation cost of their contractions to an even greater degree. Nevertheless, the compression factor of *m*! in XVCC*m* is much smaller than (*m*!)^{2} in electronic CC theory with connected *m*-fold excitations, and our preliminary implementation of the algorithm that takes into account this symmetry did not demonstrate any gain in computational speed owing to an enormous cost overhead associated with complicated data layout. Therefore, in this study, we neglected this symmetry.

### A. Strength reduction

The algebraic equations of XVCC and EOM-XVCC are sums of contractions of tensors (force constants *W*, *T* amplitudes, and *R* amplitudes). The operation cost of evaluating a contraction (product) of more than two tensors is greatly reduced by carrying out a series of binary tensor contractions. Suppose a tensor contraction *E* = *ABCD*, where *A* through *E* are tensors of some dimension (including possibly scalars). This product can be evaluated stepwise as

where *I*_{1} and *I*_{2} are the so-called *intermediate* tensors. Alternatively, the same product may be evaluated as

Both orders of contraction give the same final result (*E*) but at generally different operation and memory costs. Strength reduction determines the *sequential* contraction order with the minimal operation cost for each product by blanket search; it neglects *parallel* contraction such as (*AB*)(*CD*).

The result of this procedure is depicted in Fig. 13 for diagram *T*_{17}, which occurs in the *T*_{3} amplitude equation of XVCCSDT (see Appendix B). It is algebraically defined as

A literal evaluation of this equation involves an *O*(*N*^{7}) arithmetic operations, as it must perform *N*^{4} multiplications and additions to evaluate each of *N*^{3} elements of *T*_{17}, where *N* is the number of modes. However, the same result can be obtained at an *O*(*N*^{4}) cost by breaking the quaternary tensor contraction into three separate binary tensor contractions carried out in the following order. First, the force-constant tensor is multiplied with *τ _{l}*, which has no external index,

This step involves *N*^{4} multiplications and additions that define an intermediate, *T*_{17a}, of dimension *N*^{3}. This intermediate is then multiplied by *τ*_{jkν} at an *O*(*N*^{4}) cost to define another intermediate, *T*_{17b}, of dimension *N*^{2},

Finally, *T*_{17} is obtained by contraction of intermediate *T*_{17b} with *τ*_{iμλ} at an *O*(*N*^{4}) cost,

Therefore, the cost of evaluating diagram *T*_{17} is reduced from *O*(*N*^{7}) to *O*(*N*^{4}).

### B. Intermediate reuse

This optimization simply stores an intermediate appearing two or more times in the whole equation and reuses it without recalculating. Consider the following sum-of-product expression of tensors: *E* = *A*(*BC*) + *D*(*BC*). If the strength reduction procedure suggests the optimal contraction orders specified by the parentheses, intermediate *I*_{1} = *BC* is calculated only once and stored for reuse. Hence, the whole sum is evaluated as

The numbers of multiplications and additions as a function of *N* required to evaluate all the amplitude equations of XVCC*m* and EOM-XVCC*m* with no optimization and with strength reduction and intermediate reuse are compared in Table IV. For all methods with *m* ≥ 2, strength reduction leads to a great reduction in the operation cost from *O*(*N*^{m+4}) to *O*(*N*^{m+2}). Generally, for an *n*th-order PES, the optimal scaling of XVCC*m* and EOM-XVCC*m* is *O*(*N*^{m+⌊n/2⌋}) when *m* ≥ *n* and is no greater than *O*(*N*^{n+⌊(m−1)/2⌋}) when *m* < *n*. A proof of this statement is given in Appendix D.

Method . | No optimization . | Strength reduction . | Factorization^{a}
. |
---|---|---|---|

XVCCS | 5N^{4} + 4N^{3} + 2N^{2} | 2N^{4} + 4N^{3} + 8N^{2} | 2N^{4} + 6N^{3} + 4N^{2} |

XVCCSD | 11N^{6} + 31N^{5} + 34N^{4} | 13N^{4} + 54N^{3} + 36N^{2} | 13N^{4} + 44N^{3} + 25N^{2} |

XVCCSDT | 60N^{7} + 119N^{6} + 110N^{5} | 11N^{5} + 165N^{4} + 96N^{3} | 11N^{5} + 86N^{4} + 70N^{3} |

XVCCSDTQ | 241N^{8} + 360N^{7} + 282N^{6} | 24N^{6} + 425N^{5} + 225N^{4} | 24N^{6} + 172N^{5} + 109N^{4} |

XVCC5 | 799N^{9} + 1053N^{8} + 720N^{7} | 104N^{7} + 1061N^{6} + 488N^{5} | 82N^{7} + 342N^{6} + 194N^{5} |

XVCC6 | 2577N^{10} + 3086N^{9} + 1888N^{8} | 175N^{8} + 3079N^{7} + 1134N^{6} | 143N^{8} + 1135N^{7} + 366N^{6} |

XVCC7 | 8398N^{11} + 9246N^{10} + 5070N^{9} | 285N^{9} + 8902N^{8} + 3156N^{7} | 241N^{9} + 4392N^{8} + 1161N^{7} |

XVCC8 | 28 099N^{12} + 28 380N^{11} + 13 989N^{10} | 468N^{10} + 26 641N^{9} + 8991N^{8} | 410N^{10} + 18 021N^{9} + 4420N^{8} |

EOM-XVCCS | 5N^{4} + 3N^{3} + 2N^{2} | 2N^{4} + 4N^{3} + 7N^{2} | 2N^{4} + 4N^{3} + 7N^{2} |

EOM-XVCCSD | 30N^{6} + 70N^{5} + 54N^{4} | 21N^{4} + 97N^{3} + 63N^{2} | 21N^{4} + 57N^{3} + 52N^{2} |

EOM-XVCCSDT | 183N^{7} + 286N^{6} + 220N^{5} | 15N^{5} + 343N^{4} + 189N^{3} | 15N^{5} + 131N^{4} + 99N^{3} |

EOM-XVCCSDTQ | 729N^{8} + 916N^{7} + 598N^{6} | 34N^{6} + 965N^{5} + 474N^{4} | 34N^{6} + 311N^{5} + 170N^{4} |

EOM-XVCC5 | 2490N^{9} + 2794N^{8} + 1614N^{7} | 207N^{7} + 2649N^{6} + 1096N^{5} | 143N^{7} + 696N^{6} + 356N^{5} |

EOM-XVCC6 | 8295N^{10} + 8533N^{9} + 4432N^{8} | 348N^{8} + 8164N^{7} + 2782N^{6} | 254N^{8} + 2527N^{7} + 731N^{6} |

EOM-XVCC7 | 27 975N^{11} + 26 506N^{10} + 12 457N^{9} | 569N^{9} + 24 962N^{8} + 8309N^{7} | 439N^{9} + 11 231N^{8} + 2570N^{7} |

EOM-XVCC8 | 96 603N^{12} + 83 992N^{11} + 35 934N^{10} | 934N^{10} + 78 522N^{9} + 25 129N^{8} | 762N^{10} + 51 289N^{9} + 11 274N^{8} |

Method . | No optimization . | Strength reduction . | Factorization^{a}
. |
---|---|---|---|

XVCCS | 5N^{4} + 4N^{3} + 2N^{2} | 2N^{4} + 4N^{3} + 8N^{2} | 2N^{4} + 6N^{3} + 4N^{2} |

XVCCSD | 11N^{6} + 31N^{5} + 34N^{4} | 13N^{4} + 54N^{3} + 36N^{2} | 13N^{4} + 44N^{3} + 25N^{2} |

XVCCSDT | 60N^{7} + 119N^{6} + 110N^{5} | 11N^{5} + 165N^{4} + 96N^{3} | 11N^{5} + 86N^{4} + 70N^{3} |

XVCCSDTQ | 241N^{8} + 360N^{7} + 282N^{6} | 24N^{6} + 425N^{5} + 225N^{4} | 24N^{6} + 172N^{5} + 109N^{4} |

XVCC5 | 799N^{9} + 1053N^{8} + 720N^{7} | 104N^{7} + 1061N^{6} + 488N^{5} | 82N^{7} + 342N^{6} + 194N^{5} |

XVCC6 | 2577N^{10} + 3086N^{9} + 1888N^{8} | 175N^{8} + 3079N^{7} + 1134N^{6} | 143N^{8} + 1135N^{7} + 366N^{6} |

XVCC7 | 8398N^{11} + 9246N^{10} + 5070N^{9} | 285N^{9} + 8902N^{8} + 3156N^{7} | 241N^{9} + 4392N^{8} + 1161N^{7} |

XVCC8 | 28 099N^{12} + 28 380N^{11} + 13 989N^{10} | 468N^{10} + 26 641N^{9} + 8991N^{8} | 410N^{10} + 18 021N^{9} + 4420N^{8} |

EOM-XVCCS | 5N^{4} + 3N^{3} + 2N^{2} | 2N^{4} + 4N^{3} + 7N^{2} | 2N^{4} + 4N^{3} + 7N^{2} |

EOM-XVCCSD | 30N^{6} + 70N^{5} + 54N^{4} | 21N^{4} + 97N^{3} + 63N^{2} | 21N^{4} + 57N^{3} + 52N^{2} |

EOM-XVCCSDT | 183N^{7} + 286N^{6} + 220N^{5} | 15N^{5} + 343N^{4} + 189N^{3} | 15N^{5} + 131N^{4} + 99N^{3} |

EOM-XVCCSDTQ | 729N^{8} + 916N^{7} + 598N^{6} | 34N^{6} + 965N^{5} + 474N^{4} | 34N^{6} + 311N^{5} + 170N^{4} |

EOM-XVCC5 | 2490N^{9} + 2794N^{8} + 1614N^{7} | 207N^{7} + 2649N^{6} + 1096N^{5} | 143N^{7} + 696N^{6} + 356N^{5} |

EOM-XVCC6 | 8295N^{10} + 8533N^{9} + 4432N^{8} | 348N^{8} + 8164N^{7} + 2782N^{6} | 254N^{8} + 2527N^{7} + 731N^{6} |

EOM-XVCC7 | 27 975N^{11} + 26 506N^{10} + 12 457N^{9} | 569N^{9} + 24 962N^{8} + 8309N^{7} | 439N^{9} + 11 231N^{8} + 2570N^{7} |

EOM-XVCC8 | 96 603N^{12} + 83 992N^{11} + 35 934N^{10} | 934N^{10} + 78 522N^{9} + 25 129N^{8} | 762N^{10} + 51 289N^{9} + 11 274N^{8} |

^{a}

Factorization leads to a cost increase in XVCCS and thus was turned off in its implementation.

### C. Factorization

The operation cost of XVCC and EOM-XVCC can be further reduced by factorization. Let us consider the sum-of-product tensor equation, *E* = *A*(*BC*) + *A*(*BD*), where the parentheses indicate the optimal order of contraction for each product. This sum can be evaluated as

Two contractions by the common tensor *B* in *I*_{2} = *BC* + *BD* are factored as *I*_{2} = *B*(*C* + *D*), reducing the multiplication cost by half. Also, two contractions by *A* are also merged into one, further halving the multiplication cost. This procedure should not be confused with intermediate reuse.

Factorability of equations depends on the contraction order of each product, making strength reduction and factorization a coupled optimization problem. Since a complete solution of this coupled problem is exponentially complex, we approach it as two separate, sequential optimizations: strength reduction is performed first for each product, whereupon factorization is carried out by blanket search without altering the order of contraction. Hence, a sum of intermediates multiplied by the same tensor and the same permutation operator is factored. The resulting new intermediates are sums of tensors and they are recursively factored. In this work, the permutation operator itself is not decomposed into a product of smaller permutation operators, unlike in previous work for electronic CC.^{84}

An example of the procedure is given below and drawn in Fig. 14,

where the parentheses specify the optimal order of contraction for each product determined by the preceding strength reduction step. Factorization does not alter these orders. When evaluating the sum of these five terms, the common factor *τ*_{iν} is pulled out,

where the intermediate *D*_{6+8+9+17+23} is the sum of three intermediates,

with the subscripts indicating the diagrammatic origin of the intermediates. The first intermediate does not lend itself to any further factorization and is given by

The second intermediate is further factored by the common tensor *τ _{j}* and is evaluated as

The third intermediate is also further factored as

The operation costs of XVCC*m* and EOM-XVCC*m* after factorization (in addition to strength reduction and intermediate reuse) are listed in Table IV. It should be understood that intermediate reuse (combined with factorization) applies to common intermediates appearing two or more times that are either a product or a sum-of-product of tensors. Since factorization is more likely for the final contraction step in our implementation, it is especially effective when the final step is the most costly. Factorization does not alter the overall *scaling* of the cost function, but generally reduces the prefactor multiplying the leading-order term, while it may increase or decrease prefactors of sub-leading-order terms.

Factorization is more important for EOM-XVCC and for higher values of *m*. In particular, for *m* > 4, factorization lowers the prefactors on the leading-order term of the cost functions. For *m* = 3 or 4, only the prefactors of the second and third leading-order terms of the cost functions are reduced. For *m* = 1 and 2, the cost is dominated by the contraction of the force-constant tensor with a tensor of *T* amplitudes, which occurs first in each product (see Appendix D); therefore, factorization is found least effective there.

### D. Observed computational cost

In Fig. 15, the CPU time spent in one cycle of the iterative solution of the XVCC*m* or EOM-XVCC*m* equation (3 ≤ *m* ≤ 5) in a model QFF is plotted as a function of the number of modes *N*. Each plot falls accurately on the asymptotic cost-scaling line of *O*(*N*^{m+2}) for *N* ≥ 20, suggesting that the leading-order terms in the cost function dominate the cost. For instance, the XVCCSDTQ and EOM-XVCCSDTQ methods with the aforementioned optimizations are remarkably efficient, taking only on the order of ten seconds per iteration for *N* = 35, despite the fact that their algorithms do not exploit the index permutation symmetry. These costs are expected to be negligible as compared with the cost associated with generating all higher-order force constants at the electronic structure level appropriate for these accurate vibrational methods.

## IV. NUMERICAL TESTS

The XVCC*m* and EOM-XVCC*m* methods were applied to the vibrational zero-point energies and frequencies of the water and formaldehyde molecules. Their equilibrium structure, normal modes, and QFF’s were calculated at the MP2/aug-cc-pVTZ level using nwchem ^{86} and sindo.^{87} The results were compared with those obtained from full VCI calculations with mavi ^{88} using the 20 lowest-lying harmonic-oscillator wave functions of each mode as the basis set. The VCI results were converged to within 0.1 cm^{−1} of the exact solutions for the QFF. Our results were also compared with the XVMP2(4) and XVMP2[4] results,^{59} obtained either in the frequency-independent approximation to the self-energy or by solving the Dyson equation with the frequency-dependent self-energy.

### A. Water

Table V compiles the errors in the XVCC*m* zero-point energy (*E*_{0}) and EOM-XVCC*m* frequencies (*ν*) from the VCI results for the water molecule. These errors decrease nearly monotonically with the excitation rank *m*, reaching the accuracy of a few cm^{−1} at *m* = 3 (CCSDT) for both *E*_{0} and *ν*. At *m* = 4 (CCSDTQ), the results are generally within tenths of cm^{−1} of the VCI results with a few exceptions. The results for *m* ≥ 5 are given in Table VI and are converged with accuracy higher than necessary, e.g., within 0.001 cm^{−1} at *m* = 8.

Method^{a}
. | E_{0}
. | ν_{1}
. | ν_{2}
. | ν_{3}
. |
---|---|---|---|---|

VCI | 4642.4 | 3682.4 | 1556.6 | 3791.5 |

Harmonic | 108.4 | 139.7 | 71.7 | 156.4 |

XVSCF(4) | 104.5 | 280.2 | −17.9 | 289.4 |

XVSCF[4] | 21.5 | 105.3 | 11.7 | 116.1 |

XVCCS{4} | 21.9 | 196.5 | −1.2 | 206.7 |

XVCCS(4) | 26.6 | 197.2 | −2.4 | 207.0 |

XVCCS[4] | 21.5 | 105.3 | 11.7 | 116.1 |

XVCCSD{4} | 22.1 | 33.6 | 7.9 | 34.3 |

XVCCSD(4) | 23.8 | 38.3 | 8.5 | 38.3 |

XVCCSD[4] | 21.5 | 26.0 | 6.8 | 26.6 |

XVCCSDT{4} | 0.2 | −1.7 | 2.2 | −2.0 |

XVCCSDT(4) | 0.1 | −2.4 | 2.2 | −2.7 |

XVCCSDT[4] | 0.4 | −1.7 | 2.4 | −1.9 |

XVCCSDTQ{4} | 0.0 | −0.7 | 0.2 | −0.7 |

XVCCSDTQ(4) | −0.0 | −1.0 | 0.1 | −1.0 |

XVCCSDTQ[4] | 0.0 | −0.1 | 0.1 | −0.1 |

XVMP2(4)^{b} | 4.5 | −24.6 | 3.8 | −28.3 |

XVMP2[4]^{b} | −1.3 | −18.9 | 2.4 | −20.3 |

XVMP2(4)^{c} | 4.5 | −14.3 | 3.8 | −20.8 |

XVMP2[4]^{c} | −1.3 | −14.8 | 2.5 | −17.0 |

Method^{a}
. | E_{0}
. | ν_{1}
. | ν_{2}
. | ν_{3}
. |
---|---|---|---|---|

VCI | 4642.4 | 3682.4 | 1556.6 | 3791.5 |

Harmonic | 108.4 | 139.7 | 71.7 | 156.4 |

XVSCF(4) | 104.5 | 280.2 | −17.9 | 289.4 |

XVSCF[4] | 21.5 | 105.3 | 11.7 | 116.1 |

XVCCS{4} | 21.9 | 196.5 | −1.2 | 206.7 |

XVCCS(4) | 26.6 | 197.2 | −2.4 | 207.0 |

XVCCS[4] | 21.5 | 105.3 | 11.7 | 116.1 |

XVCCSD{4} | 22.1 | 33.6 | 7.9 | 34.3 |

XVCCSD(4) | 23.8 | 38.3 | 8.5 | 38.3 |

XVCCSD[4] | 21.5 | 26.0 | 6.8 | 26.6 |

XVCCSDT{4} | 0.2 | −1.7 | 2.2 | −2.0 |

XVCCSDT(4) | 0.1 | −2.4 | 2.2 | −2.7 |

XVCCSDT[4] | 0.4 | −1.7 | 2.4 | −1.9 |

XVCCSDTQ{4} | 0.0 | −0.7 | 0.2 | −0.7 |

XVCCSDTQ(4) | −0.0 | −1.0 | 0.1 | −1.0 |

XVCCSDTQ[4] | 0.0 | −0.1 | 0.1 | −0.1 |

XVMP2(4)^{b} | 4.5 | −24.6 | 3.8 | −28.3 |

XVMP2[4]^{b} | −1.3 | −18.9 | 2.4 | −20.3 |

XVMP2(4)^{c} | 4.5 | −14.3 | 3.8 | −20.8 |

XVMP2[4]^{c} | −1.3 | −14.8 | 2.5 | −17.0 |

^{a}

XVCC for zero-point energy and EOM-XVCC for frequencies. Suffix “{4}” means the use of a harmonic reference, “(4)” means the use of a XVSCF(4) reference, and “[4]” means the use of a XVSCF[4] reference.

^{b}

In the frequency-independent approximation.^{59}

^{c}

Solutions of the Dyson equation with the frequency-dependent self-energy.^{59}

Method^{a}
. | E_{0}
. | ν_{1}
. | ν_{2}
. | ν_{3}
. |
---|---|---|---|---|

XVCC5{4} | −0.02 | −0.05 | 0.02 | −0.06 |

XVCC5(4) | −0.03 | −0.12 | 0.01 | −0.12 |

XVCC5[4] | −0.03 | 0.00 | 0.02 | 0.00 |

XVCC6{4} | −0.001 | 0.004 | 0.004 | 0.005 |

XVCC6(4) | −0.001 | 0.005 | 0.002 | 0.006 |

XVCC6[4] | −0.001 | 0.006 | 0.003 | 0.007 |

XVCC7{4} | −0.0004 | 0.003 | 0.001 | 0.003 |

XVCC7(4) | −0.0002 | 0.004 | 0.001 | 0.005 |

XVCC7[4] | −0.0005 | −0.001 | 0.002 | −0.001 |

XVCC8{4} | 0 | 0 | 0 | 0 |

XVCC8(4) | 0 | 0 | 0 | 0 |

XVCC8[4] | 0 | 0 | 0 | 0 |

Method^{a}
. | E_{0}
. | ν_{1}
. | ν_{2}
. | ν_{3}
. |
---|---|---|---|---|

XVCC5{4} | −0.02 | −0.05 | 0.02 | −0.06 |

XVCC5(4) | −0.03 | −0.12 | 0.01 | −0.12 |

XVCC5[4] | −0.03 | 0.00 | 0.02 | 0.00 |

XVCC6{4} | −0.001 | 0.004 | 0.004 | 0.005 |

XVCC6(4) | −0.001 | 0.005 | 0.002 | 0.006 |

XVCC6[4] | −0.001 | 0.006 | 0.003 | 0.007 |

XVCC7{4} | −0.0004 | 0.003 | 0.001 | 0.003 |

XVCC7(4) | −0.0002 | 0.004 | 0.001 | 0.005 |

XVCC7[4] | −0.0005 | −0.001 | 0.002 | −0.001 |

XVCC8{4} | 0 | 0 | 0 | 0 |

XVCC8(4) | 0 | 0 | 0 | 0 |

XVCC8[4] | 0 | 0 | 0 | 0 |

^{a}

XVCC for zero-point energy and EOM-XVCC for frequencies. Suffix “{4}” means the use of a harmonic reference, “(4)” means the use of a XVSCF(4) reference, and “[4]” means the use of a XVSCF[4] reference.

Note that the VCI results in this table were obtained with 20 basis functions per mode. The number of Hartree products involved in such a calculation grows exponentially with size or as *O*(*b ^{N}*), where

*b*is the number of basis functions and

*N*is the number of modes. The numbers of the

*T*and

*R*amplitudes in XVCC

*m*and EOM-XVCC

*m*, on the other hand, increase only polynomially or as

*O*(

*N*). The values of

^{m}*b*for VCI and

^{N}*N*for XVCCSDTQ (

^{m}*N*= 3,

*b*= 20, and

*m*= 4) are 8000 and 81, respectively, implying a much severer limitation in VCI’s applicability to larger molecules. The number of unknowns in Christiansen’s VCC (Ref. 36) grows as

*O*(

*b*). Therefore, it too suffers from the consequence of having a finite number of basis functions, unlike XVCC

^{m}*m*, which is inherently an infinite-basis method.

The accuracy of *E*_{0} from XVCC*m* is relatively insensitive to the choice of the reference wave function, despite much larger differences in the results of the reference methods themselves: harmonic approximation, XVSCF(4), and XVSCF[4]. This is in line with a similar observation^{72} made in the electronic CC results for correlation energies, which is understood to be due to the single-excitation operator which corrects any deficiency in the reference wave function. Likewise, the XVCCS(4) and XVCCS[4] results of *E*_{0} are close to each other, while the XVSCF(4) results differ considerably from those of XVSCF[4]. This difference is due to the anharmonic shift in the geometry, neglected in XVSCF(4) but included in XVSCF[4]. Diagrammatically, this contribution is largely represented by *E*_{1} of Fig. 3. The fact that this diagram is found among the XVCCS energy equation means that the effect of anharmonic geometry shift is included in XVCCS regardless of the reference wave function.

XVSCF[4] determines and adopts the first-order Dyson geometry^{50} as the center of coordinates, where *W _{i}* = 0 for all

*i*. The

*T*

_{1}amplitude equation of XVCCS[4], Eq. (39), then has a trivial solution of

*τ*= 0 for all

_{i}*i*. This renders the XVCCS correction, $E0(1)\u2212E0(0)$ in Eq. (28), zero, making the XVCCS[4] and XVSCF[4] values of

*E*

_{0}identical to each other. This expected behavior is observed in Table V. Furthermore, substituting

*τ*= 0 into the

_{i}*R*

_{1}amplitude equation of EOM-XVCCS [Eq. (70)], we find the EOM-XVCCS[4] frequencies to be the eigenvalues of the $Wij+$ matrix, whose eigenvectors are known to define the first-order Dyson coordinates.

^{50,52}These eigenvalues are usually, but not always, close to the XVSCF[4] frequencies,

^{50}the latter being the diagonal elements ($Wii+$) of the same matrix. This is also observed in Table V. Therefore, XVCCS can take into account the leading-order effects in

*E*

_{0}and

*ν*from anharmonic corrections both in geometry and coordinates.

It may be noted that XVMP2 performs considerably better than XVCCSD for *E*_{0}. This is in striking contrast with the electronic counterparts; CCSD has been observed to consistently outperform MP2.^{1,2} This difference can be traced to the difference in the rank of the Hamiltonian in the two cases: the electronic Hamiltonian is a two-body operator, whereas the vibrational one is an *n*-body operator, where *n* is the rank of the Taylor-series PES (*n* = 4 in this calculation). The most important anharmonic correction to *E*_{0} is diagram *E*_{9} in Fig. 16 (in Appendix B), which involves cubic force constants and first appears at XVCCSDT. This contribution is to a great extent accounted for by the isomorphic diagram in XVMP2 labeled **2C′** in Ref. 59. This explains why only at *m* ≥ 3 does XVCC*m* outperform XVMP2. It is generally recommended to set the value of *m* equal to *n*, if not higher, because the overall cost of the calculation is dominated by the *O*(*N ^{n}*) step of the force-constant evaluation.

The rather small change in *E*_{0} from XVCCS to XVCCSD can also be understood as follows. The chief contribution of XVCCSD to *E*_{0}, diagram *E*_{5} of Fig. 4, is nearly zero in the XVSCF reference, where $Wii\u2212=0$ holds.^{52} When the normal coordinates are close to the first-order Dyson coordinates, which is usually the case, $Wij\u2212=Wij+\u22480$ is also true for *i*≠*j*, making the remaining (i.e., off-diagonal) contributions to *E*_{5} insignificant. Furthermore, the smallness of $Wij\u2212$ implies the smallness of the *T*_{2} amplitudes [cf. Eq. (53)]. Although XVCCSD uses cubic force constants in *E*_{6}, they are multiplied by these small *T*_{2} amplitudes and can hardly alter *E*_{0} of XVCCS. The frequencies from EOM-XVCCSD, on the other hand, are considerably better than EOM-XVCCS and are similar to those of XVMP2. This may be ascribed to diagram $S5\u2032$ in Fig. 20 (in Appendix C) included in EOM-XVCCSD, which can describe the same essential cubic corrections accounted for by self-energy diagrams **2e**′ and **2f**′ (Ref. 59) of XVMP2.

For the zero-point energies and frequencies that are not affected by any resonance, we observe the following general order of accuracy: harmonic ≈ XVSCF(*n*) < XVCCS{*n*} ≈ XVCCS(*n*) < XVSCF[*n*] ≈ XVCCS[*n*] < XVCCSD < XVMP2 < XVCCSDT < XVCCSDTQ ≤ XVCC5 ≤ XVCC6 ≤ XVCC7 ≤ XVCC8, which differs from the order established in electronic-structure theory.^{1,2} In a QFF (*n* = 4), XVCCSDTQ or *m* = 4 already achieves accuracy high enough for practical purposes and XVCC*m* with *m* > *n* seem to have diminishing returns. Here, XVCC stands for both XVCC and EOM-XVCC and the effects of the reference wave functions and the frequency-independent approximation in XVMP2 are insignificant because there are no resonances.

### B. Formaldehyde

The XVCC*m* and EOM-XVCC*m* results for the formaldehyde molecule are summarized in Tables VII and VIII. The calculated zero-point energy and frequencies of the fundamentals (except for *ν*_{1} and *ν*_{5}) confirm the general observations made above for the water molecule.

Method^{a}
. | E_{0}
. | ν_{1}
. | ν_{2}
. | ν_{3}
. | ν_{4}
. | ν_{5}
. | ν_{3}ν_{6}
. | ν_{2}ν_{6}
. | ν_{6}
. |
---|---|---|---|---|---|---|---|---|---|

VCI | 5813.6 | 2837.2 | 1719.8 | 1501.5 | 1160.7 | 2873.3 | 2702.8 | 2989.2 | 1237.5 |

Harmonic | 74.3 | 136.2 | 33.2 | 38.6 | 36.2 | 174.2 | … | … | 29.4 |

XVSCF(4) | 73.0 | 185.0 | 41.0 | 10.8 | −13.7 | 216.4 | … | … | −3.3 |

XVSCF[4] | 22.9 | 79.0 | 20.4 | 9.8 | 10.6 | 106.7 | … | … | 11.2 |

XVCCS{4} | 24.0 | 134.8 | 31.9 | 9.4 | −1.3 | 164.7 | … | … | 4.0 |

XVCCS(4) | 25.5 | 134.1 | 31.9 | 9.7 | −0.8 | 163.6 | … | … | 4.4 |

XVCCS[4] | 22.9 | 79.1 | 20.9 | 9.3 | 10.6 | 106.7 | … | … | 11.1 |

XVCCSD{4} | 23.5 | 30.1 | 6.7 | 5.1 | 6.7 | 38.1 | 31.6 | 34.6 | 6.4 |

XVCCSD(4) | 24.1 | 31.6 | 6.6 | 5.4 | 7.5 | 38.9 | 32.2 | 35.4 | 6.8 |

XVCCSD[4] | 22.9 | 24.2 | 6.0 | 4.2 | 6.0 | 31.0 | 35.5 | 33.0 | 5.6 |

XVCCSDT{4} | 0.6 | 0.3 | 0.5 | 1.2 | 1.7 | 1.9 | 7.3 | 2.2 | 1.4 |

XVCCSDT(4) | 0.6 | −0.0 | 0.5 | 1.2 | 1.7 | 1.8 | 7.9 | 2.4 | 1.4 |

XVCCSDT[4] | 0.8 | 0.1 | 0.5 | 1.2 | 1.9 | 1.3 | 6.1 | 1.4 | 1.4 |

XVCCSDTQ{4} | 0.0 | −0.2 | 0.0 | 0.1 | 0.1 | 0.5 | 2.4 | 0.9 | 0.1 |

XVCCSDTQ(4) | 0.0 | −0.3 | 0.0 | 0.1 | 0.1 | 0.4 | 2.4 | 0.9 | 0.1 |

XVCCSDTQ[4] | 0.1 | 0.2 | 0.0 | 0.1 | 0.0 | 0.7 | 2.3 | 1.0 | 0.0 |

XVMP2(4)^{b} | 1.3 | −190.9 | 0.7 | 0.9 | 2.7 | 33.0 | … | … | 3.0 |

XVMP2[4]^{b} | −1.0 | −15.9 | 0.1 | 1.3 | 1.6 | −414.3 | … | … | 1.8 |

XVMP2(4)^{c} | 1.3 | −6.6 | 1.4 | 1.0 | 2.6 | −0.2 | 6.6 | 29.3 | 3.0 |

XVMP2[4]^{c} | −1.0 | −8.9 | 0.3 | 1.4 | 1.5 | 2.2 | 18.0 | 24.8 | 1.8 |

Method^{a}
. | E_{0}
. | ν_{1}
. | ν_{2}
. | ν_{3}
. | ν_{4}
. | ν_{5}
. | ν_{3}ν_{6}
. | ν_{2}ν_{6}
. | ν_{6}
. |
---|---|---|---|---|---|---|---|---|---|

VCI | 5813.6 | 2837.2 | 1719.8 | 1501.5 | 1160.7 | 2873.3 | 2702.8 | 2989.2 | 1237.5 |

Harmonic | 74.3 | 136.2 | 33.2 | 38.6 | 36.2 | 174.2 | … | … | 29.4 |

XVSCF(4) | 73.0 | 185.0 | 41.0 | 10.8 | −13.7 | 216.4 | … | … | −3.3 |

XVSCF[4] | 22.9 | 79.0 | 20.4 | 9.8 | 10.6 | 106.7 | … | … | 11.2 |

XVCCS{4} | 24.0 | 134.8 | 31.9 | 9.4 | −1.3 | 164.7 | … | … | 4.0 |

XVCCS(4) | 25.5 | 134.1 | 31.9 | 9.7 | −0.8 | 163.6 | … | … | 4.4 |

XVCCS[4] | 22.9 | 79.1 | 20.9 | 9.3 | 10.6 | 106.7 | … | … | 11.1 |

XVCCSD{4} | 23.5 | 30.1 | 6.7 | 5.1 | 6.7 | 38.1 | 31.6 | 34.6 | 6.4 |

XVCCSD(4) | 24.1 | 31.6 | 6.6 | 5.4 | 7.5 | 38.9 | 32.2 | 35.4 | 6.8 |

XVCCSD[4] | 22.9 | 24.2 | 6.0 | 4.2 | 6.0 | 31.0 | 35.5 | 33.0 | 5.6 |

XVCCSDT{4} | 0.6 | 0.3 | 0.5 | 1.2 | 1.7 | 1.9 | 7.3 | 2.2 | 1.4 |

XVCCSDT(4) | 0.6 | −0.0 | 0.5 | 1.2 | 1.7 | 1.8 | 7.9 | 2.4 | 1.4 |

XVCCSDT[4] | 0.8 | 0.1 | 0.5 | 1.2 | 1.9 | 1.3 | 6.1 | 1.4 | 1.4 |

XVCCSDTQ{4} | 0.0 | −0.2 | 0.0 | 0.1 | 0.1 | 0.5 | 2.4 | 0.9 | 0.1 |

XVCCSDTQ(4) | 0.0 | −0.3 | 0.0 | 0.1 | 0.1 | 0.4 | 2.4 | 0.9 | 0.1 |

XVCCSDTQ[4] | 0.1 | 0.2 | 0.0 | 0.1 | 0.0 | 0.7 | 2.3 | 1.0 | 0.0 |

XVMP2(4)^{b} | 1.3 | −190.9 | 0.7 | 0.9 | 2.7 | 33.0 | … | … | 3.0 |

XVMP2[4]^{b} | −1.0 | −15.9 | 0.1 | 1.3 | 1.6 | −414.3 | … | … | 1.8 |

XVMP2(4)^{c} | 1.3 | −6.6 | 1.4 | 1.0 | 2.6 | −0.2 | 6.6 | 29.3 | 3.0 |

XVMP2[4]^{c} | −1.0 | −8.9 | 0.3 | 1.4 | 1.5 | 2.2 | 18.0 | 24.8 | 1.8 |

^{a}

XVCC for zero-point energy and EOM-XVCC for frequencies. Suffix “{4}” means the use of a harmonic reference, “(4)” means the use of a XVSCF(4) reference, and “[4]” means the use of a XVSCF[4] reference.

^{b}

In the frequency-independent approximation.^{59}

^{c}

Solutions of the Dyson equation with the frequency-dependent self-energy.^{59}

Method^{a}
. | E_{0}
. | ν_{1}
. | ν_{2}
. | ν_{3}
. | ν_{4}
. | ν_{5}
. | ν_{3}ν_{6}
. | ν_{2}ν_{6}
. | ν_{6}
. |
---|---|---|---|---|---|---|---|---|---|

XVCC5{4} | −0.02 | −0.03 | 0.00 | 0.01 | 0.00 | 0.03 | 0.33 | 0.06 | 0.01 |

XVCC5(4) | −0.02 | −0.06 | 0.00 | 0.01 | 0.00 | 0.01 | 0.32 | 0.05 | 0.01 |

XVCC5[4] | −0.02 | 0.02 | 0.00 | 0.01 | −0.00 | 0.04 | 0.26 | 0.06 | 0.01 |

XVCC6{4} | −0.001 | 0.000 | 0.001 | 0.001 | −0.001 | 0.005 | 0.067 | 0.011 | 0.001 |

XVCC6(4) | −0.001 | −0.002 | 0.000 | 0.001 | −0.001 | 0.003 | 0.062 | 0.010 | 0.001 |

XVCC6[4] | −0.000 | 0.002 | 0.000 | 0.000 | −0.002 | 0.006 | 0.052 | 0.008 | 0.000 |

XVCC7{4} | −0.0003 | 0.0012 | 0.0002 | 0.0004 | 0.0003 | 0.0010 | 0.0131 | 0.0021 | 0.0005 |

XVCC7(4) | −0.0002 | 0.0014 | 0.0002 | 0.0004 | 0.0003 | 0.0010 | 0.0124 | 0.0022 | 0.0005 |

XVCC7[4] | −0.0004 | −0.0011 | 0.0002 | 0.0006 | 0.0007 | −0.0005 | 0.0100 | 0.0004 | 0.0007 |

XVCC8{4} | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

XVCC8(4) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

XVCC8[4] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

Method^{a}
. | E_{0}
. | ν_{1}
. | ν_{2}
. | ν_{3}
. | ν_{4}
. | ν_{5}
. | ν_{3}ν_{6}
. | ν_{2}ν_{6}
. | ν_{6}
. |
---|---|---|---|---|---|---|---|---|---|

XVCC5{4} | −0.02 | −0.03 | 0.00 | 0.01 | 0.00 | 0.03 | 0.33 | 0.06 | 0.01 |

XVCC5(4) | −0.02 | −0.06 | 0.00 | 0.01 | 0.00 | 0.01 | 0.32 | 0.05 | 0.01 |

XVCC5[4] | −0.02 | 0.02 | 0.00 | 0.01 | −0.00 | 0.04 | 0.26 | 0.06 | 0.01 |

XVCC6{4} | −0.001 | 0.000 | 0.001 | 0.001 | −0.001 | 0.005 | 0.067 | 0.011 | 0.001 |

XVCC6(4) | −0.001 | −0.002 | 0.000 | 0.001 | −0.001 | 0.003 | 0.062 | 0.010 | 0.001 |

XVCC6[4] | −0.000 | 0.002 | 0.000 | 0.000 | −0.002 | 0.006 | 0.052 | 0.008 | 0.000 |

XVCC7{4} | −0.0003 | 0.0012 | 0.0002 | 0.0004 | 0.0003 | 0.0010 | 0.0131 | 0.0021 | 0.0005 |

XVCC7(4) | −0.0002 | 0.0014 | 0.0002 | 0.0004 | 0.0003 | 0.0010 | 0.0124 | 0.0022 | 0.0005 |

XVCC7[4] | −0.0004 | −0.0011 | 0.0002 | 0.0006 | 0.0007 | −0.0005 | 0.0100 | 0.0004 | 0.0007 |

XVCC8{4} | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

XVCC8(4) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

XVCC8[4] | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

^{a}

The *ν*_{5} fundamental is known to undergo Fermi resonances with combinations *ν*_{3}*ν*_{6} and *ν*_{2}*ν*_{6}. For this mode, XVMP2[4] in the frequency-independent approximation shows a sign of divergent perturbation series, with a large error of −414 cm^{−1} (the VMP2 result is even worse).^{59} XVMP2(4) in the same approximation, on the other hand, is reasonable (with an error of 33 cm^{−1}) but breaks down unexpectedly for nonresonant *ν*_{1} with an error of −191 cm^{−1}.^{59} Therefore, the behavior of XVMP2 in the frequency-independent approximation is rather unpredictable. This method also lacks roots for Fermi-resonant counterparts.

Solving the Dyson equation with the frequency-dependent XVMP2 self-energy rectifies the aforementioned problems of the frequency-independent approximation.^{59} With this approach, the frequency of *ν*_{5} is reproduced within 3 cm^{−1} of VCI and the roots corresponding to Fermi-resonant counterparts, *ν*_{3}*ν*_{6} and *ν*_{2}*ν*_{6}, can be located within 30 cm^{−1}. The calculated results of *ν*_{1} also regain stability with an error of less than 10 cm^{−1}. The price one pays is that a numerically cumbersome pole search needs to be performed along each mode.

Without any pole search, EOM-XVCC yields results that display stable convergence and have generally high accuracy for these strongly correlated modes. EOM-XVCCSDT seems to systematically outperform XVMP2 with the frequency-dependent self-energy for all listed modes including the Fermi resonances. EOM-XVCCSDTQ yields near-exact results with the largest error being merely 2.4 cm^{−1}. Remarkably, EOM-XVCCSD, though lacking an $R\u02c63$ operator to fully describe cubic anharmonicity, locates all three Fermi-resonant transitions (*ν*_{5}, *ν*_{3}*ν*_{6}, and *ν*_{2}*ν*_{6}) with uniform, medium accuracy. XVSCF and EOM-XVCCS, as expected, do not have roots corresponding to *ν*_{3}*ν*_{6} and *ν*_{2}*ν*_{6} and give rather poor results for *ν*_{1} and *ν*_{5}.

The combinations and overtones accessible by XVMP2 with the frequency-dependent self-energy are limited to those that are in resonance with a fundamental, which may be an advantage in some cases, but a disadvantage in others. EOM-XVCC*m* does not have such a limitation and has roots for all low-lying states involving up to *m* modes and is accurate for those involving up to (*m* − 1) modes. This is analogous to the electronic case: EOM-CCSD has roots for up to two-electron excited states and is accurate for one-electron excitation energies; for accurate results for two-electron excitation energies, one needs to use EOM-CCSDT.^{89} Table IX attests to this general trend for two- and three-mode combinations or overtones. The errors in EOM-XVCCSD for two-mode excitations and those in EOM-XVCCSDT for three-mode excitations are on the same order of magnitude as those in EOM-XVCCS for nonresonant fundamentals (*ν*_{2}, *ν*_{3}, *ν*_{4}, and *ν*_{6}). EOM-XVCCSDTQ for three-mode excitations, EOM- XVCCSDT for two-mode excitations, and EOM-XVCCSD for nonresonant fundamentals have comparable, high accuracy.

Method . | (000200) . | (000101) . | (000002) . | (001100) . | (010100) . | (002000) . | (000300) . | (000201) . |
---|---|---|---|---|---|---|---|---|

VCI | 2318.1 | 2399.6 | 2472.5 | 2660.6 | 2871.2 | 3000.8 | 3464.2 | 3555.2 |

EOM-XVCCSD[4] | 35.7 | 28.3 | 31.0 | 77.7 | 38.3 | 23.2 | … | … |

EOM-XVCCSDT[4] | 7.6 | 6.0 | 5.3 | 5.7 | 4.5 | 4.7 | 62.5 | 49.6 |

EOM-XVCCSDTQ[4] | 2.1 | 2.0 | 2.1 | 2.0 | 1.8 | 2.0 | 14.2 | 10.9 |

EOM-XVCC5[4] | 0.1 | 0.1 | 0.2 | 0.1 | 0.1 | 0.2 | 3.3 | 3.0 |

EOM-XVCC6[4] | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.5 | 0.3 |

EOM-XVCC7[4] | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.1 | 0.1 |

EOM-XVCC8[4] | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |

Method . | (000200) . | (000101) . | (000002) . | (001100) . | (010100) . | (002000) . | (000300) . | (000201) . |
---|---|---|---|---|---|---|---|---|

VCI | 2318.1 | 2399.6 | 2472.5 | 2660.6 | 2871.2 | 3000.8 | 3464.2 | 3555.2 |

EOM-XVCCSD[4] | 35.7 | 28.3 | 31.0 | 77.7 | 38.3 | 23.2 | … | … |

EOM-XVCCSDT[4] | 7.6 | 6.0 | 5.3 | 5.7 | 4.5 | 4.7 | 62.5 | 49.6 |

EOM-XVCCSDTQ[4] | 2.1 | 2.0 | 2.1 | 2.0 | 1.8 | 2.0 | 14.2 | 10.9 |

EOM-XVCC5[4] | 0.1 | 0.1 | 0.2 | 0.1 | 0.1 | 0.2 | 3.3 | 3.0 |

EOM-XVCC6[4] | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.5 | 0.3 |

EOM-XVCC7[4] | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.1 | 0.1 |

EOM-XVCC8[4] | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |

For transition frequencies to strongly correlated states involving at least two modes (Fermi resonances, combinations, and overtones), the following is the order of accuracy: XVCCSD < XVMP2 with pole search < XVCCSDT < XVCCSDTQ < XVCC5 < XVCC6 < XVCC7 ≤ XVCC8. The methods not mentioned in the list either do not have the corresponding roots or tend to display signs of nonconvergence.

## V. CONCLUSION

We have presented an accurate, stable, efficient, and widely applicable vibrational CC (XVCC) theory for both zero-point energies and transition frequencies, the latter through the EOM formalism, of a molecule with an anharmonic PES. It relies on a Hartree product of harmonic-oscillator modals as a reference wave function and a truncated Taylor-series PES. They lead to streamlined sets of equations defining XVCC and EOM-XVCC, which mirror exactly those of electronic CC and EOM-CC theories. XVCC and EOM-XVCC are diagrammatically size consistent, free of modal basis sets, and thus an inherently infinite-basis theory.

We have established the rules for enumerating Feynman– Goldstone diagrams for XVCC and EOM-XVCC theories and the ones for translating them into algebraic equations. This method of derivation is not only expedient but also insightful because the diagrammatic topology can more immediately determine size consistency of a theory or the lack thereof as well as the relationships among various diagrammatic theories. The rules are justified by the underlying normal-ordered second-quantized method of derivation, which gives the same equations as the diagrammatic derivation in each case.

XVCC and EOM-XVCC theories are general in the sense that their formalisms are unchanged for any single-reference Hartree product of harmonic-oscillator modals, including the wave functions of XVSCF(*n*), XVSCF[*n*], SCP, or the harmonic approximation. Reflecting the centrality of XVSCF (which includes SCP) as the size-consistent single-Hartree-product method in many-body vibrational theory, some diagrammatic contributions in XVCC and EOM-XVCC vanish in this reference, simplifying their equations further. Low-order XVCC and EOM-XVCC results are also the same or very close to the XVSCF results for the same reason.

The XVCC*m* and EOM-XVCC*m* methods up to *m* = 8 have been formulated and implemented by computer algebra of the aforementioned diagram logic as well as of algorithmic optimizations such as strength reduction, intermediate reuse, and factorization. The numerical tests of these methods show that the results are smoothly and rapidly convergent towards the exact values. This is in contrast to XVMP2 in the frequency-independent approximation to the self-energy, which displays signs of divergent perturbation series not only for resonant modes but also for nonresonant modes. While this problem is rectified by solving the Dyson equation with the frequency-dependent self-energy in XVMP2, EOM-XVCC may be said to achieve the same with a more well-studied diagonalization algorithm instead of a pole search. XVCC*m* and EOM-XVCC*m* with *m* = *n*, where *n* is the truncation rank of the Taylor-series PES, seem to be ones with the best cost-accuracy balance and hence the greatest practical utility.

The VCC theory of Prasad and coworkers^{27,32} can be identified as XVCC and EOM-XVCC with the XVSCF[*n*] reference. Our work is, therefore, built upon the previous studies of Prasad and coworkers,^{27,32} but it generalizes the theories to a variety of reference wave functions, establishes the normal-ordered, diagrammatic, and automated derivations, and implements higher-order methods up to connected octuple excitations using computerized code optimization and synthesis. Christiansen’s VCC theory,^{36} on the other hand, relies on VSCF and modals that are numerically defined and differs from XVCC; unlike XVCC, VCC likely retains many terms that are not size consistent and is expected to be fundamentally less efficient for large molecules and solids despite the exponential parameterization of its wave function. In VCC, the rank of disconnected excitations is also limited by the number of basis functions, which is not the case with basis-set-free XVCC.

VSCF is exact for any one-dimensional potential in the infinite basis-set limit, whereas XVSCF captures only a small portion of its anharmonic effect. For solids, VSCF reduces to XVSCF, the former giving the same results as the latter, but at a much greater computational cost. This suggests that Christiansen’s VCC may also perform well for small molecules with particularly strong anharmonicity, while reducing to more streamlined XVCC or the VCC theory of Prasad and coworkers in the thermodynamic limit. Our numerical results show that, for Morse-like anharmonicity, overtones and combinations, and anharmonic resonances, XVCC can nevertheless rectify the deficiency of XVSCF for small molecules. However, XVSCF and XVCC are expected to not work well for a double-well potential, for instance.

Furthermore, rectilinear coordinates such as normal coordinates, also underlying XVCC, may become awkward for large-amplitude vibrations. The methods that are suitable for such *strong correlation* and handle different coordinates flexibly include the multiconfiguration time-dependent Hartree (MCTDH) method,^{90} which is variational, but may lack size consistency.^{91} XVMP2 and XVCC are designed primarily for *weak correlation* in the ground state and are nonvariational and diagrammatically size consistent.^{91} Excited states are usually strongly correlated but can be handled by either class of methods.^{59} Therefore, it may be said that XVCC is not well suited for strong correlation in the ground state, but it is an excellent—arguably the best—size-consistent method for weak correlation in the ground state and for strong and weak correlation in excited states.

## Acknowledgments

This work is supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Grant No. DE-FG02-11ER16211. J.A.F. is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1144245. The authors thank Matthew R. Hermes for useful discussions and supplying the QFF’s and VCI results used in this work.

### APPENDIX A: RIGOROUS DERIVATION OF THE CONNECTED FORM OF THE XVCC EQUATIONS

As before, we demand the XVCC wave function, $eT\u02c6|\Phi 0\u3009$, satisfy the vibrational Schrödinger equation,

but only in the space of Hartree products reachable by $T\u02c6$ from |Φ_{0}〉. Prior to projecting this equation onto this space, we first act $e\u2212T\u02c6$ on both sides of the equation,

The projection of this equation then yields the energy and amplitude equations, which read

for all *ν*_{1}, …, *ν _{n}* and each

*n*(1 ≤

*n*≤

*m*). It may be understood that the projection space 〈Φ

_{0}|∪{〈Φ

_{ν1ν2⋯νn}|} used in the derivation in the main text spans the same space $\u3008\Phi 0|e\u2212T\u02c6\u222a{\u3008\Phi \nu 1\nu 2\cdots \nu n|e\u2212T\u02c6}$ as above and, therefore, the two derivations lead to one and the same theory.

The similarity-transformed Hamiltonian, $e\u2212T\u02c6H\u02c6NeT\u02c6$, can be rewritten using the Baker–Campbell–Hausdorff expansion^{1,2} as

which terminates with the *n*-fold commutator, where *n* is the truncation rank of the Taylor-series PES. The commutators immediately guarantee the diagrammatically connected structure of each surviving term. Therefore, we may write

### APPENDIX B: XVCCSDT

Here, we document the algebraic and diagrammatic forms of the XVCCSDT energy and amplitude equations in a QFF. The diagrams in the energy equation that are not already listed in the XVCCSD energy equation are shown in Fig. 16, and the corresponding algebraic expressions are given in Table X. As the numerical results illustrate (see above), the inclusion of connected triple excitations is critical for highly accurate zero-point energies. This is attributed to the importance of cubic force constants and diagram *E*_{9}. The *T*_{1} and *T*_{2} amplitude equations are given by

where $T1(2)$ and $T2(2)$ are defined by Eqs. (51) and (53), respectively. The additional diagrams in the *T*_{1} and *T*_{2} amplitude equations are shown in Figs. 17 and 18. The algebraic interpretation of the diagrams are given in Tables XI and XII. The *T*_{3} amplitude equation is written as

The diagrammatic form of this equation is shown in Fig. 19 with the corresponding algebraic interpretation in Table XIII.

$\Delta E0(3)=\Delta E0(2)+13!Wijk\tau ijk+13!Wijkl\tau i\tau jkl$ |

$\Delta E0(3)=\Delta E0(2)+13!Wijk\tau ijk+13!Wijkl\tau i\tau jkl$ |

$T1(3)\u2009\u2261\u2009T1(2)+12!Wij\u2212\tau ij\nu +13!W\nu ijk\tau ijk+12!Wijk\tau i\tau jk\nu $ |

$+12!2Wijkl\tau ij\tau kl\nu +13!Wijkl\tau i\nu \tau jkl+12!2Wijkl\tau i\tau j\tau kl\nu $ |

= 0 |

$T1(3)\u2009\u2261\u2009T1(2)+12!Wij\u2212\tau ij\nu +13!W\nu ijk\tau ijk+12!Wijk\tau i\tau jk\nu $ |

$+12!2Wijkl\tau ij\tau kl\nu +13!Wijkl\tau i\nu \tau jkl+12!2Wijkl\tau i\tau j\tau kl\nu $ |

= 0 |

$T2(3)$ | ≡ | $T2(2)+P\u02c6(\nu \mu )(12!Wi\tau i\nu \mu +12!W\nu ij\tau ij\mu +12!Wij\u2212\tau i\tau j\nu \mu $ |

$+12!W\nu ijk\tau i\tau jk\mu +12!Wijk\tau i\nu \tau jk\mu +12!2Wijk\tau ij\tau k\nu \mu $ | ||

$+12!2Wijk\tau i\tau j\tau k\nu \mu +12!3!Wijkl\tau ijk\tau l\nu \mu +12!3Wijkl\tau ij\nu \tau kl\mu $ | ||

$+12!2Wijkl\tau i\tau jk\tau l\nu \mu +12!Wijkl\tau i\tau j\nu \tau kl\mu +12!3!Wijkl\tau i\tau j\tau k\tau l\nu \mu )$ | ||

= | 0 |

$T2(3)$ | ≡ | $T2(2)+P\u02c6(\nu \mu )(12!Wi\tau i\nu \mu +12!W\nu ij\tau ij\mu +12!Wij\u2212\tau i\tau j\nu \mu $ |

$+12!W\nu ijk\tau i\tau jk\mu +12!Wijk\tau i\nu \tau jk\mu +12!2Wijk\tau ij\tau k\nu \mu $ | ||

$+12!2Wijk\tau i\tau j\tau k\nu \mu +12!3!Wijkl\tau ijk\tau l\nu \mu +12!3Wijkl\tau ij\nu \tau kl\mu $ | ||

$+12!2Wijkl\tau i\tau jk\tau l\nu \mu +12!Wijkl\tau i\tau j\nu \tau kl\mu +12!3!Wijkl\tau i\tau j\tau k\tau l\nu \mu )$ | ||

= | 0 |

$T3(3)$ | ≡ | $P\u02c6(\nu \mu \lambda )(13!W\nu \mu \lambda +13!W\nu \mu \lambda i\tau i+12!W\nu \mu i\tau i\lambda +12!W\nu i+\tau i\mu \lambda $ |

$+12!2W\nu \mu ij\tau ij\lambda +12!W\nu \mu ij\tau i\tau j\lambda +12!W\nu ij\tau i\mu \tau j\lambda +12!W\nu ij\tau i\tau j\mu \lambda $ | ||

$+12!Wij\u2212\tau i\nu \tau j\mu \lambda +12!W\nu ijk\tau i\mu \tau jk\lambda +12!2W\nu ijk\tau ij\tau k\mu \lambda $ | ||

$+12!2Wijk\tau ij\nu \tau k\mu \lambda +12!W\nu ijk\tau i\tau j\mu \tau k\lambda +13!Wijk\tau i\nu \tau j\mu \tau k\lambda $ | ||

$+12!2W\nu ijk\tau i\tau j\tau k\mu \lambda +12!Wijk\tau i\tau j\nu \tau k\mu \lambda +12!2Wijkl\tau i\tau jk\nu \tau l\mu \lambda $ | ||

$+12!2Wijkl\tau i\nu \tau j\mu \tau kl\lambda +12!2Wijkl\tau ij\tau k\nu \tau l\mu \lambda +13!Wijkl\tau i\tau j\nu \tau k\mu \tau l\lambda $ | ||

$+12!2Wijkl\tau i\tau j\tau k\nu \tau l\mu \lambda )$ | ||

= | 0 |

$T3(3)$ | ≡ | $P\u02c6(\nu \mu \lambda )(13!W\nu \mu \lambda +13!W\nu \mu \lambda i\tau i+12!W\nu \mu i\tau i\lambda +12!W\nu i+\tau i\mu \lambda $ |

$+12!2W\nu \mu ij\tau ij\lambda +12!W\nu \mu ij\tau i\tau j\lambda +12!W\nu ij\tau i\mu \tau j\lambda +12!W\nu ij\tau i\tau j\mu \lambda $ | ||

$+12!Wij\u2212\tau i\nu \tau j\mu \lambda +12!W\nu ijk\tau i\mu \tau jk\lambda +12!2W\nu ijk\tau ij\tau k\mu \lambda $ | ||

$+12!2Wijk\tau ij\nu \tau k\mu \lambda +12!W\nu ijk\tau i\tau j\mu \tau k\lambda +13!Wijk\tau i\nu \tau j\mu \tau k\lambda $ | ||

$+12!2W\nu ijk\tau i\tau j\tau k\mu \lambda +12!Wijk\tau i\tau j\nu \tau k\mu \lambda +12!2Wijkl\tau i\tau jk\nu \tau l\mu \lambda $ | ||

$+12!2Wijkl\tau i\nu \tau j\mu \tau kl\lambda +12!2Wijkl\tau ij\tau k\nu \tau l\mu \lambda +13!Wijkl\tau i\tau j\nu \tau k\mu \tau l\lambda $ | ||

$+12!2Wijkl\tau i\tau j\tau k\nu \tau l\mu \lambda )$ | ||

= | 0 |

### APPENDIX C: EOM-XVCCSD

The *R*_{1} and *R*_{2} amplitude equations of EOM-XVCCSD in a QFF are written as

where $R1(1)$ is defined by Eq. (70). The additional diagrams in the *R*_{1} amplitude equation are depicted in Fig. 20. Their corresponding algebraic interpretation is given in Table XIV. The diagrammatic form of the *R*_{2} amplitude equation is drawn in Fig. 21, while the algebraic form is given in Table XV. EOM-XVCCSD provides the first sizable correction to the vibrational frequencies. This is attributed to the inclusion of diagrams such as $S5\u2032$, which account for the important cubic force constants that are not multiplied by near-zero *T*_{2} amplitudes. In particular, using a XVSCF[4] reference, EOM-XVCCS yields frequencies that are equivalent to the eigenvalues of the $Wij+$ matrix, completely neglecting cubic force constants. This is rectified by EOM-XVCCSD giving results with accuracy comparable to XVMP2.

$R1(2)$ | ≡ | $R1(1)+Wiri\nu +12!W\nu ijrij+Wij\u2212\tau j\nu ri+Wij\u2212\tau jri\nu $ |

$+12!W\nu ijk\tau jkri+12!W\nu ijk\tau krij+12!Wijk\tau jkri\nu +12!Wijk\tau k\nu rij$ | ||

$+Wijk\tau j\tau k\nu ri+12!Wijk\tau j\tau kri\nu +12!Wijkl\tau jk\tau l\nu ri$ | ||

$+12!Wijkl\tau k\tau l\nu rij+12!Wijkl\tau j\tau klri\nu +12!Wijkl\tau j\tau k\tau l\nu ri$ | ||

$+13!Wijkl\tau j\tau k\tau lri\nu $ | ||

= | $\omega k(2)r\nu $ |

$R1(2)$ | ≡ | $R1(1)+Wiri\nu +12!W\nu ijrij+Wij\u2212\tau j\nu ri+Wij\u2212\tau jri\nu $ |

$+12!W\nu ijk\tau jkri+12!W\nu ijk\tau krij+12!Wijk\tau jkri\nu +12!Wijk\tau k\nu rij$ | ||

$+Wijk\tau j\tau k\nu ri+12!Wijk\tau j\tau kri\nu +12!Wijkl\tau jk\tau l\nu ri$ | ||

$+12!Wijkl\tau k\tau l\nu rij+12!Wijkl\tau j\tau klri\nu +12!Wijkl\tau j\tau k\tau l\nu ri$ | ||

$+13!Wijkl\tau j\tau k\tau lri\nu $ | ||

= | $\omega k(2)r\nu $ |

$R2(2)$ | ≡ | $P\u02c6(\nu \mu )(12!W\nu \mu iri+W\nu i+ri\mu +12!2W\nu \mu ijrij+12!W\nu \mu ij\tau jri$ |

$+W\nu ij\tau j\mu ri+W\nu ij\tau jri\mu +Wij\u2212\tau j\mu ri\nu +12!W\nu ijk\tau jkri\mu $ | ||

$+12!W\nu ijk\tau k\mu rij+W\nu ijk\tau j\tau k\mu ri+12!Wijk\tau j\nu \tau k\mu ri$ | ||

$+12!W\nu ijk\tau j\tau kri\mu +Wijk\tau j\tau k\mu ri\nu +12!Wijkl\tau jk\tau l\mu ri\nu $ | ||

$+12!2Wijkl\tau k\nu \tau l\mu rij+12!Wijkl\tau j\tau k\nu \tau l\mu ri+12!Wijkl\tau j\tau k\tau l\mu ri\nu )$ | ||

= | $\omega k(2)r\nu \mu $ |

$R2(2)$ | ≡ | $P\u02c6(\nu \mu )(12!W\nu \mu iri+W\nu i+ri\mu +12!2W\nu \mu ijrij+12!W\nu \mu ij\tau jri$ |

$+W\nu ij\tau j\mu ri+W\nu ij\tau jri\mu +Wij\u2212\tau j\mu ri\nu +12!W\nu ijk\tau jkri\mu $ | ||

$+12!W\nu ijk\tau k\mu rij+W\nu ijk\tau j\tau k\mu ri+12!Wijk\tau j\nu \tau k\mu ri$ | ||

$+12!W\nu ijk\tau j\tau kri\mu +Wijk\tau j\tau k\mu ri\nu +12!Wijkl\tau jk\tau l\mu ri\nu $ | ||

$+12!2Wijkl\tau k\nu \tau l\mu rij+12!Wijkl\tau j\tau k\nu \tau l\mu ri+12!Wijkl\tau j\tau k\tau l\mu ri\nu )$ | ||

= | $\omega k(2)r\nu \mu $ |

### APPENDIX D: COST SCALING OF XVCC*m* AND EOM-XVCC*m* IN AN *n*TH-ORDER PES

The peak operation cost of XVCC*m* or EOM-XVCC*m* in an *n*th-order PES can be expressed as *O*(*N ^{C}*), where

*N*is the number of vibrational degrees of freedom and

*C*is an integer. The lowest value of

*C*(

*C*

_{min}) is

*m*+ ⌊

*n*/2⌋ for

*m*≥

*n*and is not greater than

*n*+ ⌊(

*m*− 1)/2⌋ for

*m*<

*n*.

The proof is given below separately for *m* ≥ *n* and for *m* < *n*. In both cases, we first show that the aforementioned values of *C* are an upper bound; for any tensor contraction appearing in XVCC*m* or EOM-XVCC*m*, there is at least one contraction order whose cost function rank does not exceed the respective value. For *m* ≥ *n*, we additionally show that the aforementioned value of *C* is also a lower bound; there is always one tensor contraction whose minimal cost function rank is exactly this value. Together they prove that the expression of *C* given above is minimal for *m* ≥ *n*.

In what follows, we call an *R* amplitude of EOM-XVCC a *T* amplitude because they need not be distinguished in this context. We also only need to consider sequential contraction orders in which the force-constant tensor is always contracted first without losing rigor of the proof. Let *T*^{(k)} be the *k*th *T* amplitude to be contracted, having $Text(k)$ external lines and $Tint(k)$ internal lines. Let *I*^{(k)} be the *k*th intermediate. The symbol *I*^{(k)} also denotes its size, i.e., the total number of lines. *I*^{(0)} is the force-constant tensor, i.e., *I*^{(0)} = *n*.

In the contraction of *I*^{(k−1)} with *T*^{(k)}, the number of summation indexes is $Tint(k)$, because two *T* amplitudes never share common internal lines and, therefore, all internal lines of *T*^{(k)} must connect with *I*^{(k−1)}. The size of the product, i.e., the *k*th intermediate, is $I(k\u22121)+Text(k)\u2212Tint(k)$ and the cost function rank of the contraction is $I(k\u22121)+Text(k)$.

The theorem applies not only to XVCC*m* and EOM-XVCC*m* but also to other CC and EOM-CC methods for electrons and nucleons.

Let *h* be the total number of *T* amplitudes in a tensor contraction. Given intermediate *I*^{(i)} (0 ≤ *i* ≤ *h*), there is at least one contraction order of remaining (*h* − *i*) *T* amplitudes in which none of the subsequent intermediates has size exceeding max(*I*^{(i)}, *I*^{(h)}).

*k*th intermediate is given by

*T*

^{(k)}amplitudes be ordered as follows:

*T*

^{(k)}’s with negative values of the difference in the external and internal lines ($Text(k)\u2212Tint(k)$) are contracted first, next by

*T*

^{(k)}’s with zero difference, and finally by

*T*

^{(k)}’s with positive values of the difference. Then,

*I*

^{(k)}as a function of

*k*for this contraction order is a convex function and thus cannot exceed the greater of their values at the end points,

*I*

^{(k)}≤ max(

*I*

^{(i)},

*I*

^{(h)}).

*m*≥

*n*.

First, we show that *C* = *m* + ⌊*n*/2⌋ is an upper bound. We consider the following two cases: (1) there is at least one *T*^{(k)} with $Tint(k)\u2265\u2308n/2\u2309$; (2) otherwise.

In case 1, the *T* amplitude with $Tint(1)\u2265\u2308n/2\u2309$ is contracted first with the force-constant tensor. The cost function rank of this contraction is given by

where we use $Text(1)+Tint(1)\u2264m$ in the first inequality. The size of the resulting intermediate is

where Eq. (D2) is used in the first inequality.

Next, we consider case 2 and the second and subsequent contractions in case 1 together. In the former, $Tint(k)\u2264\u230an/2\u230b$ for all *k*. In the latter, $Tint(k)\u2264\u230an/2\u230b$ for all *k* ≥ 2. Here, we contract *T* amplitudes in one of the orders of the lemma. As per the lemma, this guarantees *I*^{(k)} ≤ *m* in both cases [see Eq. (D3) in case 1]. The sizes of the *k*th and (*k* − 1)th intermediates are related to each other by

Therefore, the cost function rank of the *k*th contraction is

Second, we show that the same value of *C* is a lower bound. There is always a binary product of the force-constant tensor, *I*^{(0)}, and a *T _{m}* amplitude with ⌈

*n*/2⌉ internal lines. Note that ⌈

*n*/2⌉ cannot exceed

*n*or

*m*because

*m*≥

*n*≥ 1. Furthermore, the number of external lines is

*m*+

*n*− 2⌈

*n*/2⌉ ≤

*m*. Hence, this is a valid tensor contraction. The cost function rank of this contraction is

which proves *C*_{min} ≥ *m* + ⌊*n*/2⌋.

Together, we obtain *C*_{min} = *m* + ⌊*n*/2⌋ for *m* ≥ *n*.

This is an example of the optimization problem in computer science known as the matrix chain multiplication.^{84,92,93} In this case (*m* ≥ *n*), we have shown that there is an analytical solution that minimizes the exponent of the leading-order cost function and its value.

*m*<

*n*.

We contract *T* amplitudes in one of the orders of the lemma. As per the lemma, this ensures *I*^{(k)} ≤ *n*. We consider the following three cases in each contraction: (1) $Text(k)<Tint(k)$, (2) $Text(k)>Tint(k)$, and (3) $Text(k)=Tint(k)$.

In case 1, since $Text(k)+Tint(k)\u2264m$, we have $Text(k)\u2264\u230a(m\u22121)/2\u230b$. The cost function rank of the *k*th contraction is

In case 2, similarly, $Tint(k)\u2264\u230a(m\u22121)/2\u230b$ holds. The size of the *k*th intermediate is given by

which implies

In case 3, *m* must be even and $Text(k)=Tint(k)=m/2$. The intermediate resulting from the contraction of this *T* amplitude has the size,

The first inequality follows from the fact that *I*^{(l)} is a convex function with its minimum occurring at *l* = *k*, where $Text(k)\u2212Tint(k)=0$. The cost function rank of this step is

Equations (D7), (D9), and (D11) prove *C*_{min} ≤ *n* + ⌊(*m* − 1)/2⌋ for *m* < *n*.

This upper-bound expression for *m* < *n* happens to agree with the optimal cost function ranks for all cases of *m* < *n* in Table IV. However, we have identified many other combinations of *m* < *n* whose optimal ranks are lower.