Complex absorbing potentials within EOM-CC family of methods: Theory, implementation, and benchmarks

Metastable electronic states are important in diverse areas of science and technology ranging from high-energy applications (plasmas, attosecond and X-ray spectroscopies) to electron-molecule collisions (interstellar chemistry, radiolysis, DNA damage by slow electrons). These states (called resonances) can be accessed when molecules are excited above their ionization threshold, via electron attachment to closed-shell species, or by core ionization.

A concise and pedagogical introduction to the topic as well as references to earlier reviews can be found in Ref. 1. From the quantum mechanical point of view, resonance states belong to the continuum part of the spectrum and, therefore, their wave functions are not L²-integrable. Yet, their wave functions bear certain resemblance to the bound states within the interaction region (i.e., close to the nuclei). Using boundary conditions for the outgoing wave (Siegert or Gamow formalism, see Ref. 1), one arrives at the following form of the resonance wave function:

where the phase-isolated part ϕ_R(x) resembles a bound-state wave function in the interaction region and E_R and Γ (real and imaginary parts of the complex energy E = E_R − iΓ/2) determine resonance position and width. The latter is inversely proportional to the resonance lifetime. Thus, the resonances appear as solutions of the Schrödinger equation with complex energy.^1–4 One can arrive at the same concept of complex energy via a completely different formalism (Feshbach approach) based on a separation of the Hamiltonian into coupled bound and continuum parts; in this approach, the resonance is described as a bound state coupled to the continuum, and the complex energy emerges from solving a non-Hermitian eigenproblem with an effective Hamiltonian.⁵ 

One can avoid the inconveniences of working with continuum functions or fiddling with boundary conditions by reformulating the problem using complex variables.^2–4 The most rigorous approach is the complex-scaling formalism^2,6 in which all coordinates are scaled by a complex number e^−iθ; however, practical applications of this method are limited by its extreme sensitivity to the one-electron basis set^7–11 as well as conceptual difficulties regarding the separation of nuclear and electronic motions of the scaled Hamiltonian.^12–15 

These problems are avoided in an alternative approach in which the original (non-scaled) Hamiltonian is augmented by a complex potential

$−iηŴ$ devised to absorb the diverging tail of the resonance wave function.^16–18 In the complete basis set limit, these complex absorbing potential (CAP) methods yield exact resonance positions and widths in the limit of zero CAP strength η.¹⁹ It can be shown that CAP methods are related to exterior complex scaling methods.^20,21

The use of CAPs in practical calculations is complicated by possible reflections leading to false resonances, sensitivity of the results to the form of the CAP W, and a strong basis set dependence.^19,22 Furthermore, one has to determine an optimal value for η, which is usually achieved by calculating trajectories E(η) and requiring |η dE/dη| = min.

In terms more familiar to electronic structure practitioners, reflections can be described as perturbations of the resonance wave functions and, consequently, energies caused by a finite-strength CAP. Finite basis sets give rise to additional reflections. Several reflection-free CAPs have been proposed;^19,23 such CAPs are either energy-dependent or non-local.

In our previous paper,²⁴ we introduced a simple density-matrix based correction to the energy that removes the perturbation due to the CAP. The correction was derived based on energy decomposition analysis and response theory. Our starting point was the observation that the energy's dependence on η becomes roughly linear beyond some critical value of η. By analyzing the response equations, we also proposed an alternative criterion for finding an optimal value for η. Physically, our approach is grounded in the behavior of the resonance wave function and, ultimately, the one-particle density matrix. It was shown²⁴ that when the CAP is sufficiently strong, both real and imaginary parts of the density become near-stationary indicating that the resonance is stabilized. Then the perturbation to the resonance position by the CAP can be eliminated by subtracting the term η Tr[γW] from the energy. The optimal η is found by considering the de-perturbed resonance energies; moreover, we argued that η_opt is not the same for real and imaginary parts.²⁴ Preliminary benchmarks illustrated that this approach results in a computationally more robust scheme in which the dependence on the onset of the CAP is significantly reduced compared to the straightforward application of a CAP along with the original energy-based criterion for finding the optimal η, as was done in most CAP applications.^25–27 We note that Moiseyev et al.^28,29 also observed a linear energy dependence on η beyond some critical value of η and proposed the use of Padé approximants to extrapolate the energy of the stabilized resonance to the zero η limit. If the energy depends strictly linearly on η, their approach yields results identical to our first-order correction.

One of the difficulties of understanding the capabilities and limitations of different approaches is that a method that has shown excellent performance for a small model problem may fail when applied to a realistic system. In the context of electronic structure, the results of calculations of resonances will also be affected by the quality of standard approximations such as the incompleteness of one- and many-electron basis sets.³⁰ Thus, it is important to test different methods for meta-stable states within robust and accurate ab initio approaches. For bound states, the coupled-cluster (CC) and equation-of-motion (EOM) hierarchies of methods^31–35 provide a reliable and predictive set of theoretical model chemistries.³⁶ These methods can be systematically improved to approach the exact solution, are size-extensive (size-intensive for excitation energies), describe dynamical and non-dynamical correlation in one computational step, and do not involve system-dependent parameterization. The CC hierarchy of methods works best for wave functions dominated by a single Slater determinant, however, the EOM-CC approach extends this single-reference formalism to tackle various open-shell and multi-configurational cases.^33,37

In EOM-CC, the target-state wave function is described by an excitation operator

$R̂$ acting on the reference-state CC wave function

where

$|0⟩$ is the reference Slater determinant – usually satisfying the Hartree-Fock (HF) equations – and

$\hat{T}$

$T̂$ is the coupled-cluster operator. Different choices of

$\hat{R}$

$R̂$ provide access to different target states, e.g., in EOM-EE-CC

$\hat{R}$

$R̂$ is electron and spin-conserving thus enabling the description of various excited states. Open-shell electron-attached states (such as temporary anions) can be described by EOM-EA-CC in which the reference state is again a well-behaved closed shell state and the operator

$\hat{R}$

$R̂$ changes the number of electrons. Likewise, ionized states can be described by EOM-IP-CC with the operator

$\hat{R}$

$R̂$ removing an electron. Thus, EOM-CC is a natural choice for extending the excited-state methodology to resonances via complex scaling and CAP approaches.

Recently, we presented an implementation of complex-scaled EOM-CC singles and doubles (EOM-CCSD) methods and illustrated their performance by considering several atomic systems (He, H⁻, Be). Here, we present an implementation of CAPs within the EOM-CC family of methods. Our main focus is on the EOM-EA-CC variant; however, our implementation also includes EOM-EE-CC. While limited implementations of CAPs within EOM-CC have been reported before (e.g., Refs. 25 and 26), this work presents the first formally complete and production-level implementation of the method.

The main focus of the paper is on investigating the numeric performance of the method and the sensitivity of the results towards the CAP parameters and the choice of basis set. Our goal is to develop a black-box approach that could be calibrated and then applied to the calculation of resonances without any prior knowledge of the system, as advocated by John Pople.³⁶ In particular, we want to avoid the system-dependent optimization of basis sets and the CAP's shape and onset. Thus, rather than aiming at results converged with respect to all computational parameters individually for each system, we wish to establish a uniform protocol that can be applied to any system and can be characterized by error bars estimated from prior calibration studies, as routinely performed in electronic structure calculations.³⁰ 

We note that the validation of the accuracy of computed resonance lifetimes in molecular systems is difficult.³⁸ A complete theoretical description should involve coupled electronic and nuclear dynamics; this is beyond the scope of the present paper, where we only compute the lifetime of the resonance state at a fixed molecular geometry. This is appropriate for resonances whose lifetimes are short relative to nuclear motions, or when nuclear motions do not strongly affect the computed Γ values (Condon-like approximation). Thus, our focus is on the comparison with other theoretical studies and the robustness of the results with respect to the one-electron basis set as well as variations of the CAP parameters.

In this context, we add that the sensitivity of the results to the one-electron basis set is of a fundamentally different origin in CAP calculations as compared to complex scaling. In the latter case, the basis should be sufficiently flexible to describe the resonance wave function at different values of the scaling angle, whereas in the former case, one simply needs to supply a basis set of a sufficient spatial extent to represent a given CAP and a stabilized resonance wave function. This implies that the diffuseness of the basis must be coordinated with the CAP onset, e.g., in a compact basis, the CAP onset should be smaller, otherwise, the calculation will be blind to the CAP. Thus, although the basis-set dependence is a nuisance, its simpler nature in CAP calculations suggests that a solution can be found.

Originally, CAP methods were introduced to study shape resonances. Since the decay of Feshbach resonances is a two-electron process (and the CAP is a one-electron operator), one may expect difficulties in describing Feshbach resonances within the CAP formalism. Moiseyev et al.^39–41 showed that additional steps need to be taken for the construction of reflection-free CAPs in order to reliably calculate Feshbach-type resonances. The present paper focuses solely on understanding CAPs in the context of molecular shape resonances.

The article is structured as follows: Secs. II and III present the formalism of the CAP-augmented EOM-CC method and our implementation. In Sec. IV, we put forward a protocol for determining the resonance positions and lifetimes and investigate its robustness towards the choice of the one-electron basis set and the CAP's onset. In Sec. V, we subsequently apply our new scheme to a variety of molecular resonance states and compare the results to those obtained from experiment as well as using other theoretical approaches. Sec. VI provides concluding remarks.

The basic idea of the CAP method^16–19 is the addition of an artificial complex potential to the original Hamiltonian

where W aims to absorb the outgoing electron and η controls its strength. As in complex scaling,^2–4,11,42 the addition of the CAP results in a non-Hermitian complex symmetric operator H(η)¹⁶ converting resonances into square integrable (L²) wave functions. In our calculations, we choose the CAP as a quadratic potential with an unaffected region of cuboid (i.e., box) shape

with r_α denoting the three Cartesian coordinates (α = x, y, z). Thus, the CAP is controlled by 4 parameters: 3 parameters for the onset in each direction (⁠

$rx0$⁠,

$r_y^0$

$ry0$⁠,

$r_z^0$

$rz0$⁠) and the strength η. In principle, the CAP strength is unbound (η ∈ [0, ∞)), but should be chosen such that the effect is large enough to absorb the wave function over a certain range, but not too large to prevent excessive perturbation of the wave function and the resonance energy.¹⁹ 

In the complete one-electron basis set, the exact position of the resonance in the complex plane can be obtained as lim_{η → 0}E(η).¹⁶ That is, an infinitesimally weak CAP, which is represented exactly (and, therefore, goes to infinity at large r), is sufficient to stabilize the resonance without perturbing it. Working with finite Gaussian basis sets requires one to perform series of calculations for different η in order to find an optimal value of the strength parameter η_opt along the η-trajectory and the corresponding value of the resonance energy E(η_opt). A commonly used criterion for determining the optimal value of the strength parameter η is finding the minimum of the logarithmic velocity^16,19

Unfortunately, the position and width of the resonance computed using this criterion are very sensitive to the CAP onset^17,24,27; thus, this approach does not provide a black-box scheme. In our recent paper,²⁴ the first-order deperturbative correction to the raw zeroth-order resonance energies E^R and E^I was introduced as

with γ(η) as the one-particle density matrix. The correction is based on perturbation theory; it removes the explicit dependence on the CAP from the computed resonance energies. As was illustrated in Ref. 24, the corrected energies, U^R(η) and U^I(η), exhibit nearly constant behavior for large η (that is, past the stabilization point, when the resonance wave function does not change much anymore); furthermore, the corrected trajectories computed using different CAP onsets become much more similar in the asymptotic region as opposed to the uncorrected trajectories, E^R(η) and E^I(η).

By looking separately at the real and imaginary parts of the deperturbed energy (U^R and U^I) as a function of η, we showed that the energy becomes near stationary at certain values of optimal strength (⁠

$η opt R$ and

$\eta _{\rm opt}^I$

$η opt I$⁠) giving the position and lifetime of the resonance. Our results showed that this recipe leads to values for the resonance position and lifetime that are less sensitive to the CAP onset and thus more robust than the zeroth-order values E^R and E^I.²⁴ 

In our scheme, the CAP is introduced at the HF level of theory, where we obtain a set of complex molecular orbitals (MOs) as the solution for a given strength η. Hence, a complex Koopmans' theorem holds for the virtual orbitals, i.e., they can be interpreted as zeroth-order approximations of the resonance and rotated continuum states.

As the next step, we solve the CCSD equations for the reference state^{31,32,43–45} using H(η)

with Φ_μ denoting the excited determinants. The resulting amplitudes t^η are also complex.

To compute electronically excited and electron-attached resonance states, we use EOM-EE-CCSD and EOM-EA-CCSD^33,46–50 methods that provide accurate and predictive descriptions for such target states. The wave function of the resonance state is found by solving a non-Hermitian eigenvalue problem for the right eigenvectors

which yields a set of complex amplitudes R^η and complex excitation energies Ω^η. The latter are the raw, η-dependent resonance energies (which are equal to the difference between the total energy of the excited/attached EOM-CCSD state and the reference CCSD energy for a given η).

We note that for moderate CAP strengths and CAP onsets comparable with the spatial extent of the electron density of the reference state, the perturbation to the reference CCSD energy is small (10⁻⁵ a.u.) in contrast to complex-scaled calculations.¹¹ To compute the first-order correction to the raw resonance energies, the one-electron density matrix needs to be calculated. We employ an unrelaxed one-electron EOM-CCSD density matrix containing no amplitude- or orbital-response terms⁵¹ 

where L^η and R^η are the left and right EOM-CCSD eigenvectors, respectively.

Because of the non-Hermitian nature of

$H¯$ left eigenvectors have to be computed and biorthogonalized against the right eigenvectors in order to compute the density matrix

where i and j denote different electronic states. Once the density matrix is computed, the energy of the resonance state is corrected according to Eqs. (7) and (8).

The equations for CAP-CCSD and CAP-EOM-EE/EA-CCSD are identical to the original CCSD and EOM-EE/EA-CCSD equations except that all quantities such as the Fock matrix, the two-electron integrals, the MO matrix C, the T and R/L amplitudes are now complex and η-dependent. There is no need to add the CAP explicitly to the CCSD or EOM-CCSD equations since it is already included at the HF level.

In our paper on complex-scaled EOM-CC,¹¹ we considered several variants of implementation, including one in which the HF and CCSD equations for the reference state were solved for the unscaled Hamiltonian and the scaling was introduced only at the EOM-CC level. This required significant reformulation of the EOM-CC equations. By analogy, one may also consider an implementation of CAP-EOM-CC in which the CAP is introduced only at the EOM-CC level; this will be the subject of future work.

Due to the CAP, the Hamiltonian becomes non-Hermitian and complex symmetric, which necessitates using a different metric, the so-called complex symmetric scalar product (c-product),^16,42,52,53 such that the variational principle is maintained

The difference to the regular scalar product is that the bra-vector is not complex conjugated. Mathematically, the c-product is a pseudoscalar product which does not induce a valid metric norm.^42,52 However, one can still define the c-norm (f|f) which is, contrary to the regular norm, complex in general and might become zero for a non-zero function f (“self-orthogonality”)

where ⟨|⟩ is a regular scalar product which is equivalent to the c-product for real functions. Thus, the normalization of all vectors (for example, left and right EOM-CC eigenvectors) is done by multiplying by a complex number⁴² 

resulting in

$f̂$ being a normalized vector. As mentioned above, use of the c-norm may lead to “self-orthonormality” ((f|f) = 0 for f ≠ 0), but this has not been observed in practice. Orthogonality is defined for the c-product in the same manner as for the scalar product as (f|g) = 0 (c-orthogonality).

The suite of CAP-EOM-CC methods has been implemented in the Q-Chem electronic structure package^54,55 and was released in version 4.2. All complex CCSD and EOM-EE/EA-CCSD equations were implemented using the libtensor library⁵⁶ for high-performance tensor operations.

The calculations begin by solving the CAP-augmented restricted HF equations (CAP-RHF). CAP-RHF has been implemented as an extension of regular RHF using the object-oriented SCF library SCFman⁵⁷ in Q-Chem that employs the Armadillo linear algebra library⁵⁸ for matrix computations. We add that an adaptation of our implementation for CAP-UHF will be straightforward. The CAP is introduced as an additional term in the regular Fock matrix

The molecular orbitals must satisfy the following orthonormalization condition in the c-product metric:

where S is the overlap matrix in the atomic orbital (AO) basis and (C^η)^T is transposed but not conjugated. Since the augmented Fock matrix F^η is non-Hermitian, the orbitals obtained using standard linear algebra routines for the diagonalization of general matrices are not normalized. In order to satisfy the orthonormality condition [Eq. (18)], the MOs are orthogonalized by using a modified Gram-Schmidt procedure with projections calculated using the c-product.

The Fock matrix

$Fμνη$ and the two-electron integrals (μν|λσ) are transformed into the MO basis by applying the complex orbital transformation matrix,

$C_{\mu p} ^{\eta }$

$Cμpη$⁠; thus, these quantities become complex in the MO basis. The calculations proceed by solving the CCSD amplitude equations using a DIIS procedure⁵⁹ adapted for complex algebra with c-product. Once the complex t-amplitudes are converged, we find the excited-state energies and right eigenvectors by using Davidson's procedure⁶⁰ generalized for non-Hermitian complex matrices. Note that the original

$\bar{H}$

$H¯$ matrix is also non-Hermitian but real, thus, one only needs to modify the procedure to make it work with complex quantities and the c-product. We observe that for large values of η the convergence of Davidson's procedure is sometimes problematic, likely due to more pronounced non-Hermiticity. However, we were always able to converge a reasonable number of roots (2-10) by tweaking the parameters of Davidson's procedure such as subspace size, residual inclusion threshold, etc. Since the one-electron density matrix is needed for the calculation of the first-order correction to the energies, we also solve for the left eigenvectors using Davidson's procedure, as well as for left and right eigenvectors together to ensure their c-biorthogonality [Eq. (13)].

The CAP is evaluated in the AO basis through numerical quadrature using a Becke-type grid⁶¹ of (99, 590) points (99 radial points and 590 angular points per radial point). Currently, we have implemented a shifted quadratic potential (r − r₀)² for a rectangular cuboid [Eqs. (4) and (5)], but our implementation allows for an easy extension of the shape of the unaffected region as well as the type of potential, e.g., higher order monomials (r − r₀)⁴, (r − r₀)⁶, etc. Our implementation also includes the optional addition of a real potential to the CAP, as was advocated in Ref. 41.

Relative to the conventional CCSD and EOM-CCSD methods, the addition of the CAP does not change the scaling of the computational cost [O(N⁶) for CCSD and EOM-EE-CCSD, O(N⁵) for EOM-EA-CCSD] or the memory requirements [O(N⁴)]. However, because we need to work with complex numbers, the computational cost increases roughly by a factor of 4 and the storage requirements increase by a factor of 2. Furthermore, computation of the first-order energy correction requires the left eigenvectors, which increases the computation time relative to EOM-CCSD energy calculations. Another possible issue arises when the resonance state is lying high in energy, so that the standard Davidson procedure has to find all lower roots that might require a lot of iterations to converge. To solve this issue, we have implemented iterative solvers for interior eigenstates for conventional EOM-CCSD methods. Implementation of the interior eigenvalue solvers for CAP-EOM-CCSD methods is a subject of future work. Finally, the necessity to compute η-trajectories requires to run calculations for different values of η, but these calculations can be performed independently and can therefore be run in parallel.

The necessity to find optimal values for the strength and the onset of the CAP as well as a pronounced basis set dependence of the results have precluded routine applications of CAP-based methods so far. Hence, investigating the numeric performance of CAP-EOM-CCSD, in particular with respect to the two aforementioned issues, is crucial for it to become a useful tool for studying resonances.

As for the one-electron basis set, it has been established that the straightforward application of standard basis sets yields poor results. Additional diffuse functions need to be incorporated for two reasons: (i) to obtain a sufficiently good basis-set representation of the CAP and (ii) to describe the outgoing electron correctly.^16,19,62 Owing to these requirements, CAP-based computations were often carried out using non-standard basis sets.^{26,27,62–67} While such approaches are able to provide results that agree with experiment for some resonance states, a treatment based on standard basis sets not involving any optimization procedure would be superior as it is of black-box type, has predictive power, and is also computationally less demanding. Since the shape of the resonance wave function is similar to a bound-state wave function in the interaction region, it should be possible to lessen the basis-set dependence to the degree observed in regular EOM-CCSD calculations.

Concerning the choice of the CAP onset, one has to realize that the artificial nature of the CAP implies that it is impossible to deduce from basic physical laws a universally applicable procedure for finding optimal parameters for the CAP strength, onset, and shape. While this is unsatisfying from a formal point of view, a pragmatic approach is to mitigate the dependence of the physically meaningful results on the artificial parameters as much as possible. Along these lines, we introduced a first-order correction²⁴ that was shown to desensitize resonance positions and widths to the choice of onset parameters by removing the perturbation due to the CAP. However, as this dependence cannot be removed completely, one has to develop a protocol for the unique and system-independent choice of the CAP onset to make CAP-EOM-CCSD applicable in a routine manner.

In the following, we examine the π^* shape resonances of CO⁻ and

$C2H4−$ using different basis sets and CAP onsets. Both states arise from adding an electron to the lowest unoccupied valence MO (LUMO) of the respective neutral molecule. Bond lengths and angles were chosen as R(CO) = 2.1316 a.u. for carbon monoxide and R(CC) = 2.5303 a.u., R(CH) = 2.0522 a.u., and ∠(CCH) = 121.2° for ethylene. All electrons were active in the correlated calculations. The CAP strength η was varied with a step size between 0.0001 a.u. and 0.001 a.u. Optimal values for η were determined according to the criterion from Eq. (6) as well as using the procedure outlined in Ref. 24. The respective results are referred to as zeroth-order and first-order in all tables and in the discussion below.

The basis sets used in our calculations were derived from the aug-cc-pVXZ (X = D, T, Q, 5) series⁶⁸ through augmentation by additional even-tempered basis functions. In all cases, the exponents for the first additional basis functions were obtained as one half of the exponent of the most diffuse basis function with the same angular momentum in the parent aug-cc-pVXZ basis set. The exponents for the remaining additional basis functions were calculated as one half of the exponent of the preceding function.

We explored two different series of basis sets, namely, one where we augmented the basis sets for all atoms except hydrogen (denoted as (A) below) and one where we placed only one set of diffuse functions with averaged exponents in the center of the molecule (denoted as (C) below). Since the second approach is computationally less demanding, it is especially preferable when targeting larger systems. To ensure that this strategy does not give rise to artifacts, we computed EOM-EE-CCSD excitation energies for a number of bound excited states of CO and C₂H₄ using the aug-cc-pVTZ+3s3p3d(C) and aug-cc-pVTZ+3s3p3d(A) bases. The results are reported in Table I together with the corresponding values for ⟨R²⟩, which are helpful in distinguishing valence states from Rydberg states. As apparent from Table I excitation energies computed with the two basis sets differ by not more than 0.04 eV for the very diffuse 4 ¹A₁ state of CO and by just 0.001 eV for all other states. This shows that placing the additional diffuse functions in the center does not lead to inferior results for the bound excited states and thus suggests the application of this scheme to resonances. We note that a similar scheme has been employed before in the context of stabilization techniques.^69,70

As for the choice of the CAP onset, we employed the square roots of the expectation values ⟨α²⟩ (α = x, y, z) for the ground states calculated at the CCSD level of theory as a starting point and considered the impact of small variations. The values used are

$rx0=ry0=2.76$ a.u. and

$r^0_z = 4.97$

$rz0=4.97$ a.u. for CO and

$r^0_x = 7.10$

$rx0=7.10$ a.u.,

$r^0_y = 4.65$

$ry0=4.65$ a.u., and

$r^0_z = 3.40$

$rz0=3.40$ a.u. for C₂H₄. The orientation of the molecules was as follows: For CO, the z-axis formed the molecular axis, whereas the C₂H₄ molecule was placed in the xy-plane with the CC bond oriented along the x-axis.

Table II compiles resonance positions E_R and widths Γ for the ²Π resonance of CO⁻ and the ²B_2g resonance of

$C2H4−$ obtained using the aug-cc-pVTZ+3s3p3d(C) basis and different CAP onsets. We started with the aforementioned values for

$r^0_\alpha$

$rα0$ based on the spatial extent of the ground-state wave function and then varied

$r^0_x$

$rx0$⁠,

$r^0_y$

$ry0$⁠, and

$r^0_z$

$rz0$ independently. In addition to E_R and Γ, we report optimal CAP strengths as well as values for the norm of the CAP in the AO representation. As expected, the representation of the CAP becomes more complete for smaller values of

$r^0_\alpha$

$rα0$⁠. In addition, the results for ‖W‖, the Frobenius norm of the CAP matrix, show that the onset parameters are not all of the same importance: Consistent with the π^* character of the resonance states,

$r^0_z$

$rz0$ makes the largest impact for

${\rm C}_2{\rm H}_4^-$

$C2H4−$ and

$r^0_x=r^0_y$

$rx0=ry0$ for CO⁻. This trend is also reflected in all values for η_opt, E_R, and Γ: Varying the pivotal onset parameter by ±0.5 a.u. shifts zeroth-order values for E_R by up to 0.028 eV and zeroth-order values for Γ by up to 0.048 eV, while the impact of the remaining onset parameters is roughly one order of magnitude smaller. Thus, in the remaining discussion we will focus on the onset parameter with the most pronounced influence.

Table II shows that both zeroth-order and first-order resonance positions and widths become smaller when increasing

$rα0$⁠, but we emphasize that these fluctuations are mitigated, especially for the width, when considering first-order results: Here, E_R and Γ are both shifted by at most 0.025 eV upon variation of

$r^0_\alpha$

$rα0$⁠. Also, it is apparent from Table II that smaller values for the CAP onset lead to smaller η_opt and that the first-order correction always entails larger values for η_opt. However, the relevance of the latter trends is debatable as the CAP strength η is not a physically meaningful quantity. One can argue along the same lines regarding the quantity η · dE/dη that needs to be minimized to find the optimal value for η when using the conventional criterion from Eq. (6). As can be seen from Table II, trends in η · dE/dη are weakly pronounced and not uniform.

To analyze the convergence of resonance positions and widths with respect to the addition of diffuse basis functions, we studied the π^* resonances of CO⁻ and

$C2H4−$ using different augmentations. The results are summarized in Table III. The crucial role of the angular momentum of the additional basis functions becomes clear at the first glance: In the case of CO⁻, a jump of almost 0.5 eV is observed for the zeroth-order resonance position when going from aug-cc-pVTZ+3s(C) to aug-cc-pVTZ+3s3p(C), while an additional augmentation by three sets of d-functions leads to a change of 0.06 eV. Three sets of f-functions on top of aug-cc-pVTZ+3s3p3d(C) shift E_R by just 0.007 eV. For the resonance width, d-functions play a more important role: Going from the 3s(C) to the 3s3p(C) augmentation changes Γ by 0.08 eV and the next step to 3s3p3d(C) changes Γ by 0.10 eV, but the value for 3s3p3d3f(C) differs from that for the preceding augmentation by just 0.004 eV.

For

$C2H4−$⁠, the changes are in general of similar magnitude as for CO⁻, but the big jump is observed when adding d-functions for both E_R and Γ. We add that similar trends are found for the first-order E_R and Γ of both resonance states. Also, we note that the use of basis sets with an augmentation including just s-functions (for CO⁻) or just s- and p-functions (for

${\rm C}_2{\rm H}_4^-$

$C2H4−$⁠) entails much larger η_opt values and sizable differences between zeroth-order and first-order values. Finally, the importance of the angular momentum of the diffuse basis functions is also reflected in the E(η) trajectories for

${\rm C}_2{\rm H}_4^-$

$C2H4−$ displayed in Figure 1. Their shape is altered considerably when d-functions are added, but is very insensitive to the addition of p-functions or f-functions. The role of angular momentum can be easily related to the spatial symmetry of the resonance states thus allowing us to choose the augmentation scheme based on symmetry considerations prior to the actual computations. In addition, we note that the values of ‖W‖ show that basis functions with angular momentum higher than ℓ = 2 do not significantly improve the basis-set representation of the CAP.

We also investigated the effect of adding more than three additional diffuse s, p, and d-functions. The corresponding results in Table III show that, while values for ‖W‖ become considerably larger indicating a more complete basis-set representation of W, the impact on resonance positions and widths does not exceed 0.035 eV except for one case: For

$C2H4−$⁠, the zeroth-order E_R and Γ calculated using the augmentation schemes 3s3p3d(C) and 6s6p6d(C) differ by more than 0.1 eV. However, this discrepancy disappears when the first-order correction is applied. These results suggest that an accurate basis-set representation of the CAP near the interaction region is crucial for obtaining correct resonance positions and widths, whereas regions further away do not need to be covered by the basis set. Furthermore, since an increased dependence of E_R and Γ on the CAP onset is found in some cases, we conclude that it is neither necessary nor advisable to employ more than three additional sets of diffuse basis functions with the required angular momentum.

Table III also reports the results from calculations with the aug-cc-pVTZ+3s3p(A) and aug-cc-pVTZ+3s3p3d(A) bases. Contrary to what we observed for bound states (cf. Table I), the differences between values for E_R and Γ obtained with the two augmentation schemes “C” and “A” are not negligible. Zeroth-order values differ by up to 0.23 eV and first-order values still by up to 0.15 eV. In one case, namely,

$C2H4−$/aug-cc-pVTZ+3s3p, discrepancies of 0.4 eV are found; this is probably related to the poor performance of the 3s3p augmentation scheme for

${\rm C}_2{\rm H}_4^-$

$C2H4−$ discussed earlier. We consider results obtained with the scheme “C” superior for several reasons: From the trajectories shown in Figure 2, one can see that the first-order quantities U^R and U^I enter the region of near-stationarity for smaller η, i.e., the resonance wave function shows faster convergence with respect to η, which is reflected in smaller η_opt values obtained in calculations using the augmentation scheme “C.” Also, an increased dependence on the CAP onset is found in some cases when using scheme “A.”

Besides the impact of additional diffuse functions, variations in the valence basis set also influence E_R and Γ. This is illustrated by Table IV, which reports values for resonance positions and widths of the ²Π resonance of CO⁻ and the ²B_2g resonance of

$C2H4−$ computed using the aug-cc-pVXZ (X = D, T, Q, 5) bases augmented according to the 3s3p3d(C) scheme. Concerning the resonance position, one can see that for both zeroth-order and first-order values the basis-set dependence is more pronounced than for excitation energies corresponding to bound states of the neutral molecules. The position of the resonance state in CO⁻ changes by 0.06 eV when going from aug-cc-pVQZ to aug-cc-pV5Z, whereas the largest shift observed for a bound state is less than 0.03 eV. We also note that the positions of the resonance states become smaller with increasing basis-set size, while the opposite is true for the excitation energies of the bound Rydberg states.

Concerning the resonance width, trends are less clear. For CO⁻, both zeroth-order and first-order values show a non-monotonic behavior with respect to the basis-set size and no convergence is observed. The magnitude of the variations in Γ is however comparable to that in E_R. For

$C2H4−$ in contrast, the dependence of Γ on the basis-set size is much less pronounced and convergence seems to be reached. We also see that the dependence of both E_R and Γ on the CAP onset is somewhat mitigated when increasing the size of the valence basis set: For the aug-cc-pVDZ basis set, a decrease of

$r^0_x$

$rx0$ and

$r^0_y$

$ry0$ by 0.5 a.u. increases E_R and Γ of the ²Π resonance of CO⁻ by 0.032 eV and 0.073 eV, respectively, whereas the same decrease leads to changes of just 0.025 eV and 0.026 eV when using the aug-cc-pVQZ basis. This should be contrasted with the opposite impact of additional diffuse functions discussed before. Table IV also shows that the values for η_opt decrease for larger basis sets, which is in line with that η_opt should be zero in the complete basis set limit.¹⁶ 

To gain further insight into the dependence of E_R and Γ on the size of the valence basis set, we performed an energy decomposition analysis for the ²Π resonance of CO⁻ and the ²B_2g resonance of

$C2H4−$ based on the following partition of the electronic Hamiltonian:

where the CAP is considered as a part of F_pq and

$⟨pq‖rs⟩$ stands for the two-electron integrals in MO basis. The expectation value of the one-electron part is then interpreted as one-electron energy, whereas the expectation value of the remainder represents the contribution from the EOM-EA-CCSD two-particle density matrix. The results are compiled in Table V. For the real part of the energy, this illustrates that the one-electron part converges significantly faster to the complete basis-set limit than the two-electron part, which suggests that the overall slow convergence of the total energy is mainly driven by an incomplete treatment of electron correlation. In contrast, for the imaginary part of the energy, the one-electron and two-electron parts seem to diverge in opposite directions with a basis set increase. This holds true for both molecules, but whereas the trends roughly cancel out for

${\rm C}_2{\rm H}_4^-$

$C2H4−$⁠, this is not the case for CO⁻ leading to a seemingly different behavior for the overall resonance width of the two systems.

One might be tempted to relate the basis set dependence of Γ to an insufficient description of the interaction of the resonance state with the continuum, but we point out that the addition of further diffuse functions has only little impact on E_R and Γ (cf. Sec. IV C), which suggests the opposite. In sum, we feel that the behavior of Γ requires further investigation in order to develop a scheme for the extrapolation to the complete basis set limit, but such an extension is beyond the scope of the present article. We point out, however, that the variations in Γ observed for CO⁻ do not exceed 0.15 eV so that reliable computations are still possible based on our current approach.

In this section, we will report resonance positions and widths for a number of shape resonances of small to medium-size molecules and compare the results from our CAP-EOM-EA-CCSD scheme to those obtained using other theoretical approaches or through experiment. Systems included in this study are

$N2−$⁠, CO⁻,

${\rm C}_2{\rm H}_2^-$

$C2H2−$⁠,

${\rm C}_2{\rm H}_4^-$

$C2H4−$⁠, CH₂O⁻,

${\rm CO}_2^-$

$CO 2−$⁠, and C₄

${\rm H}_6^-$

$H6−$⁠. All these resonance states except for the last one are derived by electron attachment to the π^* LUMO of the corresponding neutral molecules. C₄

${\rm H}_6^-$

$H6−$ (1,3-butadiene) is a special case as its π system extends over more than a single double bond, which results in two low-lying resonance states.

We add that only the real part of the resonance wave function has a well defined single-attachment character, while the imaginary part has a considerably different form. Most often, it exhibits multiconfigurational character represented by several single attachments to very diffuse orbitals. Also, its dependence on the CAP strength is more pronounced than that of the real part. A detailed investigation of this phenomenon will be conducted in future work.

In all calculations, we employed the scheme developed in Sec. IV, i.e., we chose the CAP onset based on the spatial extent of the ground-state wave function and used the aug-cc-pVXZ+3s3p3d(C) bases. Computational details are compiled in Table VI.

The scattering of slow electrons by the N₂ molecule has been studied experimentally several times^71–78 so that the ²Π_g resonance of

$N2−$ is rather well characterized. Consequently, this resonance has served as a testing ground for numerous theoretical approaches including stabilization techniques,^69,70,79,80 methods based on complex scaling,^10,81,82 CAP-based schemes,^{16,26,27,62–66} as well as other approaches.^83–85 Among various aspects, the impact of electron correlation on the resonance position and width⁶² as well as their basis-set dependence^26,62,64 have been investigated in detail. In addition, the resonance wave function has been studied over a wide range of different bond lengths and adiabatic excitation energies have been determined.^10,80,86,87 The potential interplay of the ²Π_g ground state of

${\rm N}_2^-$

$N2−$ with other resonance states has been also investigated.^10,88

CAP-EOM-EA-CCSD results obtained for the resonance position and width are compiled in Table VII together with several values available from the literature. The optimal CAP strengths found in our calculations are 0.0015 (0.0037) a.u. for the zeroth-order result and 0.0119 (0.0025) a.u. and 0.0148 (0.0071) a.u. for the real and imaginary part of the first-order result obtained with the aug-cc-pVTZ+3s3p3d(C) (aug-cc-pVQZ+3s3p3d(C)) basis set. We refrained from including experimental results in Table VII except for the fixed-nuclei estimate by Berman et al.⁸⁹ (E_R = 2.32 eV, Γ = 0.41 eV), which is usually considered as the reference value in previous theoretical studies. We add that this value was derived through a fit to the experimental data using Feshbach's projection operator formalism.

Table VII illustrates that a confusing plethora of values for the resonance position and width of the ²Π_g state of

$N2−$ have been reported. One can see that our CAP-EOM-EA-CCSD results overestimate E_R relative to the value from Ref. 89, while Γ is underestimated. For our highest-level calculation (first-order CAP-EOM-EA-CCSD/aug-cc-pVQZ+3s3p3d(C)), the deviation is about 0.15 eV for E_R and 0.13 eV for Γ. For E_R such a deviation is generally considered acceptable for EOM-EA-CCSD when dealing with bound states. An assessment of the deviation in Γ is more difficult as a comparison to bound states cannot be made. The differences between our highest-level value and the rest of the results show however that the basis-set size and the correction for the CAP potential both make a sizable impact on E_R and Γ, but whereas these effects work in the same direction for the resonance position, the change in Γ is more involved. As for the first-order correction, the results obtained within the static-exchange approximation¹⁶ exhibit trends similar to those observed in the present study.

Table VII also shows that most approaches yield too high values for the resonance position, whereas the reference value for the resonance width from Ref. 89 was often surprisingly well reproduced. The impact of electron correlation is illustrated by the comparison of high-level correlated methods to lower levels of theory. The position of the resonance state is consistently calculated to be above 3 eV with HF, density functional theory (DFT), and CIS (configuration interaction singles) based methods, regardless of whether stabilization techniques, complex scaling, or CAPs are employed, while the use of correlated methods leads to a significantly better agreement with the reference value. The only notable exception is the complex-scaled MRCI (multireference configuration interaction) result (1.38 eV) from Ref. 10, which is almost 1 eV below the reference value. Interestingly, CAP-augmented MRCI calculations⁶³ with a rather similar basis set yielded a quite different resonance position (2.97 eV).

Compared to its effect on the resonance position, the role of electron correlation for the width Γ is less clear. For example, the complex-scaled HF and CAP-HF calculations from Refs. 81 and 27 agree with the reference value within 0.03 eV and 0.015 eV, respectively, while some correlated calculations led to deviations of more than 0.2 eV.⁶³ In fact, it was stated explicitly⁶² that several values reported in the literature might have benefitted from error cancellation. We note that most authors reported values for Γ that were higher than the reference value with some low-level approaches overestimating the width by a factor of more than two, whereas our calculations underestimated Γ.

Regarding the basis-set dependence, a comparison between different schemes is hampered by the fact that a variety of different basis sets has been used in previous studies. Bearing in mind our findings from Sec. IV D, it seems justified to conclude for CAP-based methods that basis-set effects may account at least partly for the differences between values for Γ reported by different authors. In addition, the results from Ref. 80 suggest that, as compared to our CAP-augmented EOM-EA-CCSD scheme, the stabilization method combined with EOM-EA-CCSD leads to a somewhat faster convergence with respect to basis-set size. However, the highest-level (aug-cc-pV5Z+3p) results obtained with the latter method (E_R = 2.49 eV, Γ = 0.496 eV) still show a sizable deviation from the reference values.

While we have employed the ²Π resonance of CO⁻ as a test system in Sec. IV, we will consider it in this section with a different focus, i.e., we will compare our results to those obtained using other methods. In contrast to

$N2−$⁠, which has been studied frequently, comparatively few results have been reported for the isoelectronic CO⁻.^{24,26,70,82,90–93} Table VIII compiles the values available from the literature together with some representative values from Sec. IV.

Regarding the resonance position, Table VIII shows that most theoretical values are higher than the experimental value (1.50 eV)^78,94 as in the case of

$N2−$⁠. Also, similar to

${\rm N}_2^-$

$N2−$⁠, the resonance position is clearly overestimated in the static exchange approximation, in which electron correlation is neglected. Furthermore, CAP-EOM-EA-CCSD results obtained using different bases vary by up to 0.75 eV, which demonstrates the sizable impact of the basis set. We note, however, that our results obtained using the aug-cc-pVXZ+3s3p3d(C) bases approach the experimental value with growing basis-set size and that the first-order correction improves the resonance position with respect to experiment. Our highest-level result (1.762 eV, first order, aug-cc-pV5Z+3s3p3d(C) basis set) deviates from the experimental value by less than 0.3 eV.

The available values for the resonance width of CO⁻ differ by more than an order of magnitude as again illustrated by Table VIII. The largest value reported (1.65 eV) was obtained in the static exchange approximation, while the narrowest width (0.08 eV) was computed from the Σ² decouplings of the electron propagator, a pattern that is again similar to

$N2−$⁠. The CAP-EOM-EA-CCSD results for the resonance width from the present work and Refs. 24 and 26 vary by up to 0.68 eV, which illustrates once more the great influence of the basis set. The first-order correction improves the resonance width considerably with our highest-level result (0.604 eV, first-order, aug-cc-pV5Z+3s3p3d(C)) differing by 0.204 eV from the experimental value. However, we finally note that the experimental values for the resonance position and width of CO⁻ from Ref. 94 are not strictly comparable to the fixed-nuclei extrapolation for

${\rm N}_2^-$

$N2−$ from Ref. 89, which further complicates a rigorous assessment of the accuracy of theoretical approaches.

$C2H2−$ is a relatively well studied system and a number of theoretical^{26,80,82,95–97} as well as experimental^98–103 values for the resonance position and width are available from the literature. Those values as well as the results we obtained using our new CAP-EOM-EA-CCSD approach are compiled in Table IX. The optimal CAP strengths corresponding to our results are 0.0036 a.u. for the zeroth-order value and 0.0071 a.u. and 0.0058 a.u. for the real and imaginary part of the first-order value.

Table IX shows that our results for the resonance position are in qualitative agreement with those obtained from most experiments as well as from other theoretical approaches. Only when using the trapped electron method^98,99,103 considerably lower (∼0.7 eV) resonance positions were found. Our zeroth-order CAP-EOM-EA-CCSD result (2.655 eV) agrees within 0.05 eV with the experimental values from Refs. 99 and 100, and 102. We note that the first-order correction lowers the CAP-EOM-EA-CCSD result for E_R by about 0.2 eV bringing it closer to the experimental value from Ref. 101, but basis-set effects may have an impact of similar magnitude.

Theoretical values for the resonance width vary between 0.19 eV and 1.11 eV and only two rough estimates of 0.8 eV⁹⁹ and >1.0 eV¹⁰⁰ are available from experiment. Our calculations qualitatively confirm these two estimates with the zeroth-order result (0.979 eV) being closer to one value and the first-order result (0.831 eV) being closer to the other value. As mentioned for the resonance position, basis-set effects may have a sizable impact so that an ultimate decision between the two experimental values cannot be made.

Similar to CO⁻, the ²B_2g resonance in

$C2H4−$ has been chosen as a benchmark system in Sec. IV and here we will compare it with previously reported values. The ²B_2g resonance in

${\rm C}_2{\rm H}_4^-$

$C2H4−$ has been studied quite extensively by both experimental^104–106 and theoretical methods.^{66,67,80,107–110} Experimental measurements by electron scattering and electron impact techniques^104–106 located the position of the ²B_2g resonance around 1.8 eV with a width of Γ = 0.7 eV. Theoretically, this resonance has been studied by a wide variety of methods including complex scaling,^107,109 CAP-based approaches,^66,67 stabilization,⁸⁰ as well as other techni-ques.^108,110 The reported theoretical values vary from 1.77 to 2.62 eV for the position of the resonance and from 0.11 to 1.32 eV for the width, i.e., by more than an order of the magnitude in the latter case. All values along with the results from our method are summarized in Table X.