Convergence of coupled cluster perturbation theory

In Møller-Plesset (MP) perturbation theory,^1–3 the Hartree-Fock (HF) one-determinant wave function is perturbatively corrected for the effect of single (S), double (D), etc., excitations out of the state (up to the level of N excitations, where N is the number of electrons in the system at hand). The target energy for the MP series is thus that of the full configuration interaction (FCI) model. The electronic Hamiltonian is in MP theory partitioned into the Fock operator, $f^$ (a zeroth-order part), and the fluctuation potential, $Φ^$ (the perturbation), and at lowest (second) order, the resulting MP2 model marks a successful approach to the correlation problem in quantum chemistry, providing a useful, size-extensive doubles correction to the HF energy at a low computational cost. The higher-order corrections of the MP series, however, represent less successful treatments of the electron correlation problem, in part due to the fact that these models offer significantly worse compromises between cost and accuracy than does the MP2 model, and in part because of convergence problems that reflect inherent problems with the actual partitioning of the Hamiltonian. We will now briefly revisit the reason why the MP series might occasionally diverge.

Made possible by an efficient FCI implementation,^4,5 a study appearing some twenty years ago numerically showed how, upon augmenting standard correlation-consistent basis sets of double-$ζ$ quality by diffuse functions,⁶ the MP series was at the risk of becoming oscillating and, ultimately, divergent for a few prototypical, single-reference dominated examples.⁷ These results were in stark contrast to the behavior of the series within non-augmented standard basis sets, for which the convergence behavior of the MP series had previously been found to reflect the dominance of the reference state in the exact solution^8,9 (at least for MP expansions formulated upon restricted HF (RHF).^10,11 Thus, since the convergence behavior of the MP series was found to exhibit an extreme dependence on the choice of one-electron basis set, even for simple examples dominated by a single determinant, for which no a priori reason for expecting divergences seemed to exist, the usefulness of the higher-order models of the MP series was immediately brought into question, for total energies as well as properties derived from these.¹² 

The numerical findings mentioned above were subsequently complemented by an explicit determination of the radius of convergence for the MP series.¹³ In that study, the electronic Schrödinger equation was solved for an electronic Hamiltonian expressed as $H^(z)=f^+zΦ^$⁠, where z is a complex perturbational strength parameter, such that the zeroth-order (HF) and physical (FCI) wave functions were recovered by the values z = 0 and z = 1, respectively. When regarded in this way, the MP series is recognized as a power series expansion in z of the FCI eigenvalues. This expansion will then be divergent if the ground state becomes degenerate with some excited state at a point of degeneracy, $z=ξ$⁠, within the complex unit circle defined by |z| = 1.^14,15 In practice, the authors of Ref. 13 replaced the above search for branch points of the ground state energy function by a search for avoided crossings between the ground state and excited state on the real axis only (i.e., they searched for real-valued points of avoided crossings within the interval $z∈[−1;1]$⁠), since such a search is far less complex than a search within the entire unit circle while still enabling the identification of degeneracies. In this terminology, an excited state is denoted a back-door intruder if an avoided crossing is observed in the interval from $−1<z<0$⁠, while it is denoted a front-door intruder if the crossing falls within the interval from $0<z<1$⁠. A few years later, the analysis in Ref. 13 was supplemented by a simple two-state model capable of explaining the divergence of the MP series.¹⁶ In that study, it was illustrated how divergences are bound to occur whenever a basis set is sufficiently flexible to give a reasonable description of highly excited and diffuse back-door intruders, which couple only weakly to the ground state (and, as such, are nearly invisible in the energy spectrum at either z = 0 ∨ z = 1), leading the authors to state that divergence is the rule rather than exception in MP theory, which converges only in small basis sets. This point has since then been repeatedly stressed and extended in related studies.^17–23 

We have recently presented a range of perturbation series^24–26 that are all based on a coupled cluster²⁷ (CC) wave function rather than an HF wave function as the zeroth order state, while maintaining the MP partitioning of the Hamiltonian. In the present study, we will extend the type of investigation outlined for the MP series above to these perturbative CC series as well as truncated MP series. Our investigation is centered around the bivariational CC perturbation series introduced in Ref. 24, while a few results for the recently introduced energy-based CC perturbation series²⁶ will also be presented for comparison. For the sake of brevity, we will introduce a general notation, CC[m], for the CC model resulting from a truncation of the cluster operator at a specific level, m (i.e., CC[2] is the CC singles and doubles²⁸ (CCSD) model, CC[3] is the CC singles, doubles, and triples²⁹ (CCSDT) model, etc.). In particular, we want to compare the convergence radii of (i) a family of so-called HF([$mQ$]–n) expansions, in which the correlation energy of a target CC[$mQ$] model is expanded in the fluctuation potential, with (ii) the bivariational CC[$mP$]([$mQ$]–n) expansions,²⁴ which describe the energetic difference between a parent model CC[$mP$] and a target model CC[$mQ$], again in orders of the fluctuation potential, cf. Figure 1. Thus, the HF([$mQ$]–n) series are MP-like expansions, which theoretically converge from the HF energy, not towards the exact FCI energy, but rather towards the energy of a CC target model, e.g., the target energy of the HF(SDTQ–n) series is the CCSDTQ³⁰ (CC singles, doubles, triples, and quadruples) energy. These series are thus equivalent to the standard MP series whenever the latter is restricted to at most $[mQ]$ excitations. As examples of the second type, the CCSD(T–n), CCSD(TQ–n), and CCSDT(Q–n) series have all recently been proposed,^31–34 forming order expansions that converge from either the CCSD or CCSDT energies onto the CCSDT or CCSDTQ energies, subject to the same partitioning of the electronic Hamiltonian as that used in MP theory, i.e., $H^=f^+Φ^$⁠. Through lowest (second and third) orders, these series encompass (or are related to) a number of established perturbation models developed within different theoretical frameworks;^35–38 we here refer the reader to Refs. 24, 31, and 33 for recent reviews and comparisons of the various different triples and quadruples models that exist in the literature.

For both families of expansions (based on either an HF or a CC reference state), we will here probe for possible back- and front-door intruders in the real-valued interval $z∈[−1;1]$ by searching for those z-values at which the ground state becomes degenerate with an excited intruder state. We consider four examples, all of which have previously been shown to have slowly convergent or divergent MP series—the Ne atom, the F^– anion, and the singlet CH₂ and HF molecules at equilibrium and distorted geometries, respectively.^7,13,16 Furthermore, numerical results for the CC-based perturbation series for different choices of parent and target models will be reported in order to assess to what degree potential divergences might influence these. In this respect, we wish to distinguish between two types of convergence, namely, (i) formal convergence, i.e., does a given expansion indeed converge, and (ii) practical convergence, i.e., are the lowest-order corrections of a given series physically meaningful? Thus, while the formal convergence may be assessed by explicitly determining the radius of convergence, the practical convergence may be assessed simply by inspecting the lowest-order energy corrections. For this purpose, we have implemented arbitrary-order energy corrections for all members of the CC[$mP$]([$mQ$]–n) family of CC perturbation series.

We will now outline the theory behind the intruder state scan and the computation of CC[$mP$]([$mQ$]–n) energy corrections to arbitrary order. In Section II A, we discuss how the radii of convergence for the perturbation expansions are defined by the largest perturbation parameter for which the CC Jacobian remains non-singular. In Section II B, we use this insight to scrutinize the theoretical origin of possible intruder states from a comparison of the Jacobians for HF- and CC-based perturbation series. Finally, in Section II C, we review our implementation of general-order CC[$mP$]([$mQ$]–n) energy corrections.

The original analysis of the convergence of MP perturbation theory discussed in Section I may not directly be applied to the present CC context, as the energies and wave function parameters are now obtained from a non-linear set of equations, rather than as the solutions to a standard Hermitian eigenvalue problem for the Hamiltonian operator,

In CC theory, the energy is in general an analytic algebraic function of the cluster parameters, {t}, and the perturbation strength, z, i.e., $ECC=ECC(z,𝐭)$⁠. In searching for intruder states, the task boils down to determining the range of z within which the energy may be expressed in terms of an analytic function of z, or, equivalently, as a convergent Taylor series in z. In turn, this range of z is identical to the one where the cluster amplitudes may be formulated as analytic functions of z. In order to actually determine this range, we write the non-linear set of equations that define the cluster amplitudes on the generic form

with one equation for each individual cluster amplitude. Assuming that for z = z₀, the amplitudes {t₀} satisfy the equations v(z₀,t₀) = 0, the implicit function theorem³⁹ will ensure that the amplitudes are indeed analytic functions of z in a neighborhood of z₀, provided that the CC Jacobian, i.e., the derivative of Eq. (2.2) with respect to {t}, is non-singular in the point (z₀, t₀). By examining the Jacobian as a function of z, we may determine a value $ξ$ as that value of z that has the smallest norm while still giving rise to a singular Jacobian. For $z<|ξ|$⁠, the cluster amplitudes—and hence the energy—is thus analytic, implying that $ECC(z,𝐭)$ may be expanded as a convergent Taylor expansion around z = 0.

In standard CC response theory, the eigenvalues of the Jacobian are identified as excitation energies, which implicates that a singularity in the Jacobian corresponds to a degeneracy between the ground state and an excited state. Thus, although the convergence criteria for the standard MP series and the present CC perturbation series are obtained in different manners, the radius of convergence is in both cases identified as the lowest value of |z| at which such a degeneracy occurs. Phrased differently, the determination of the radius of convergence requires the determination of a complex variable z that produces a zero-eigenvalue of the Jacobian. As discussed in Ref. 16, the imaginary part of this branch point is in general small for back-door intruders, whereas it may be large for front-door intruders. In both cases, its approximate location may be determined from a scan over the eigenvalues of the Jacobian for real values of z, and on par with the procedure in Refs. 13 and 16, we will therefore only examine the lowest eigenvalue of the Jacobian as a function of real values of z. Thus, instead of exact degeneracies, we will probe for avoided crossings where the real value of the lowest Jacobian eigenvalue reaches a minimum. For a back-door intruder, the smallness of the imaginary part of the branch point will often lead to very pronounced avoided crossings with the lowest eigenvalue approaching zero, whereas for front-door intruders, the large imaginary part of the branch point tends to produce less pronounced avoided crossings.

Generally speaking, the CC Jacobian is a non-symmetric matrix. This fact leads to two main differences compared to the standard approach of FCI, for which eigenvalues of a Hermitian matrix are determined (assuming a vanishing imaginary component, $Im(z)=0$⁠). First, a non-symmetric matrix need not have a complete set of eigenvalues and eigenvectors. Thus, whereas the singularity of the Jacobian matrix at the branch point trivially corresponds to a zero eigenvalue, it cannot in general be ensured that this eigenvalue changes into another very small eigenvalue for a value of z infinitesimally close to the branch point. However, we note that such abrupt “dissolutions” of zero eigenvalues were never encountered in the course of the present study. Second, the eigenvalues of a non-symmetric matrix may be complex. Thus, whereas, for the Hermitian case, the occurrence of an intruder state lower in energy than the reference state is a clear sign of a singularity (and hence a branch point), this may not be the case for a non-Hermitian problem. However, since complex eigenvalues of a real matrix always occur in pairs with identical real parts and opposite imaginary parts, and due to the analyticity of the eigenvalues of a general non-degenerate matrix as functions of its elements, a real eigenvalue may only turn into a complex eigenvalue at a branch point. The immediate consequence of this is the practical convenience of being able to determine branch points in much the same way for CC perturbation theory as was previously realized for MP theory in Refs. 13 and 16.

In formulating either of the HF- or CC-based series discussed in Section I as order expansions in a perturbational strength parameter, z, we partition the electronic Hamiltonian according to Eq. (2.1). Next, we define the amplitude equations (Eq. (2.2)) needed for the evaluation of the energy corrections in the series ({t(z)} and ${𝐭∼(z)}$ amplitudes, respectively, for the HF- and CC-based series)

In Eqs. (2.3), the T(z), $T∼(z)$⁠, and $T*$ cluster operators are all generically defined as $T=∑i∑μitμiτ^μi$ (with $tμi$ and $τ^μi$ being the amplitude and excitation operator, respectively, for excitation $μi$ within the manifold at level i), and the ${*𝐭}$ amplitudes are those of parent model CC[$mP$] (for an HF parent state, these trivially vanish). Furthermore, the P- and Q-spaces define the so-called primary and secondary (complementary) excitation manifolds, such that the ${*𝐭}$ amplitudes have components only within the primary space and the {t(z)} amplitudes exist in the complementary space only (i.e., the primary space is empty for HF-based series), while the ${𝐭∼(z)}$ amplitudes have components in both of the two subspaces. For example, for the MP-like HF(SDT–n) series, the P-space will be empty (as is the case for all HF([$mQ$]–n) series) and the Q-space will be the complete manifold of all single, double, and triple excitations, while for the CCSD(T–n) series, the P-space will contain all single and double excitations whereas the Q-space will be restricted to the triples manifold only. A schematic illustration of some of these models is given in Figure 1.

Since we will search for intruder states by probing for avoided crossings between the ground state and an excited state in the interval $z∈[−1;1]$⁠, the degeneracies will show up as zero-valued excitation energies. As mentioned in Section II A, excitation energies are determined within CC response theory as eigenvalues of the CC Jacobian, which is defined as the derivative of the amplitude equations in Eq. (2.3) with respect to either of the {t(z)} or ${𝐭∼(z)}$ sets of amplitudes. The elements of the Jacobians for the HF([$mQ$]–n) and CC[$mP$]([$mQ$]–n) series are thus given as

where i,j may refer to any excitation level up to level [$mQ$]. As is evident upon comparing the two Jacobians for the HF([$mQ$]–n) and CC[$mP$]([$mQ$]–n) series in Eq. (2.4), these have the same structure, differing only in the amplitudes, {t(z)} and ${𝐭∼(z)}$⁠. Using the Baker-Campbell-Hausdorff expansion,^3,40 the difference between the matrix elements $Jμi,νj$ and $J∼μi,νj$ becomes

which, in general, makes the eigenvalues of the two Jacobians, i.e., the excitation energies, differ if the {t(z)} or ${𝐭∼(z)}$ amplitudes themselves differ notably. For z = 0 ∧ z = 1, however, we note how the excitation energies obtained from a diagonalization of Eqs. (2.4a) and (2.4b) will be identical; for z = 0, both of the J and $𝐉∼$ matrices are diagonal with orbital energy differences along the diagonal, and for z = 1, both return the excitation energies of target model CC[$mQ$].

To a first approximation and for a general z, the Jacobians in Eq. (2.4) are both dominated by the common contribution

which is independent of either set of amplitudes. In particular, we note from Eq. (2.6) how for $z<0$⁠, the unphysical action of the fluctuation potential will be to attract, rather than repel the electrons.^18,23 To see what ramifications this effect might have on the resulting excitation energies, we focus on the dominating diagonal elements of $𝐉com$ in Eq. (2.6), which may be written as

where $ϵμi$ is the (positive) difference in energy between the virtual and occupied orbitals of excitation $μi$⁠. For the diagonal elements in Eq. (2.7), the second term is the difference in energy between the determinant for excited state $μi$ and that for the HF ground state. Thus, $(⟨μi|H^(1)|μi⟩−⟨HF|H^(1)|HF⟩)$ will always be positive, and it will be large for a system with a very diffuse excited state (high excited state energy) and/or for an electron-rich system with a dense HF ground state (low HF ground state energy). Hence, assuming that the eigenvalues of $𝐉com$ will be dominated by the diagonal elements in Eq. (2.7), we may expect these to approach zero for general combinations of electron-rich ground states and diffuse excited states at $z<0$⁠, if the difference in electron repulsion energy counterbalances the orbital energy difference, i.e., if the following condition is met:

where $ϵμi$ increases as well when the virtual states involved are located high up in the energy spectrum. We will test this hypothesis in Section III by calculating excitation energies, not only from the two Jacobians in Eq. (2.4) but also from the dominant contribution to these, i.e., $𝐉com$ in Eq. (2.4). In passing, however, we note that the above reasoning is in line with that of Ref. 16, in which the presence of intruder states was also linked to differences in first-order MP corrections for the ground state and diffuse excited states, i.e., $(⟨μi|Φ^|μi⟩−⟨HF|Φ^|HF⟩)$⁠.

In calculating energy corrections to arbitrary order, the simplest procedure is to evaluate these using the standard n + 1 rule of Rayleigh-Schrödinger perturbation theory, which allows for the energy correction at order n + 1 to be calculated from the nth-order amplitude corrections. For the MP series, the resulting recursive algorithm, as formulated within an FCI program, has previously been described in, e.g., Refs. 7 and 9, and will not be repeated here. For the CC[$mP$]([$mQ$]–n) family of perturbation series, on the other hand, which was originally derived using Wigner’s 2n + 1 and 2n + 2 rules, arbitrary-order corrections are here computed using the same n + 1 rule as used for the MP series. Furthermore, we will also report numbers for the so-called E-CC[$mP$]([$mQ$]–n) family of series,²⁶ which too form order expansions of the difference in energy between a parent and a target CC model in orders of the fluctuation potential. For a given choice of CC[$mP$] and CC[$mQ$], the E-CC[$mP$]([$mQ$]–n) and CC[$mP$]([$mQ$]–n) series can both be derived from a bivariational energy Lagrangian for the target model CC[$mQ$], by using information exclusively of the right-hand CC state (E-CC[$mP$]([$mQ$]–n)) or a combination of information on the right- and left-hand (Λ) states (CC[$mP$]([$mQ$]–n)) of the parent CC[$mP$] model.

For the physical system (z = 1), the nth-order amplitude corrections, from which energy corrections to order n + 1 may be evaluated in either of the two series, are given as²⁶ 

where $Φ^T^*=exp(−*T)Φ^exp(*T)$ is the CC[$mP$] similarity-transformed fluctuation potential. For the sake of notational brevity, we have partitioned the correction amplitudes in Eq. (2.9) into two components, i.e., the corrections to the parent state amplitudes, ${*𝐭}$⁠, in the P-space (⁠$[P]$⁠) and the full amplitudes in the Q-space (⁠$[Q]$⁠),

In the CC[$mP$]([$mQ$]–n) and E-CC[$mP$]([$mQ$]–n) series, the difference in energy between the CC[$mP$] and CC[$mQ$] states is next expanded in terms of the E and $E¯$ corrections, respectively, according to the following relations:

where the nth-order energy corrections are defined as

In Eq. (2.12a), $⟨*Λ|=⟨HF|+⟨*λ|$⁠, where the states $⟨*λ|=∑p∈[P]∑νpλνp*⟨HF|τ^νp†$are spanned in terms of CC[$mP$] multipliers, and in Eq. (2.12b) and summation range for the second contribution is limited by the fact that amplitude corrections within the P-space vanish at first order, cf. Eq. (2.9a).

Finally, we note that the expansion of the cluster amplitudes in orders of the MP fluctuation potential is the same in the E-CC[$mP$]([$mQ$]–n) and CC[$mP$]([$mQ$]–n) series, see Eq. (2.9), and the two series therefore share the same Jacobian in Eq. (2.4b). Hence, the E-CC[$mP$]([$mQ$]–n) and CC[$mP$]([$mQ$]–n) series have identical convergence radii, albeit significantly different rates of convergence, as exemplified by the comparison of the E-CCSD(T–n) and CCSD(T–n) triples series in Ref. 26. This point will also be highlighted by the numerical results to follow.

In the present section, we will report the position of the nearest avoided crossing along the real axis within (or outside) the unit circle, as calculated using either J in Eq. (2.4a) for the HF([$mQ$]–n) series or $𝐉∼$ in Eq. (2.4b) for the CCSD([$mQ$]–n) series ([$mP$] = 2). As mentioned in the closing paragraph of Section I, the investigations will be made for four prototypical closed-shell examples, all of which are known to have slowly convergent or even divergent MP series upon moving to higher orders in the perturbation. The four examples are as follows: the Ne atom in the aug-cc-pVDZ basis⁴¹ (Section III A), HF (at a distorted geometry of twice the equilibrium bond length, $re=91.6$ pm) and CH₂ (⁠$A11$⁠) (⁠$re=110.7$ pm, $∠(HCH)=102.0°$⁠) in the cc-pVDZ basis⁶ (Sections III B and III C, respectively), followed by the F^– anion in the aug-cc-pVDZ basis (Section III D). The frozen-core FCI level is [$mQ$] = 6 for CH₂ and [$mQ$] = 8 for Ne, F^–, and HF. The main results are summarized in Table I, which reports whether a given series is convergent (con.) or divergent (div.) for a given truncation level, [$mQ$]. Furthermore, we will report total deviations from frozen-core FCI results for all four systems. The MP and CC[$mP$]([$mQ$]–n) series start at second order, while for the E-CC[$mP$]([$mQ$]–n) series, the leading-order correction is of third order, cf. Eq. (2.12). All of the involved calculations have been performed using the lucia program⁴² with the verification of the low-order results for the simplest expansions done using the aquarius program.⁴³ 

In Figure 2, we consider the convergence of the neon atom for various perturbation series. In spite of having a substantial weight of the HF determinant in the FCI state (with a resulting HOMO-LUMO energy gap in excess of 30 eV), the MP series for Ne was previously in Ref. 7 found to be divergent in the aug-cc-pVDZ basis, due to a weakly coupled back-door intruder with a real component, $Re(ξ)$⁠, of about –0.8, cf. Table I. If this intruder is decomposed into contributions from individual excitation levels, it is observed (in the MP series, i.e., with [$mQ$] = 8) to have more than 70% of the weight of its wave function assigned to hextuple and higher-level (seven- and eightfold) excited determinants. Thus, the intruder is clearly an unphysical and diffuse state, which will only appear in augmented basis sets that allow for such continuum states. However, by truncating either of the HF- and CC-based families of series at the level of hextuple excitations, the position of the avoided crossing is observed to move towards a lower z-value (⁠$Re(ξ)≃−0.9$ for [$mQ$] = 6), and by further truncating them at the level of pentuple (or lower-level) excitations, the crossing eventually falls outside the unit circle (i.e., $Re(ξ)<−1.0$ for [$mQ$] $<6$⁠). Thus, while the MP series does not converge for Ne in the aug-cc-pVDZ basis, for [$mQ$] $<6$⁠, any of its truncated HF([$mQ$]–n) series will, as will the corresponding CC[$mP$]([$mQ$]–n) (and E-CC[$mP$]([$mQ$]–n), cf. Section III A 2) series. The equivalence between the HF- and CC-based perturbation series in this case is further confirmed by calculations of excitation energies from the $𝐉com$ Jacobian in Eq. (2.6), for which the crossing is again observed to enter the unit circle upon an inclusion of hextuple excitations in the total excitation manifold, i.e., whenever [$mQ]≥6$⁠.

In the present section, we want to directly assess the effect of higher-level excitations and basis set diffuseness on the convergence behaviour of the CC[$mP$]([$mQ$]–n) series. First, in order to quantify that the description of hextuple excitations is indeed at the root of the divergence problems for Ne in the aug-cc-pVDZ basis, we probe for avoided crossings in the CCSDTQ5(6–n) series, that is, a series with a CC[5] parent state and a CC[6] target state, for which the difference in energy is a mere 0.0005 kcal/mol. These results are presented in Figure 3. In summary, an identical crossing with the same intruder state as in the case of the CCSD-based series in Figure 2 is observed within the unit circle (at $Re(ξ)≃−0.9$⁠), making the series divergent in spite of the minuscule magnitude of the energy difference which the series aims at expanding in orders of the perturbation.

Next, we assess whether other intruders but the hextuples-dominated state in the aug-cc-pVDZ basis might enter the unit circle upon an augmentation of the basis by an additional set of diffuse functions. In Figure 4, the convergence of the CCSD(TQ–n) series is shown in the aug-cc-pVDZ and d-aug-cc-pVDZ basis sets. As is clear, the convergent behaviour of the series in the former of the two basis sets is deteriorated in the latter, as a new intruder dominated by quadruple excitations now enters the unit circle. Thus, in general, a convergent HF- or CC-based series will likely become divergent if the excitation level of the target state is increased and/or the diffuseness of the basis set is increased.

As discussed in Section II C, theCC[$mP$]([$mQ$]–n) and E-CC[$mP$]([$mQ$]–n) series will exhibit identical convergence radii, due to their common Jacobian in Eq. (2.4b).However, as initially shown in Ref. 26 through results for thelowest-order corrections of the CCSD(T–n) and E-CCSD(T–n) series, the rate of convergence towards the CC[$mQ$] target energy will be different. In Figure 5, this is exemplified for Ne not only for the CCSD(T–n) and E-CCSD(T–n) series through higher orders but also for differing choices of CC[$mQ$] target states. From these results, we note how the oscillations through the lowest orders are considerably dampened in the CCSD([$mQ$]–n) series over the corresponding E-CCSD([$mQ$]–n) series, a direct consequence of the fact that the order expansions of the former family of series are markedly more balanced than those of the latter, cf. the discussions in Section II C and, in particular, Ref. 26. Thus, in terms of practical convergence, the CC[$mP$]([$mQ$]–n) variants will always be preferable.

In Ref. 16, the convergence of the MP series for hydrogen fluoride was investigated for combinations of two different geometries and two different basis sets. At the equilibrium geometry, the MP series was found to be rapidly convergent in the cc-pVDZ basis set and excessively oscillating (ultimately divergent) in the aug-cc-pVDZ basis. Upon distorting the geometry, however, the MP series was found to diverge in the standard cc-pVDZ basis as well, with an avoided crossing between the ground state and a back-door intruder well within the unit circle. As opposed to the intruder state in the case of the neon atom (Section III A), the back-door intruder for distorted HF is not dominated by an exorbitant amount of high-level excitations, but rather by a mixture of double, triple, and, in particular, quadruple excitations. In fact, the quadruples contribution to the description of the intruder state is so pronounced that if these are left out of the cluster operators in the CC-based series, the position of the avoided crossing is observed to shift from just within the unit circle, as for the CCSD(TQ–n) series in Table I and the CCSDT(Q–n) series in Figure 6, to a value slightly below z = –1 for the CCSD(T–n) series, also shown in Figure 6. Thus, whereas both HF- and CC-based expansions are in general divergent for distorted hydrogen fluoride, it is possible to form truncated expansions that do converge by choosing a target state that does not include quadruple (and higher-level) excitations. This conclusion is thus equivalent to the one for the neon atom in Section III A 1, for which the omission of hextuple and higher-level excitations in the cluster operator leads to a convergent series. For both of the CC-based series in Figure 6, however, a practical convergence towards the respective target energy (e.g., CCSDT or CCSDTQ) is observed, cf. Figure 7. Specifically, the divergence of the CCSDT(Q–n) series through the first 70 corrections is manifested only in very small energy oscillations of the order of 0.001 kcal/mol, which are negligible for all practical purposes and not even visible on the energy scale of Figure 7.

Turning next to the case of singlet methylene, we note how the radius of convergence of the MP series was previously found in Ref. 16 to be larger than 1.0 in the cc-pVDZ basis (as well as in the aug-cc-pVDZ basis), thereby validating a convergent MP series. However, an avoided crossing between the ground state and a low-lying doubly excited state of the same spin and spatial symmetry was found at a positive z-value of about 1.2 (that is, in the vicinity of, but still outside the boundary of the unit circle), which is confirmed for all of the tested perturbation series in Table I. In particular, the position of this crossing is not observed to shift significantly with changes in the truncation level, [$mQ$]. However, despite the fact that this excited state does not represent a front-door intruder state, the presence of the avoided crossing near the boundary of the unit circle is seen in Figure 8 to cause a monotonic but distinctly slow rate of convergence for all of the tested series, regardless of the choice of target state. In particular, we note how the shapes of the energy curves for singlet methylene in Figure 8 are completely different from the ones for Ne and HF, cf. Figures 2 and 6, differences which are intimately related to those of the energy profiles of front- and back-door intruders as described in, for instance, Ref. 13.

Finally, we close the present section with the most challenging of the four examples studied here, namely, the fluoride anion in the aug-cc-pVDZ basis set, for which the MP series was previously found, in Refs. 7 and 16, to diverge from the very onset of the expansion (from the MP3 model and onwards), essentially exploding up through higher orders. Although F^– is isoelectronic with Ne, and the weight of the HF determinant is approximately the same for the two species in the cc-pVDZ and aug-cc-pVDZ basis sets (93%–97%), being an anion, the need for augmented basis sets in the accurate description of F^– is of utmost importance. This is also recognized from the fact that the HOMO-LUMO gap in F^– is considerably lowered upon augmenting the cc-pVDZ basis by diffuse functions, from 55 eV in the cc-pVDZ basis to less than 20 eV in the aug-cc-pVDZ basis. In Ref. 16, an avoided crossing between the ground state and a relatively strongly coupled back-door intruder was observed at $Re(ξ)≃−0.6$ for the MP series, and the position of the crossing closer to z = 0—as well as the larger coupling between the two states—was used to explain the more rapid divergence for F^– over that for Ne.

In analogy with the studies of the neon atom in Section III A and distorted HF in Section III B, the wave function for the intruder state may be decomposed into contributions from individual excitation levels. By doing so, we find that more than 60% of the wave function consists of singly, doubly, and triply excited determinants. Furthermore, for an HF(SD–n) order expansion of the CCSD correlation energy, no back-door intruder is observed. Collectively, this indicates how the inclusion of triple excitations is the primary cause for the divergence of the HF(SDT–n) series, see Figure 9. However, in contrast to this result, we note from Figure 9 how the CCSD(T–n) triples series converges. This difference between the HF(SDT–n) and CCSD(T–n) series hence indicates that the choice of parent state may also impact the convergence behavior of CC perturbation series, despite being a minor overall factor compared to the influence of basis set diffuseness and level of target state. Relating this back to the analysis in Section II B, this implies that if pronounced differences exist between the amplitudes of T and $T∼$ (i.e., the difference in Eq. (2.5) is substantial), this may too impact the radius of convergence. The actual convergence of the CCSD(T–n) series for F^–, though, is observed to be highly irregular and almost pulse-like, as may be realized from the detailed view of the order corrections in Figure 10.

For the sake of completeness, it should also be noted that the lowest Jacobian eigenvalue obtained during the CCSD(T–n) z-scan for F^– contained a very small imaginary component for z-values in the interval [–1.0; –0.9]. However, based on the CCSD(T–n) energy corrections in Figures 9 and 10 alongside the fact that this imaginary component was found to be a full six orders of magnitude smaller than the (positive) real component, we conclude that the CCSD(T–n) series for F^– is indeed convergent. It should also be mentioned that the z-scan of the Jacobian for the CCSD(TQ–n) series was the only case in which the eigenvalue with the smallest real part had a large imaginary component. In this case, it is therefore not possible to use the z-scan on the real axis to determine whether the series is convergent or divergent. However, the perturbation expansion (not shown) was observed to diverge rapidly.

Of all the examples considered in the present work, we only encountered a single example—the HF(SDT–n) and CCSD(T–n) series for F^–—for which the choice of parent state could be attributed as having an effect on the convergence behaviour, while the choice of target state is generally a much more important factor (see Sections III A and III B). This contrast in significance is also illustrated by the fact that the CCSD(TQ–n) and CCSDT(Q–n) series for F^– both diverge, see Table I and Figure 9. It is, however, worth noting—by pragmatically comparing the corrections of the various series in Figure 9—how the use of the lowest-order corrections of any of the CC-based series is not invalidated in the case of F^–, as opposed to the situation for HF-based analogues, for which only the second-order (MP2) correction is typically considered meaningful.

Now, the special case of the convergent CCSD(T–n) series aside, the results in Figure 9 highlight how the majority of the problems for F^– originate from the underlying MP partitioning of the electronic Hamiltonian used in both the HF- and CC-based series, that is, the division of $H^$ into the zeroth-order Fock operator and a sizeable perturbation in the form of the fluctuation potential. This is so because the combination of an electron-rich ground state and an exceedingly diffuse excited state (archetypal for anions in augmented basis sets) makes for dominant contributions to either of the J and $𝐉∼$ Jacobians that induce avoided crossings between such states for negative values of z. This notion is further substantiated by probing for negative excitation energies using the $𝐉com$ Jacobian in Eq. (2.6), for which the same crossing as in the HF(SDT–n) and CCSDT(Q–n) series is again observed.

The radius of convergence of various HF- and CC-based MP perturbation expansions with different CC target states has been determined for a selection of prototypical closed-shell examples that are known from the literature to have slowly convergent or even divergent MP series. In particular, we have probed for the presence of potential intruder states in truncated MP series as well as in the recently proposed CCSD(T–n) triples and CCSD(TQ–n)/CCSDT(Q–n) quadruples series, in turn by probing for zero- or negative-valued excitation energies as calculated from the CC Jacobian for each of the series. The similarities between the intruder states encountered in the HF- and CC-based perturbation theory have been interpreted in terms of structural similarities between the resulting Jacobians for the different series. By focusing on the common and dominating contribution to these, we have been able to confirm and supplement the main conclusions of the original MP convergence study in Ref. 16.

We have detailed how perturbation theory formulated around an MP partitioning of the electronic Hamiltonian, be that based on an HF or a CC ground state wave function, will be prone to back-door intruder states whenever a description of diffuse continuum states is possible, as, for instance, is the case when an augmented basis set is used or the target state contains highly excited determinants. Thus, in general, a divergent HF- or CC-based series for a given combination of molecular structure and basis set can likely be made convergent by reducing the excitation level of the target state and/or reducing the diffuseness of the basis set. On the other hand, the choice of parent state has typically (but not always) no effect on the convergence behaviour.

Furthermore, we have reiterated the statement saying that back-door intruders are indeed prevalent for electron-rich systems with dense ground state wave functions. From a comparison of the Jacobians of HF- and CC-based perturbation theory, both of these conclusions have been traced back to the fundamental MP partitioning of the Hamiltonian, namely, the fact that the perturbation within this framework, i.e., the fluctuation potential, is substantial in size, in contrast with the general premise of perturbation theory where the perturbation in itself is assumed small. Thus, we argue that the observed divergences are artefacts of the MP partitioning and, as such, inherent to any perturbation expansion that makes use of this. Perhaps somewhat counterintuitive, we have also illustrated how the magnitude of the actual energy difference, which is expanded in orders of the perturbation in any given perturbation series, is completely irrelevant for the actual convergence behaviour of said series.

However, divergences aside, the CC-based perturbation series have all been observed to behave significantly more stably than the MP series or any of its truncated variants. Furthermore, we have found the bivariational CCSD(T–n) and CCSDT(Q–n) series to be capable of remedying the disordered oscillations at lower orders in the perturbation—which is traditionally one of the premier characteristics of divergent perturbation expansions—even in cases where these series formally diverge at higher orders. As such, despite the fact that the present analysis has been limited to minimal systems (atoms and diatomics), we find that even when potential divergences exist, as, for instance, in the case of the CCSDT(Q–n) series for stretched hydrogen fluoride and the fluoride anion, these are not found to invalidate the use of lower-order models such as the CCSDT(Q–3) and CCSDT(Q–4) models. This is indeed an important point to stress, as formal divergences remain the norm rather than the exception for perturbation expansions in orders of the MP fluctuation potential, irrespective of whether these are based on an uncorrelated HF state or a correlated CC state.

J.J.E., K.K., and P.J. acknowledge the support from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 291371, D.A.M. acknowledges the support from the U.S. National Science Foundation (NSF) under Grant No. ACI-1148125/1340293 and from the Arnold and Mabel Beckman Foundation as an Arnold O. Beckman postdoctoral fellow, and J.O. acknowledges the support from the Danish Council for Independent Research, Grant No. DFF-4181-00537.