Variational coupled cluster for ground and excited states

In single-reference coupled-cluster (CC) methods, one has to solve a set of non-linear polynomial equations in order to determine the so-called amplitudes which are then used to compute the energy and other properties. Although it is of common practice to converge to the (lowest-energy) ground-state solution, it is also possible, thanks to tailored algorithms, to access higher-energy roots of these equations which may or may not correspond to genuine excited states. Here, we explore the structure of the energy landscape of variational CC (VCC) and we compare it with its (projected) traditional version (TCC) in the case where the excitation operator is restricted to paired double excitations (pCCD). By investigating two model systems (the symmetric stretching of the linear \ce{H4} molecule and the continuous deformation of the square \ce{H4} molecule into a rectangular arrangement) in the presence of weak and strong correlations, the performance of VpCCD and TpCCD are gauged against their configuration interaction (CI) equivalent, known as doubly-occupied CI (DOCI), for reference Slater determinants made of ground- or excited-state Hartree-Fock orbitals or state-specific orbitals optimized directly at the VpCCD level. The influence of spatial symmetry breaking is also investigated.


I. COUPLED CLUSTER AND STRONG CORRELATION
Single-reference (SR) coupled-cluster (CC) methods offers a reliable description of weakly correlated systems through a well-defined hierarchy of systematically improvable models. [1][2][3][4][5] On top of this hierarchy stands full CC (FCC), which is equivalent to full configuration interaction (FCI), and consequently provides, at a very expensive computational cost, the exact wave function and energy of the system in a given basis set. Fortunately, more affordable methods have been designed and the popular CCSD(T) method, which includes singles, doubles and non-iterative triples, is nowadays considered as the gold standard of quantum chemistry for ground-state energies and properties. 6,7 Despite its success for weakly correlated systems, it is now widely known that CCSD(T) flagrantly breaks down in the presence of strong correlation as one cannot efficiently describe such systems with a single (reference) Slater determinant. This has motivated quantum chemists to design multi-reference CC (MRCC) methods. [8][9][10][11][12] However, it is fair to say that these methods are computationally demanding and still far from being black-box.
Because SRCC works so well for weak correlation, it would be convenient to be able to treat strong correlation within the very same framework. This is further motivated by the fact that one can compensate the poor quality of the reference wave function by simply increasing the maximum excitation degree of the CC expansion. However, this is inevitably associated with a rapid growth of the computational cost, and hence one cannot always afford this brute-force strategy. The development of SR-based methods for strong correlation is ongoing and some of them (usually based on the "addition-by-subtraction" principle) have shown promising results. A non-exhaustive list includes pair coupledcluster doubles, [13][14][15][16][17][18][19][20][21][22] singlet-paired CCD, 23,24 the distinguishable cluster methods, [25][26][27][28][29][30][31][32][33][34] CCD-based variants involving a welldefined subset of diagrams, [35][36][37][38][39] the nCC hierarchy, 40,41 and parametrized CCSD. 42 Each of these methods sheds new light on the failures of SRCC to treat static correlation. For the sake a) Electronic mail: loos@irsamc.ups-tlse.fr of brevity, we omit the single-reference prefix hereafter.
The CC wave function |Ψ CC is obtained by applying a wave operator onto a single Slater determinant reference |Ψ 0 as In CC theory, the wave operator is defined as the exponential of the cluster operatorT which is the sum of the kth-degree excitation operator up to k = N (where N is the number of electrons). In second quantized form, we haveT a i and a † a being the second quantization annihilation and creation operators, respectively, which annihilates (creates) an electron in the spin-orbital i (a). The cluster amplitudes t ab... i j... are the quantities of interest in order to compute the CC energy (see below).
Throughout the paper, p, q, r, and s denote general spinorbitals, i, j, k, and l refer to occupied spin-orbitals (hole states) and a, b, c, and d to unoccupied spin-orbitals (particle states).
In quantum mechanics, one convenient way to determine the parameters of a wave function ansätz is to minimize the energy with respect to its parameters. The Rayleigh-Ritz variational principle ensures that the energy thus obtained is an upper bound to the exact ground-state energy. Following this strategy, the variational CC (VCC) energy [43][44][45][46][47][48][49][50][51] is thus minimized with respect to the cluster amplitudes which ensures min t ab... i j...
Unfortunately, independently of the excitation rank ofT , this procedure is not tractable in practice. Indeed, because the series expansion of the exponential in Eq. (4) does not truncate before the Nth-order term, VCC has an inherent exponential scaling with respect to system size. Usually, one sacrifices the attractive upper bound property of the variational principle in exchange for computational tractability. To do so, the similarity-transformed Schrödinger equation is projected onto the reference determinant |Ψ 0 , which gives This energy can be seen as the expectation value of a similaritytransformed HamiltonianH = e −TĤ eT for the reference determinant |Ψ 0 . One can expandH thanks to the Baker-Campbell-Hausdorff formula and show that this series naturally truncates after the fourth-order term. This truncation is due to the twoelectron nature of the Hamiltonian and is responsible for the affordable polynomial scaling of this method (contrary to the exponential cost of the variational approach). In such a case, the cluster amplitudes are no longer determined by minimization of the VCC energy functional (4) but via the amplitude equations which are the projection of the similarity-transformed Schrödinger equation (6) onto excited determinants. In Eq. (8), the determinant |Ψ ab... i j... is obtained by promoting the electrons occupying the orbitals i, j, . . . in |Ψ 0 to the vacant orbitals a, b, . . . . One usually refers to this type of methods as traditional CC (TCC).
As reported in Refs. 48,52, and 53, VCC has been shown to give correct results in situations where TCC fails. These benchmark studies evidenced that the breakdowns of TCC cannot be explained solely by its single-reference nature, as part of the problem actually originates from its non-variational character. Unfortunately, because of the exponential scaling of VCC, it is computationally cumbersome and cannot be applied in practice except for small molecules in small basis sets. This drawback has motivated the search for approximate methods that retain the advantages of VCC but at a polynomial cost. Because VCC inherits its exponential scaling from the lack of truncation of its energy functional [see Eq. (4)], some authors have designed ingenious truncation schemes. 43,44 The quasivariational CC (QVCC) method from Knowles' group has been designed along these lines. 50,[54][55][56][57] This method, which is an improvement of the former linked pair functional, 49 can be seen as an infinity summation of a given subset of diagrams of the VCC energy functional. This method has most of the desirable properties of an approximate VCC theory [see Ref. 54 for an exhaustive discussion of these properties] but is not an upper bound to the exact energy. Yet, QVCC has been proved to be much more robust than TCC in cases where the latter exhibits non-variational collapse below the FCI energy like, for example, in the symmetric bond stretching of the nitrogen and acetylene molecules. [55][56][57] Using VCC instead of TCC has its advantages but its computational complexity is very nettlesome. It would be simpler if one could describe strong correlation while retaining the projective way of solving the equations and its associated polynomial cost. Surprisingly, restricting the cluster operator to paired double excitations (pCCD), which is a simplification with respect to CC with doubles (CCD), 58 can give qualitatively good results for strongly correlated systems. [13][14][15][16][17][18][19][20][21][22]24 This can be understood thanks to the concept of seniority number which is defined as the number of unpaired electrons in a determinant. 59 Indeed, the seniority-zero subspace (i.e., the set of all closedshell determinants) has proven to give a good description of static correlation. 60 Unfortunately, doubly-occupied configuration interaction (DOCI), which is a CI calculation in the seniority-zero subspace, inherits the exponential scaling of its FCI parent. [60][61][62][63][64][65][66] However, benchmark results 15,16,18,67 have shown that pCCD provides ground-state energies which are almost indistinguishable from the DOCI ones but at a mean-field computational cost, hence providing a tractable way to qualitatively describe strongly correlated systems. Note that pCCD is equivalent to the antisymmetric product of 1-reference orbital geminals (AP1roG) 13,14,21,[68][69][70][71][72][73][74][75][76][77] which has been designed as a computationally tractable approximation to the antisymmetric product of geminals (APG), 78,79 a method that has been recently further explored by the group of Scuseria. [80][81][82][83][84][85] Because the seniority-zero subspace is not invariant to orbital rotations, one must energetically optimize the orbitals to obtain the optimal pairing scheme, (i.e., the orbital set that minimizes the energy in the seniority-zero subspace). 60 In Ref. 13, Limacher et al. determined this pairing scheme by optimizing the orbitals at the DOCI level and then using these orbitals for their geminal wave function methods. Later, Henderson et al. designed an orbital-optimized pCCD (oo-pCCD) procedure which provides a more straightforward route to obtain this optimal pairing scheme. 15

A. TCC for excited states
Excited-state energies and properties can be computed within the TCC paradigm through the well-established equation-of-motion (EOM) formalism. [86][87][88][89][90] In EOM-CC, one applies a suitably chosen (linear) excitation operator on a ground-state CC wave function to compute excited states. This procedure can be conveniently recast as a non-Hermitian eigenvalue problem involving the similarity-transformed Hamilto-nianH in a space of excited determinants. 5 Like in ground-state TCC, one can systematically expand the excitation space to form a well-defined hierarchy of EOM methods. As an example, EOM-CCSD restricts the set of excited determinants to singles and doubles. EOM-CCSD is known to accurately describe single excitations 91,92 but dramatically fails to describe double excitations because of the lack of triples and higher excitations. 93,94 This shortcoming can be corrected by the in-clusion of these higher excitations but this is not without a steep increase of the computational cost. [95][96][97][98][99][100] Albeit being by far the most popular excited-state formalism, EOM is not the only route to excited states within CC theory. Indeed, the amplitude equations (8) constitute a set of nonlinear polynomial equations and consequently possess many solutions. These solutions, sometimes labeled as "non-standard", can be non-physical or correspond to genuine excited states. 101 Therefore, performing a first standard ground-state CC calculation and a second one converging towards a given excited state provides an alternative way to obtain excitation energies. 102,103 Lee et al. refers to this type of methods as ∆CC 103 by analogy with the ∆SCF methods where one basically follows the same procedure but at the self-consistent field (SCF) level. Indeed, the use of Hartree-Fock (HF) or Kohn-Sham higher-energy solutions corresponding to excited states is becoming more and more popular and new algorithms designed to target such solutions, like the maximum overlap method (MOM) [104][105][106][107] or more involved variants, [108][109][110][111][112][113][114][115][116][117][118][119][120][121][122] are being actively developed. Besides providing a qualitatively good description of excited states, 122 these solutions can also be very helpful for ∆CC methods, as we shall illustrate below (see also Ref. 103).
The set of orbitals used, particularly the orbitals that constitute |Ψ 0 , strongly influences the performance of ∆CC methods. Importantly, the use of state-specific orbitals plays the role of a magnifying glass and facilitates the convergence towards a given CC solution by enlarging the associated basin of attraction. In addition to the orbital set, two other factors influence in a significant way the solutions that can be reached: the guess amplitudes for the CC equations and the algorithm employed for solving these equations. Even if the chosen orbitals can enlarge or shrink the basin of attraction of a given solution, one still has to pick an appropriate starting point within this basin to be able to converge to the desired solution. Moreover, the type of iterative algorithms (usually based on Newton-Raphson method and/or supplemented by Pulay's DIIS method [123][124][125] must also be carefully chosen so as to target, for example, saddle points or maxima instead of minima. For example, as shown in Ref. 126, the usual CC iterative algorithm is inappropriate to converge towards excited states. Because of the non-linearity of the CC equations [see Eq. (8)], the number of solutions can be higher that the physically meaningful number. However, claiming that a given solution corresponds to a genuine electronic state (or not) is a rather tricky task as the overall picture behind the structure of the CC solutions is still far from being thoroughly understood. Zivkovic and Monkhorst were the first to tackle this outstanding problem with their seminal work on the existence conditions of the higher roots of the CC equations. 127,128 However, their model was too simplistic and most of the pathological solutions that they found or predicted were due to this unrealistic model as argued later by Jankowski et al. who investigated the CCD solutions of 1 A 1 symmetry in the H 4 molecule. [129][130][131] Still, they evidenced that some non-standard solutions may be non-physical. They also showed that the CC solution structure highly depends on the reference. 131 Few years later, the introduction of the homotopy method (which gives all the solutions of a set of non-linear equations) in the CC paradigm enabled the first systematic study on the structure of the CC energy landscape. 132,133 In particular, these studies showed that, in practice, the number of CC solutions is much lower than the theoretical upper bound known as Bézout's number. We refer the interested reader to the series of papers by Jankowski et al. [134][135][136][137] and the book chapter of Piecuch and Kowalski 101 for an extensive discussion about the homotopy method and the higher-energy solutions of the CC equations. We should also mention that the homotopy method has been employed to locate the CC solutions of the PPP model for some cyclic polyenes, 138,139 as well as in the context of MRCC and the Bloch equation formalism. [140][141][142] More recent studies have further improved our understanding of the CC energy landscape from which multiple solutions emerge. 103,143 As pointed out by Mayhall in their study of the CCSD solutions in the NiH molecule, the problem of the CC solution structure still needs to be addressed for more realistic systems. 143 Lee et al. showed that ∆CC can provide fairly accurate double excitation and double core-hole energies. 103 Recently, we have pursued along these lines by analyzing the non-standard solutions of the pCCD equations. We have shown that the agreement between pCCD and DOCI holds for excited states on the condition that state-specific optimized orbitals are employed. 126 Moreover, Ref. 126 brought some answers to Mayhall's open question as we have shown that ∆oo-pCCD provides double excitation energies that are comparable in terms of accuracy to the more expensive EOM-CCSDT method 97,98,[144][145][146] for a set of small (yet realistic) molecules. It is worth mentioning again that all the studies mentioned above deal with TCC methods. Note also that Ref. 21 discusses, in particular, alternative ways to find multiple excited states for AP1roG (and related methods).

B. VCC for excited states
For the sake of clarity, from here on, we restrict ourselves to VCCD (i.e.,T =T 2 ) but the procedure presented below is general and can be applied to higher-order variants. To the best of our knowledge, the present study is the first one to investigate excited states at the VCC level.
Because saddle points and maxima of the HF energy functional represent excited states, one can genuinely wonder if the same holds for the VCC energy functional (4). Thus, we seek for its stationary points, i.e., the different sets of cluster amplitudes t with elements t ab i j satisfying where the VCCD residuals r ab i j are the elements of the tensor r. The ground-state variational solution obtained via the minimization of Eq. (4) is also a solution of the more general equations (9) which provide all the stationary solutions of the VCCD equations. In this study, we restrict ourselves to solutions with real cluster amplitudes. The explicit expressions of the residual equations under this assumption are derived in Appendix A. Of course, stationary points of the VCC energy functional associated with complex cluster amplitudes may also exist. Indeed, the hermiticity of eT †Ĥ eT ensures that E VCC is real for any set of amplitudes. 44 Because VCC has an inherent exponential scaling, one can take advantage of the more convenient FCI representation to implement VCC algorithms. 48,52 Following Van Voorhis and Head-Gordon, we represent the (unnormalized) CC wave function as a CI vector (i.e., in the Slater determinant basis) by N successive applications of the cluster operatorT on the reference wave function: Using this CI representation, the action of second quantized operators on the CC wave function is quite straightforward, and one can evaluate the energy (4) and the residuals (A1) by simple matrix products. Note that the coefficients of the resulting CC wave function [see Eq. 10] are equal to the cluster analysis of the CI coefficients. 87,[147][148][149] In their VCCD benchmark study, Van Voorhis and Head-Gordon 52 relied on the standard TCCD iterative procedure (where one computes an approximate diagonal Jacobian matrix based on the difference of the Fock matrix elements f q p ) to solve Eq. (9): However, this approximate form of the Jacobian matrix cannot be employed to target excited states as it systematically converges towards the ground state or eventually diverges (see Ref. 126 for an exhaustive discussion on this point). If one aims at excited states, one should be aware that they generally are saddle points of the energy landscape. However, local minima could also correspond to physical excited states but it was not the case for the two model systems considered here. Therefore, to target saddle points, one should take into account the curvature of the energy landscape. To do so, we consider the Jacobian matrix J with elements which is then used to update the amplitudes according to the usual Newton-Raphson algorithm, i.e., The general expression of the Jacobian matrix elements is given in Appendix A. Note that this updating scheme of the amplitudes [see Eq. (13)] is more computationally demanding than the usual one [see Eq. (11)] as it requires to compute the entire Jacobian matrix and invert it. We should nonetheless mention that alternative algorithms are available to target such solutions. For example, approximate Newton-Raphson schemes which preserve the information about the energy curvature contained in the exact Jacobian matrix, or large-scale iterative solvers where one does not need to construct the full Jacobian, could also be employed. In difficult cases, it can be useful to damp the Newton-Raphson steps. However, one has to ensure that the structure of the Jacobian matrix is preserved during this process. This can be done by diagonalizing the Jacobian and adding a positive/negative constant to the positive/negative eigenvalues, similarly to what we have recently done for orbital optimization at the pCCD level. 126 To fully specify our algorithm, we still need to choose our reference |Ψ 0 as well as the starting values of the cluster amplitudes. In this study, we rely on both ground-and excited-state HF wave functions as references in order to study the influence of state-specific references. The orbitals employed to construct these excited-state HF wave functions have been obtained using initial MOM (IMOM). [104][105][106][107] State-specific orbitals optimized at the correlated level are also considered, as discussed below. Regarding the starting values of the cluster amplitudes t, once again we have taken advantage of the FCI representation by obtaining these via a cluster analysis of the corresponding CI eigenvectors. 87,148

C. Orbital optimization for excited states
The solutions obtained via this iterative process [see Eq. (13)] are stationary points of the VCCD energy functional with respect to the cluster amplitudes but not with respect to the orbital coefficients. Indeed, the orbitals have usually been obtained at the HF level and no longer represent a stationary point when electron correlation is introduced. The next step is thus to optimize the orbitals at the corresponding correlated level to find solutions that are stationary with respect to both the cluster amplitudes and the orbital coefficients.
As usually done, 150,151 we introduce a unitary operator eˆκ into the VCCD energy functional, to account for orbital rotations. Now, Eq. (14) can be minimized with respect to the cluster amplitudes t ab i j and to the orbital rotation parameters κ pq of the one-electron anti-Hermitian operatorκ. For a given set of cluster amplitudes, we search for the stationary points with respect to the orbital rotation parameters using the second-order Newton-Raphson method. We then expand the VCC energy around κ = 0, where g is the orbital gradient and H is the orbital Hessian, both evaluated at κ = 0, i.e., The orbitals are then updated following the usual Newton-Raphson step where C is the orbital coefficient matrix. Then, one finds the solution of Eq. (9) for this new set of orbitals and the procedure is repeated until convergence. To compute the gradient and the Hessian, one must compute the one-and two-body density matrices, 15 with respective elements where the orbital index refers to spatial orbitals, and σ and σ to spin indexes. Once again, we take advantage of the CI representation of the VCCD wave function to compute these quantities. We express the string of second quantized operators in Eqs. (18a) and (18b) as a matrix in the Slater determinant basis, and then evaluate the elements of the one-and two-body density matrices by simple matrix products.
In the present study, we restrict the cluster operator to a pair double excitation operator (with P † q = a † q↑ a † q↓ ) and investigate the properties of ground and excited states at the traditional pCCD (TpCCD) and variational pCCD (VpCCD) levels. This choice is motivated by the two following arguments. Firstly, our aim is to compare the VpCCD solution structure with its TpCCD counterpart (which has received our attention recently 126 ) in order to provide new insights into the multiple solutions of the VCC equations. Secondly, this restriction of the cluster operator significantly lowers both the computational cost and the complexity of the energy landscape, hence simplifying the present analysis. The VpCCD equations are easily obtained from their VCCD analogs (see Appendix B for their explicit expressions). We refer the interested reader to Ref. 15 for a complete list of equations and an exhaustive discussion of the orbital optimization algorithm in the case of ground-state TpCCD and to Ref. 126 for the case of excited-state TpCCD.
In the following, taking the symmetric dissociation of the linear H 4 molecule as a first case study, ground-and excitedstate energies obtained at the TpCCD and VpCCD levels are compared to DOCI for three different sets of orbitals: groundstate HF orbitals, state-specific HF orbitals and state-specific orbitals optimized at the VpCCD level. In a second stage, we look at the various TpCCD, VpCCD and DOCI electronic states in the presence of strong correlation (i.e., near degeneracies) by examining the continuous deformation of H 4 from a square to a rectangular arrangement.

III. COMPUTATIONAL DETAILS
The computational methods investigated here (HF, MOM, TpCCD, VpCCD, DOCI, and FCI) have been implemented as standalone mathematica modules, 152 which makes them easily interconnectable and modifiable depending on the actual purpose.
These are provided in an accompanying notebook available for download from Zenodo at http://doi.org/10.5281/zenodo.4971904. All the calculations have been performed in the restricted formalism. The only required input is the one-and two-electron integrals which are usually computed with a third-party software like quantum package. [153][154][155] The convergence threshold (based on the DIIS commutator) was set to 10 −10 a.u. for the restricted HF (RHF) calculations, while the convergence thresholds (based on the maximum absolute value of the gradient) for the cluster amplitude and orbital optimization procedures were both set to 10 −6 a.u.

A. Influence of the orbital set: the linear H 4 molecule
As a first example, we consider the symmetric stretching of the linear H 4 molecule in a minimal basis (STO-6G 156 ). This corresponds to a system with 4 electrons in 4 spatial orbitals with respective symmetries σ g , σ u , σ * g , and σ * u (in ascending energies). Linear chains of hydrogens are prototypical examples of left-right correlation and, therefore, have been widely studied in order to probe electronic structure methods in presence of such correlation. 13,15,25,57,60,[157][158][159][160][161][162][163][164] Hereafter, the distance between two consecutive hydrogens is denoted by R.
The first stage of this study consists in investigating the quality of the TpCCD and VpCCD ground-and excited-state energies in the case where the reference wave function is chosen as the ground-state RHF determinant, a choice that obviously induces a bias towards the ground state. The VpCCD energies (solid lines) are plotted alongside the DOCI ones (markers) in the left-hand-side of Fig. 1. Thanks to the simplicity of this model, one can access, via mathematica's implementation of the Jenkins-Traub algorithm, 165,166 the entire set of solutions (with real cluster amplitudes) associated with the system of polynomial equations (9). These VpCDD solutions are represented as thin solid lines in Fig. 1, while the thick parts of the curves correspond to the energies that we have been able to obtain using the Newton-Raphson algorithm described earlier [see Sec. II C]. Figure 1 also shows the TpCCD energies (dashed lines) which are also determined with the Jenkins-Traub algorithm applied to Eq. (8). In addition, the difference between TpCCD (VpCCD) and DOCI energies are also plotted in the top (bottom) right panel of Fig. 1.
Considering the ground-state RHF determinant as reference wave function, the convergence towards the VpCCD ground state, (σ g ) 2 (σ u ) 2 , is numerically straightforward all along the potential energy curve (PEC). On the other hand, converging excited states with the Newton-Raphson algorithm has been found to be much more challenging. We have not been able to converge the two lowest VpCCD excited states, (σ g ) 2 (σ * g ) 2 and (σ u ) 2 (σ * g ) 2 , further than R = 2.0 a 0 for this set of orbitals. Even worse, the two other doubly-excited states, (σ g ) 2 (σ * u ) 2 and (σ u ) 2 (σ * u ) 2 , have been reached only for R = 1.0 a 0 . This is not the case for the (σ * g ) 2 (σ * u ) 2 quadruply-excited state for which one can converge VpCCD calculations fairly easily all along the PEC with the Newton-Raphson algorithm. This might be due to the fact that the corresponding stationary points are maxima for the quadruply-excited state whereas doublyexcited states correspond to saddle points (see below). Despite such numerical difficulties, the complete set of solutions could be obtained thanks to the Jenkins-Traub algorithm.
Overall, the VpCCD method provides a fairly good approximation to the DOCI energies. At small R (i.e., in the weak correlation regime), VpCCD is in much closer agreement with DOCI than TpCCD, most noticeably for the (σ * g ) 2 (σ * u ) 2 quadruply-excited state. At large R (i.e., in the strong correlation regime), this comparison is trickier. Yet, the difference between VpCCD and DOCI seems more regular (see the bottom-right panel of Fig. 1) whereas the behavior of TpCCD is more erratic (top-right panel). Thus, one can state that, if one considers the ground-state RHF determinant as reference wave function, the main difficulties associated with VpCCD calculations concern its convergence, as the energies compare very favourably with the DOCI reference (at least in the weak cor-relation regime). At large R, the agreement between VpCCD and DOCI is less obvious as we shall see below.
Thanks to previous investigations, we know that some of the TpCCD excited-state solutions can be labeled as nonphysical. 101,126,[129][130][131][134][135][136][137] For example, in the case of the linear H 4 molecule using the ground-state RHF determinant as reference wave function, the lowest-lying DOCI excited state (blue markers in Fig. 1) can be represented by two TpCCD solutions (dashed blue curves). 126 These solutions eventually merge for R > 1.7 a 0 and become a complex conjugate pair of solutions with real components represented as black dashed lines in Fig. 1. The same phenomenon occurs for the fourth doubly-excited state, but the complex conjugate pair of solutions exists up to R = 3.4 a 0 before splitting into two real solutions (dashed green curves). In the case of VpCCD, there are only six real-amplitude solutions in the weak correlation regime. However, for R 3.5 a 0 , two additional real solutions appear as one can see in the inset of Fig. 1. The fact that these spurious solutions appear as a pair indicate that they may exist for smaller R as a pair of solutions with complex conjugate amplitudes. Because all the solutions are energetically close in this region of the PECs, it is hard to tell whether a solution is unphysical or not, and which one models better the corresponding DOCI solution. This is why the curve corresponding to the difference between VpCCD and DOCI for the fourth doubly- excited state stops at R = 3.2 a 0 in the bottom right panel of Fig. 1. The same unpredictability occurs between the green and purple TpCCD curves in the strong correlation regime. Therefore, in the weak correlation regime, the problems caused by unphysical solutions seem to be less severe in VpCCD than in TpCCD. Yet, when the correlation becomes strong, VpCCD is also plagued by unphysical solutions. Note that unphysical solutions at the VpCCD level are due to the approximate nature of the method which originates from the truncation of the cluster operatorT . On the other hand, unphysical TpCCD solutions can originate from this same truncation and/or from the projection step of Eq. (7). The stability analysis of the various stationary points via the computation of the eigenvalues of the Jacobian matrix (12) provides useful information on the presence of additional solutions. 167,168 For example, a change in the number of negative eigenvalues (the saddle point index) indicates the appearance of additional solutions. 169 For R < 3.4 a 0 , the index of the VpCCD solutions, in ascending energies, increases smoothly (0, 1, 2, 2, 3, and 4) up to the cyan curve which is an index-4 stationary point (i.e., a maximum). At R = 3.4 a 0 , the index associated to the green VpCCD solution decreases by one unit, this solution becoming an index-1 saddle point. We see in the inset of Fig. 1 that two additional solutions appear right after this index variation, these two spurious solutions being index-3 saddle points.
Because the agreement between DOCI and both TpCCD and VpCCD ground-state energies is very satisfying when one employs as reference the ground-state RHF wave function, one can reasonably wonder if the same similarity holds in the case of excited states by considering excited-state RHF wave functions as state-specific references. We have recently shown that this holds for TpCCD when one uses state-specific orbitals optimized at this correlated level, 126 but it remains to be seen whether this still applies with state-specific mean-field orbitals. Using IMOM, 170 we have obtained five additional restricted solutions of the RHF equations, corresponding to the five possible non-Aufbau closed-shell determinants (see Fig. 2). Note that each excited-state solution corresponds to a different set of orthogonal orbitals, but these sets are, a priori, not orthogonal to each others because they originate from distinct Fock operators. Of course, spatially symmetry-broken RHF solutions do exist but we have not considered them here. For R > 2.4 a 0 , an additional solution, plotted as a black line in Fig. 2, has been found by systematically occupying the two highest-energy molecular orbitals at each SCF cycle. The molecular orbitals associated with this H + H -H -H + ionic configuration have a more localized character than the orbitals constituting the (σ * g ) 2 (σ * u ) 2 determinant (the orbitals associated with these two solutions are available in the supplementary material). A stability analysis of these RHF solutions [171][172][173] shows that the cyan curve is a maxima at small R but, for R > 2.4 a 0 , it becomes a saddle point whereas the ionic configuration (i.e., the black curve) corresponds to a maxima.
The MOM excited states represented in Fig. 2 can be used as reference wave functions at both the TCC (see Ref. 103) and VCC levels. The energies at the three different correlated levels (namely, DOCI, TpCCD, and VpCCD) using these statespecific excited-state RHF reference wave functions are plotted in Fig. 3 and are labeled as MOM-DOCI, MOM-TpCCD, and MOM-VpCCD in the following. As one can see, these energies are visually indistinguishable, except at large R in the case of the (σ * g ) 2 (σ * u ) 2 state. Still, the right panel of Fig. 3 shows that MOM-VpCCD is closer to MOM-DOCI than MOM-TpCCD by roughly one order of magnitude all along the PEC. Also, we can see in the top-right panel that the difference between MOM-TpCCD and MOM-DOCI is less erratic than its analog using ground-state RHF orbitals (see Fig. 1).
As expected, using state-specific RHF determinants as reference rather than the ground-state one significantly improves the description of excited states at the TpCCD level. Therefore, if one wants to target excited states at the TpCCD level, it is worth investing in the design of proper state-specific references in order to make the key projection step in Eq. (7) more effective. Even if it is less pronounced at the VpCCD level, state-specific RHF reference determinants also improve the accuracy of the excited-state energies (with respect to DOCI). The most noticeable positive side effect of these state-specific references on VpCCD is the greater ease of convergence. Indeed, as shown in Fig. 3, we have been able to converge almost all the states up to R 3.5 a 0 . Therefore, we argue that using state-specific references enlarge the basin of attraction of the associated solution and consequently facilitates the convergence towards it.
The logical next step is to compare DOCI, TpCCD, and VpCCD at the orbital-optimized (oo) level (as described in Sec. II C). For ground states, DOCI and pCCD have already been shown to perform best when one relaxes the spatial symmetry constraint as it allows the orbitals to be fully localized. 14,70,71 On the other hand, relaxing this constraint also considerably increases, in principle, the number of attainable solutions. For example, multiple solutions corresponding to the ground state have already been observed. 14,70,71 However, in the case of excited states, we have only obtained symmetryadapted solutions even if the orbital optimization algorithm could, in principle, target symmetry-broken solutions. 15 It may be due to the lack of flexibility associated with the minimal basis. Indeed, we have shown, at the TpCCD level, that for larger molecules in larger basis sets one could also break the spatial symmetry to improve the description of excited states. 126 We expect analog symmetry-broken excited-state solutions for larger molecules in the case of VpCCD.
As shown by Henderson et al., 15 DOCI and TpCCD optimized orbitals are virtually indistinguishable in molecular systems. Here, we have observed that the VpCCD optimized orbitals are also virtually indistinguishable from the two other sets. Therefore, as expected, oo-TpCCD and oo-VpCCD energies are also highly similar, so that we only report oo-VpCCD energies in Fig. 4. The right panel of Fig. 4 evidences that the accuracy of oo-TpCCD and oo-VpCCD is similar to their MOM-TpCCD and MOM-VpCCD counterparts (see Fig. 3), at least in the weak correlation regime (always having DOCI as the reference results). However, in the strong correlation regime, the scenario is rather different. The orbital optimization at the correlated level allows them to strongly localize when the bond is stretched, hence the PECs exhibit the correct dissociation limits. As a direct consequence, the agreement between VpCCD (and TpCCD) and DOCI is improved at large R, as compared to MOM-VpCCD (and MOM-TpCCD). We can then conclude that, in the absence of strong correlation effects, state-specific RHF determinants provide robust and cheaper alternatives to determinants made of optimized orbitals at the correlated level. To further illustrate this, we provide the VpCCD optimized orbitals as well as the MOM orbitals in the supplementary material.

B. A strong correlation model: the ring H 4 molecule
We now turn our attention to another widely known model for strong correlation, namely the ring H 4 molecule where the four hydrogen atoms lie on a circle of diameter d = 6.569 a 0 . As represented in Fig. 5, varying the angle θ with respect to θ = 90°connects two equivalent D 2h rectangular geometries with non-degenerate molecular orbitals. At θ = 90°though, the D 4h square-planar geometry has strictly degenerate orbitals and strong multi-reference effects are at play. Therefore, giving an accurate description of this system as a function of θ has been shown to be a real challenge for CC methods. 14,25,52,[55][56][57]174,175 In the following, we restrict ourselves to the minimal STO-6G basis set in which the D 2h symmetry-adapted molecular orbitals are determined solely by symmetry. The resulting four molecular orbitals, ordered by ascending energy, have the symmetry a 1g , b 2u , b 3u , and b 1g . At θ = 90°, the b 2u and b 3u orbitals are degenerate and form a pair of orbitals of e g symmetry in the D 4h point group. Although one can choose to break spatial symmetry to gain flexibility, in a first stage, we restrict ourselves to the symmetry-adapted RHF molecular orbitals. In such situation, excited-state RHF wave functions correspond to non-Aufbau determinants built with this set of symmetry-pure orbitals, hence freeing ourselves from the orbital optimization issue to focus solely on the optimization of the cluster amplitudes. Because we deal with the seniority-zero subspace, the RHF determinants are made, by definition, of two doubly-occupied orbitals. For example, for the ground-state determinant at θ = 90°, the lowest-energy a 1g orbital and one of the doubly-degenerate b 2u and b 3u orbitals are doubly-occupied. Of course, in this case, the seniority-zero determinant is a poor approximation of the exact wave function as it tries to model an inherently multi-reference wave function with a single Slater determinant.
We start by looking at the description of the ground state using the symmetry-adapted orbitals. The DOCI (dashed lines), VpCCD (thick solid lines), and TpCCD (thin solid lines) energies are plotted in the bottom-left panel of Fig. 6. The configuration of the reference determinant is a 2 1g b 2 2u for θ < 90°a nd a 2 1g b 2 3u for θ > 90°. Hereafter, the electronic configuration of the RHF wave function is given for θ < 90°; the corresponding configuration for θ > 90°is obtained by simply swapping b 2u and b 3u . The first interesting fact is that TpCCD does not closely match DOCI for this system. On the other hand, VpCCD provides energies in fairly good agreement with DOCI. Therefore, the hermiticity property of VpCCD leads to a notable improvement over TpCCD. Yet, VpCCD exhibits a cusp at θ = 90°which is not present in DOCI. The derivative of the TpCCD PEC is also discontinuous at θ = 90°with an upside-down cusp compared to VpCCD. The comparison of the two pCCD variants and DOCI for the ground-state PEC provides similar insights to those reported in Ref. 52   state and the lowest-lying excited state form a conical intersection. This is a drawback of the HF approximation as FCI produces an avoided crossing (HF and FCI energies are given in the supplementary material). In the seniority-zero subspace, the avoided crossing between these two states is not present. Indeed, as shown in the bottom-left panel of Fig. 6, the DOCI dashed curves are smooth. Yet, they do not form an avoided crossing.
Then, we turn to the description of the excited states. The simplicity of the ring H 4 model in a minimal basis allows us to access the entire set of solutions using the Jenkins-Traub algorithm (see Sec. IV A). All the VpCCD solutions (with real cluster amplitudes) obtained with a ground-state RHF reference made of symmetry-adapted orbitals are represented in the center panel of Fig. 6, except for the quadruply-excited state which is plotted in the top-left panel. The convergence towards the ground state and the quadruply-excited state is fairly straightforward using the Newton-Raphson algorithm presented in Sec. II. However, the VpCCD solutions represented by the blue and cyan curves are the only two other solutions that we have been able to get for all θ values using the Newton-Raphson algorithm. In addition, we have obtained some parts of the three highest VpCCD PECs of Fig. 6, but the iterative algorithm was highly oscillatory and we have not been able to get any of these solutions for all values of θ.
The agreement between the VpCCD and DOCI excited states is less evident than for the linear H 4 model studied in Sec. IV A. The first important point to mention here is that there are more VpCCD than DOCI solutions, for all values of θ. More importantly, as we shall discuss below, it is challenging to tell which of these solutions are unphysical. We believe that three VpCCD solutions can be assigned to DOCI states with certainty: the ground state as well as the pink (quadruply-) and cyan (doubly-)excited states. However, even if the VpCCD solution corresponding to the quadruply-excited state (top-left panel) is attainable for all θ, it is a poor approximation to its DOCI counterpart, exhibiting a cusp and a local maxima at θ = 90°. Meanwhile, the two highest-lying DOCI doublyexcited states could correspond to some of the three VpCCD solutions (see top-right panel of Fig. 6). One could argue that the brown curve should be associated with the dashed brown curve, but for the two other VpCCD solutions it is hard to tell which one is unphysical. Finally, one can see in the bottom-left panel of Fig. 6 that the lowest-lying doubly-excited state could be associated with two VpCCD solutions. However, these solutions eventually disappear around θ = 80°and θ = 100°. Similarly to the spurious solutions in the linear H 4 model, it is possible that beyond this region the two solutions acquire complex-valued cluster amplitudes. Even for the part of the PECs where the solutions are real, the agreement with DOCI is quite poor (except around θ = 90°for the lower VpCCD solution), the PECs having the wrong topology. Alternatively, one could argue that the green VpCCD solution is associated with the lowest-energy DOCI excited state because it has the same topology as the first RHF excited state. Moreover, this solution exists for all geometry in contrast with the blue and yellow ones. We think that, at this stage, it would be arbitrary to assign a particular VpCCD solution to this DOCI state.
We now compare the previous set of VpCCD solutions with the ones obtained using non-Aufbau reference determinants made of the same set of symmetry-adapted orbitals. More specifically, we consider the lowest-lying RHF excited state of configuration a 2 1g b 2 3u , i.e., the other adiabatic state involved in the conical intersection with the a 2 1g b 2 2u RHF ground state (see the supplementary material). At θ = 90°, these two configurations become degenerate, and the choice of the e g orbital to occupy is arbitrary. Therefore, it seems logic to compare these two closely related references as function of θ. This may shed light on the meaning of the VpCCD solutions observed in Fig. 6. The four lowest VpCCD solutions of Fig. 6, i.e., the solutions obtained using the ground-state RHF determinant of configuration a 2 1g b 2 2u as reference, are also reported in Fig. 7 alongside the four lowest solutions obtained using the lowest-lying RHF excited state (a 2 1g b 2 3u ) as reference. One can see that three a 2 1g b 2 3u -VpCCD solutions (dashed) are connected with three a 2 1g b 2 2u -VpCCD solutions (solid) at θ = 90°, while the remaining pair of solutions coincide between θ = 80°and θ = 100°. For other θ values, the a 2 1g b 2 2u -VpCCD solution disappears.
As already mentioned earlier, the two lowest diabatic RHF states, a 2 1g b 2 2u and a 2 1g b 2 3u , intersect at θ = 90°. Hence, the lowest adiabatic RHF state has a cusp at θ = 90°. Cusps observed in ground-state CC PECs are often claimed to be unphysical. 23,52 However, it has been pointed out by Burton and Thom that these cusps are not unphysical but are consequences of the RHF reference used to construct the corresponding CC wave functions. 174 They further argued that these cusps indicate crossing of solutions at the CC level. This is indeed what we observe in Fig. 7 where the cusps on the VpCCD PECs are actually formed by two VpCCD solutions obtained with distinct reference RHF wave functions. In short, the inherent single-reference character of pCCD prevents it from correctly describing the FCI avoided crossing. On the other hand, a non-orthogonal CI (which is inherently multi-reference) between the two RHF states reproduces the correct shape of the PEC. 174 We should also mention that the projected CC method introduced by Qiu et al., in which one constructs a CCSD wave function on top of a projected HF reference, does not exhibit such a cusp. 175 As stated earlier, the D 2h molecular orbitals are fully determined by the spatial symmetry of the system. The corresponding set of symmetry-adapted (sa) orbitals is represented in Fig. 8 and ordered by ascending energies. However, one may wonder if there exists solutions associated with (spatial) symmetry-broken orbital sets. A stability analysis in the space of real RHF solutions reveals that the symmetry-pure groundstate RHF solution is a minimum with respect to occupiedvirtual rotations. [171][172][173] Thus, there is, a priori, no spatially symmetry-broken RHF solution lower in energy.
Next, we study the influence of orbital rotations on the VpCCD energy for the ground state of the ring H 4 model. The diagonalization of the orbital Hessian [see Eq. (16)] associated with this solution shows that this stationary point is an index-2 saddle point. Therefore, there is at least one additional state below the sa-VpCCD ground-state PEC. Indeed, by following the direction provided by the eigenvectors associated with the two distinct negative eigenvalues, we have been able to locate two additional VpCCD solutions. The first one, labeled as "sb 1 " (symmetry-broken) in Fig. 8, results from occupied-occupied and virtual-virtual rotations. Let us recall that, although the HF energy is invariant under occupied-occupied and virtual-virtual rotations, it is not the case for pCCD as the seniority-zero subspace depends on the orbital basis used to define it. 60 The direction associated with the second negative eigenvalue involves occupied-virtual rotations. Going downhill following the eigenvector associated with this eigenvalue leads to a different spatially symmetry-broken solution (see the "sb 2 " orbital set in Fig. 8). The point group of the electron density associated with each solution and the irreducible representation of each orbital are also given in the right panel of Fig. 8. A stability analysis of these two additional solutions shows that they are minima with respect to orbital rotations. Note that, for each set of symmetry-broken orbitals, four of the six possible reference determinants yield the same correlated energies. The two other references, 1a 2 1 1b 2 2 and 2a 2 1 2b 2 2 for the sb 1 set and 1a 2 g 2b 2 u and 1b 2 u 2a 2 g for the sb 2 set, yield energies close to the quadruply-excited state obtained with symmetry-adapted orbitals.
The left panel of Fig. 8 shows that the agreement between VpCCD and DOCI is much better when one considers the two sets of symmetry-broken orbitals. In addition, we emphasize that the sb 1 -DOCI/VpCCD PECs (which exhibit a cusp at θ = 90°) are fairly good approximations of the ground-state FCI PEC while containing only seniority-zero determinants. Likewise, the cuspless sb 2 -DOCI/VpCCD PECs are close in energy to the lowest-lying FCI excited-state one. The downside is that the corresponding wave functions do not possess the correct spatial symmetry. This is the famous Löwdin symmetry dilemma. [176][177][178] Moreover, for these two orbital sets, the TpCCD energies are in good agreement with DOCI (see the supplementary material). Yet, VpCCD is closer to DOCI than TpCCD as already seen in the linear H 4 case [see Sec. IV A]. The cusps exhibited by the sb 1 -DOCI/VpCCD PECs are also due to a crossing between two diabatic states. More specifically, it originates from a change of axis along which the D 2h symmetry breaks, leading to two different C 2v subgroups related by a rotation of π/2. Unfortunately, we have not been able to converge to the higher-lying symmetry-broken C 2v state. One would have noticed that we have not plotted the excitedstate TpCCD energies in Fig. 6. In fact, TpCCD suffers from the same issues related to additional solutions and their physical meaning. Similarly to VpCCD, projection on non-Aufbau references leads to moderate improvements. Yet, the TpCCD energy landscapes remain plagued by unphysical solutions. Consequently, it is hardly possible to assign a TpCCD solution to a given DOCI excited state, as discussed in the case of VpCCD.
Finally, we would like to mention that the improvement of VpCCD/TpCCD brought by state-specific reference wave functions is mitigated in comparison to the case of the linear H 4 molecule. Therefore, it seems that state-specific (MOM or oo-pCCD) references provide a very significant improvement for weak correlation, but does not help much in the presence of strong correlation. In such a case, if one is willing to sacrifice the spatial symmetry of the wave function, the description of the ground state (at least) can be improved. Symmetry-broken excited-state wave functions also exist at the VpCCD level but we have struggled to systematically converge towards these solutions. Hence, their performance still need to be properly assessed. This is left for future work.

V. CONCLUSION
Recently, there has been a renewed interest in singlereference methods for excited states in the context of Hartree-Fock, density-functional, and coupled-cluster theories. [103][104][105][106][107][108][109][110][111][112][113][114][115][116][117][118][119][120][121][122]126,143,169 This has been made possible thanks to the development of new algorithms specifically designed to target higher-energy solutions of these non-linear equations. These so-called non-standard solutions provide genuine alternatives to the usual linear response and equation-of-motion formalisms (which are naturally biased towards the reference ground state) for the determination of accurate excited-state energies in molecular systems. This is especially true for double excitations which are known to be difficult to model with the two latter formalisms. [91][92][93][94]99,[179][180][181] There is, therefore, a real need for a better understanding of the structure of the energy landscape associated with these methods. In this study, we have focused on the case of CC.
Due to the non-linearity of the CC equations, the topology of its energy landscape from which multiple solutions emerge is still far from being thoroughly understood. During the last decades though, several groups have been tackling this formidable problem. 101,103,126,143,182 In a recent study, 126 we have pursued along these lines by investigating the structure of the CC energy surface and the comparison between DOCI and TpCCD for excited states. More specifically, we have shown that the agreement which has been observed for ground-state energies [13][14][15][16]18,67 remains in the case of excited states only if one minimizes the TpCCD energy with respect to the orbital coefficients. In the present study, we have investigated the solution structure of the VCC method, a version of CC where the cluster amplitudes and the energy are determined variationally instead of the usual projective way. To the best of our knowledge, VCC excited states have never been investigated before.
Restricting ourselves to the case of pCCD (in which the cluster operator includes only pair excitations), we have looked at the VpCCD solution structure of two model systems, namely the linear and ring H 4 molecules, both in the minimal STO-6G basis. The former system has been used to investigate the influence of the orbital set on the VpCCD and TpCCD energy landscape. In contrast to TpCCD, VpCCD provides a much better approximation to the DOCI solution structure when one builds the reference determinant with ground-state RHF orbitals, at least in the weak correlation regime. When the correlation becomes strong, i.e., when the hydrogen chain is stretched, additional spurious solutions appear due to the truncation of the excitation operatorT , and VpCCD does not seem to improve with respect to TpCCD. In either regime, however, these excited-state solutions are hardly attainable using an iterative Newton-Raphson algorithm. Therefore, we replaced the ground-state RHF reference wave function by state-specific excited-state RHF references computed with MOM and targeted the corresponding VpCCD solution for each of them. We have observed that these state-specific references enlarge the basin of attraction of their associated solution, hence easing the convergence of the Newton-Raphson algorithm towards the targeted VpCCD solution. In addition, considering statespecific RHF orbitals greatly improves the TpCCD results for excited states. However, the difference between TpCCD and DOCI energies remains roughly one order of magnitude larger than the one between VpCCD and DOCI.
Then, we have turned our attention to the situation where the reference orbitals are optimized at the correlated level. In the weak correlation regime, the agreement between DOCI and the two variants of pCCD (TpCCD and VpCCD) is only slightly better than with MOM orbitals. However, in the strong correlation regime, the orbital optimization procedure allows the orbitals to localize further (while keeping their spatial symmetry), hence improving the accuracy of VpCCD and TpCCD (with respect to DOCI) at large internuclear separation. The take-home message of this first part is that TpCCD energies computed with state-specific RHF orbitals provide a good balance between robustness and computational cost to describe excited states, at least in the weak correlation regime. Of course, further studies on real molecules are required to assess the accuracy of these methods.
In a second stage, we have studied the ring H 4 molecule to investigate the influence of strong correlation on the energy landscape. We have seen that spurious VpCCD solutions, due to the truncation of the cluster operator, seem unavoidable in the presence of strong static correlation. Therefore, the description of excited states is much less accurate than in the weak correlation regime. Even worse, these spurious solutions prevent an unambiguous assignment of (some of) the excited states. This problem remains if one considers state-specific references at the VpCCD level. TpCCD suffers from the same issues, but in a more severe way. In addition, inspired by Burton and Thom, 174 we have investigated the physical origin of the cusps of the PEC at the VpCCD level. In agreement with Burton and Thom, 174 we have shown that, at the VpCCD level, these cusps are due to crossing of diabatic states obtained with distinct reference determinants.
Finally, we have investigated spatially symmetry-broken VpCCD solutions of the ring H 4 molecule. In a minimal basis set, the symmetry-adapted molecular orbitals are completely determined by the D 2h symmetry of the system. Yet, one can deliberately break this symmetry to relax the constraints imposed on the molecular orbitals. Doing so, we have shown that it is possible to locate two symmetry-broken VpCCD groundstate wave functions with energies in much better agreement with FCI, improving in the process the agreement between DOCI and VpCCD.

SUPPLEMENTARY MATERIAL
Included in the supplementary material are a standalone mathematica notebook gathering modules for the computational methods investigated here (HF, MOM, TpCCD, VpCCD, DOCI, and FCI), raw data for each figure, orbitals obtained at various levels of theory, and files containing the one-and two-electron integrals for the systems studied in the present manuscript.
hereafter we denote the CC wave function and its energy as |Ψ CC ≡ |Ψ and E VCC ≡ E. In this case, the derivative of the VCC energy functional with respect to the cluster amplitudes reads where we have used the assumptions that both the cluster amplitudes and orbital coefficients are real to simplify the second to last equality. Then, we differentiate the residual with respect to the cluster amplitudes to obtain the Jacobian matrix elements required to target saddle points. This yields the following formula: Ψ|a † a a † b a i a j a † c a † d a k a l |Ψ Ψ|Ψ . (A2)

Appendix B: VpCCD residual, Jacobian and density matrices
For the sake of completeness, we also provide the equations in the case of VpCCD. These can be easily derived from their VCCD analogs gathered in Appendix A by setting j σ = iσ, b σ = aσ, l σ = kσ, and d σ = cσ, where the indexes now refer to spatial orbitals andσ is the opposite spin of σ. Dropping the spin indexes, the matrix elements of the VpCCD residual are Ψ|a † a a † a a i a i a † c a † c a k a k |Ψ Ψ|Ψ . (B2) Following Ref. 15, in the case of pCCD, one can take advantage of the fact that |Ψ is a linear combination of seniority-zero determinants to compute the one-body density matrix γ [see Eq. (18a)] and the two-body density matrix Γ [see Eq. (18b)]. Indeed, many elements of γ and Γ are zero due to seniority considerations. Therefore, one only needs to compute the "diagonal" elements γ pp of the one-body density matrix, while for the two-body density matrix the only non-zero elements are Γ pp,qq , Γ pq,pq , and Γ pq,qp = − 1 2 Γ pq,pq .