Calculation of the Lowest Resonant States of H − and Li by the Complex Absorbing Potential Method

: The analysis of the features of the method of complex absorbing potential (CAP) is carried out for a single-channel problem with an explicit parameterization of the scattering matrix. It is shown that there can be several types of CAP trajectories depending on the choice of the initial conditions. In any case, the estimation of the resonance parameters from the position of the optimal trajectory point can lead to a systematic error or an ambiguous result. In special cases, the search for the optimal point can be replaced by the averaging over a closed section of the trajectory. The CAP trajectories constructed in the H − and Li resonance calculations correlate well with the model trajectories, which have a curl around the resonance. The averaging over a closed area of the trajectory leads to better estimates of the energy and width of the resonance in comparison with the technique of searching for the optimal point.


Introduction
The possibility of calculating the energy and width of the resonance state by introducing an imaginary potential into the Hamiltonian was first shown more than thirty years ago [1][2][3][4][5]. The procedure proposes constructing the dependence of the complex eigenvalue of the stationary Schrödinger Equation by varying the scaling factor with the subsequent estimation of the resonance position as the geometric center of the closed interval of the obtained dependence (trajectory). This technique was substantiated because using a finite basis set leads to the appearance of a trajectory section closed around the resonance, which should converge to the exact position of the resonance in a complete basis set [5].
Later, an attempt was made to explain the introduction of the absorbing potential into the Hamiltonian by compensating for the incompleteness of the basis in the complex coordinate method (CCM). As in the CCM, in this case, the resonance position was estimated from the "optimal" point of the trajectory, for which the derivative of the energy concerning the scaling factor is equal to zero [6].
The key point in the application of the absorbing potential became the work [7,8], in which a theoretical foundation was given, and the results of calculations of the metastable 2 Π g state of the anion of the N 2 molecule were presented. Demonstration of the successful application of the method called the "complex absorbing potential method" (CAP) and the relative simplicity of its implementation in comparison with the CCM [9] were the reasons for its application to the calculation of resonances in many electron systems .
The CAP method's development was aimed at multiconfiguration quantum chemical methods that provide the correct account of electron correlation effects. At present,

Theory
Let there be a radial equation Atoms 2021, 9, 72 3 of 17 where is the cyclic Planck constant, m is the particle mass, r ≥ 0 is an independent variable, V(r) is the potential, E is the eigenvalue, and ψ(E, r) is the solution satisfying the initial condition ψ(E, r = 0) = 0. Here we will assume that for the function ψ(E, r) and its derivative, the continuity conditions are satisfied at the point r = R, i.e., ∂ψ int (E, r) ∂r The int and ext indexes here are introduced for the inner (0 ≤ r ≤ R) and outer (r ≥ R) domains.
Let us assume that there are functions {φ n (r)} in the outer domain and a system φ n (r) biorthogonally conjugated to it. Then the elements of these systems satisfy the condition where δ nm is the Kronecker symbol and the function ψ ext (E, r) can be represented as series with the coefficients A n = ∞ R ψ ext (E, r) φ n dr and A n = ∞ R ψ ext (E, r)φ n dr.
Multiplying Equation (1) by the functions φ n (r) (n = 1, 2, . . . ) on the left and integrating over the range r ∈ [R, ∞), we obtain the equation where W( f , g) = f ∂g ∂r − g ∂ f ∂r is the Wronskian, resulting from integration by parts of the term containing the second derivative. This term can also be represented in the integral form provided if the Bloch operator [59,60] is introduced, where δ(r − R) is the Dirac delta function, and b = (∂ln ψ ext (E, r)/∂r)| r=R . Expand the function ψ ext (E, r) by functions φ n (r), and substitute the result into Equation (6). Then let us combine the coefficients A n , functions φ n (r), and their derivatives in the obtained equations into the row vectors A = A 1 . . . A n , F = φ 1 . . . φ n and F / = ∂ φ 1 /∂r . . . ∂ φ n /∂r whereas the integrals into the matrix H ext with elements which after reducing similar terms transforms to the form where R(E)| r=R = b −1 is the element of R-matrix. Thus, solving Equation (11) with the given boundary condition b −1 = R(E) r=R , it is possible to find a value E that satisfies the radial Equation (1). Since the eigenfunction in Equation (1) is assumed to be continuous at the point r = R, the solution of Equation (11) for the given boundary condition can be considered as one of the options for combining the R-matrix with the approximate solution (1) in the basis of functions that are squareintegrated over the range r ≥ R. The possibility of using the simplified form (11) to combine R-matrix and the complex scaling method was discussed earlier, but the practical side of the issue was reduced to solving the generalized eigenvalue problem [61,62]. Case (11) was also used to solve the model problem by CAP [7].
Let us consider the application of Equation (11) to search the resonant pole by the CAP [7]. We will assume that in the region r ≥ R, there is an absorbing potential with a complex scaling factor s (Im(s) < 0) [7,25,55]. To expand ψ ext (E, r), let us choose the Hermite functions φ n (x) = C n H n (x) exp −x 2 /2 , where H n (x) is the n-th degree Hermite polynomial x = β(r − R), β is a parameter and is a normalization factor. Functions (13) are eigenfunctions for the har- and eigenvalues E n = 2 β m (n + 1/2) = λ m (n + 1/2), where n = 0, 1, 2, . . .. In an interval x ∈ (−∞, +∞), functions φ n (x) form a complete and closed orthonormal system for β > 0 [66,67]. Using the orthonormality of the functions φ n (x) and representing the integral over the domain x ∈ (−∞, +∞) as the sum of the integrals it can be noted that for an even integrand at x ∈ [0, +∞) the condition is satisfied, i.e., orthogonality is preserved.  Thus, it follows from expressions (15) and (16) that for real r and β > 0 Hermite functions of even and odd degrees form two mutually nonorthogonal systems of orthonormal functions in the range x ∈ [0, +∞).
When passing to complex values of the parameter β properties (15) and (16) are preserved, which allows defining the biorthogonal conjugate functions by the condition φ n (x) = φ n (x) without restrictions on the values β. It should be noted that the completeness of systems of Hermite functions of even (or odd) degrees on the entire complex plane of x values has not been proven, but it has been shown that these systems are closed but not complete [69,70].
Hence it follows that for the expansion of ψ ext (E, r) that only Hermite functions of even degrees can be used. Then ∂φ n (x)/∂r| r=R = 0, and consequently F / r=R = F / r=R = 0, which leads to the equation particular cases of which were introduced in [7,60,62].
Building the trajectory E(s) using the CAP can be carried out in two ways: fix the value of the parameter λ and solve Equation (17) for all required values of the parameter s, or postulate s = λ at each point of the trajectory. The first method is used in most studies of many-electron systems since it does not require multiple calculations of one-and two-electron integrals; the second one was introduced if considering the model problem in work [7].
Substituting the Hermite functions into the expressions for the elements of the matrix H ext , we obtain where V nm = ∞ K φ n (r)(r − R) 2 φ m (r)dr. Using the recursion relation, H n+1 (x) − 2xH n (x) + 2nH n−1 (x) = 0, it can be shown that H ext is a tridiagonal matrix with nonzero entries V nn = (2n + 1/2)/β, V nn+1 = [(n + 1/2)(n + 1)] 1/2 /β and V nn−1 = [n(n − 1/2)] 1/2 /β. Suppose the value of the parameter λ is fixed. In that case, after the substitution of matrix elements (18) and the corresponding values of the Hermite functions into it, the expression (17) takes on its final form. If s = λ the matrix H ext becomes diagonal (H ext nm = E m δ nm ), which allows replacing the right-hand side of (17) with the sum which includes terms with n = 0, 2, 4, . . .. Substituting the values of E n and φ n (r) in (19), taking into account the change in the normalization ( C n → C n √ 2 ), and performing the necessary transformations, we obtain where Γ(x) is a gamma function, z = E 2 m λ = k 2 4β . k = √ 2mE/ 2 is a wave number, and l = 0, 1, 2, . . . . If the summation in Formula (20) is carried out over an infinite number of terms, then their sum can be represented in a more convenient form where z = 5/4 + n (n = 0, 1, 2, . . . ). Details of transformation (21) are given in [68,71] (see the Appendix A). The substitution of (21) into (20) followed by multiplication by −ik gives the final expression

Calculation Method
In the numerical solution of Equations (17), (20) and (22), the function ψ int (E, r) at r → R was taken in the form where χ (±) = e ±iθ , θ = kr − z 0 k ln(2kr) + argΓ 1 + i z 0 k , z 0 is the target charge, k is the wavenumber, and S(k) is the scattering matrix. For a system with isolated Breit-Wigner resonances, the scattering matrix can be represented as the sum where S 0 (k) = e i2δ 0 is the nonresonant component of the scattering matrix, δ 0 is the background part of the phase shift, E n = E 0n − iΓ n /2 is the position of the n-th pole, and T n are some complex numbers [72]. Using the unitarity of the scattering matrix at real values of the energy, (24) can be replaced by the expression Parametrization of the scattering matrix (25) assumed a single resonance with E 0 = 3.42639031 − 0.00638724i a.u., which corresponds to the energy of the lowest resonance in the potential [73], and the nonresonant component of the phase shift δ 0 = π/4. The particle mass was set equal to m = 1 a.u., R = 20 a.u. Equations (17), (20) and (22) were solved iteratively with an energy convergence threshold of 10 −10 a.u. The method [74] was used to calculate the gamma functions of a complex argument. The inversion of the tridiagonal matrix in (17) was carried out by the method proposed in [75].
In the calculations of the Li atom, we used a technique previously tested on anions H − and He − [46,83]. Accordingly, the wave function was represented as a linear combination of three-electron configuration state functions (CSF), including the correlation and asymptotic components. Both components included complete sets of CSFs built from orthonormal AOs of the form where λ is the scaling factor (λ > 0) and N nl is the normalizing factor [83][84][85][86]. The AO basis for the correlation component included the functions χ nlm and |m| ≤ l, whereas the asymptotic one χ nlm → r with l ≤ 1, 5 − l ≤ n ≤ 26 − l, and |m| ≤ l (after this referred to as the basis of AO 5s4p3d2f/22s22p). The calculations used the absorbing potential (12) with an imaginary scaling factor (−8.9 ≤ ln(−Im(s)) ≤ 8.0).
To select the point of superposition of the absorbing potential calibration, calculations were performed with λ = 3.5, corresponding to the minimum energy of the first excited state of the target (the closed channel with the lowest energy). Analysis of CAP trajectories obtained for 10 ≤ R ≤ 17 a.u. revealed that for the chosen basis, AO closed sections appear on trajectories at 14 ≤ R ≤ 16 a.u. Subsequent calculations were carried out for R = 15 a.u. at 3.1 ≤ λ ≤ 3.7. The resonance position was calculated using two methods: averaging over a closed section of the trajectory (A) and searching for the optimal point (minimum of the function |s∂E(s)/∂s|) (B) [7,8]. The minimum of the function determined the position of the optimal point where η = Im(s). For an independent assessment of energy and width of 2 P 0 resonances, we used the method of stabilization by a modified Coulomb potential [83][84][85][86][87]. CAPtrajectories for the lowest 1 S resonance of H − were calculated using an 9s8p7d/26s26p26d AO basis set. The details of calculations can be found elsewhere [46].

Results and Discussion
Regarding Hermite functions of even degrees, the expansion of Equation (17) solutions leads to three model problems. The first of them ψ ext (E, r) is a solution to (22) and is expanded in a series in an infinite number of scalable (λ = s) functions (I). The second is reduced to solving Equation (20) and uses an expansion of ψ ext (E, r) in a series in a finite number of scalable functions (II). The latter corresponds to the solution (17)    Since the E(s) is a multivalued function at s > 0, it follows that for the absorbing potential with a complex scaling factor s = ζ − iη (ζ, η > 0), with the appropriate values of ζ = Re(s), it is possible to construct several trajectories converging to the same real value E(s) at η → 0 . It was found that in the case of problem I for an uncharged target (z 0 = 0), in  (Figure 2). Trajectories of the first type (ζ < 0.1) have a quasilinear section near the resonance pole, whereas the second type (0.3 < ζ < 1) bypasses the resonance, forming a section that is closed around it. The third type of trajectory (0.1 < ζ < 0.3) has several closed sections near the pole. Note that the reasons for the appearance of quasilinear trajectory segments in the vicinity of the resonance were discussed earlier, and the corresponding explanations remain valid for the model problem I [55,58]. Trajectories of all types have optimal points (27), which can be considered as estimates of the resonance position. However, in the case with 0.3   the resonance position cannot be estimated unambiguously since the trajectories have more than one optimal point. At 0.3   the resonance position can be estimated uniquely, but the optimal point is away from the pole, leading to a significant systematic error (Table 1). When passing to expansions of finite length (problem II), the shape of the trajectories changes, but the number of optimal points remains the same (Figures 3-5). Simultaneously, the accuracy of the estimation of the resonance position is significantly reduced. (Table 2).  Trajectories of all types have optimal points (27), which can be considered as estimates of the resonance position. However, in the case with ζ < 0.3 the resonance position cannot be estimated unambiguously since the trajectories have more than one optimal point. At ζ > 0.3 the resonance position can be estimated uniquely, but the optimal point is away from the pole, leading to a significant systematic error (Table 1). When passing to expansions of finite length (problem II), the shape of the trajectories changes, but the number of optimal points remains the same (Figures 3-5). Simultaneously, the accuracy of the estimation of the resonance position is significantly reduced. (Table 2).         For problem III it was found that in the case Re (s) = 0, there are only trajectories of the first type (Figures 3-5), while for Re (s) > 0, the type of trajectory depends on the chosen value Re (s) like for the problem I. The last conclusion for problem II is also valid for problem III (Table 3). Since the solution of problems I-III for a charged target (z 0 = 1) gives practically the same results as the above case, we give only a representative part of the results (Tables 1-3). It should be added that the systematic shift of the resonance energy by~10 −5 a.u. for a charged target (Table 1) is the result of the finite-difference calculation of the minimum of function (27). Thus, two methods of calculating the resonance parameters can be used for the considered model problems. First, one can use the absorbing potential with an imaginary scaling factor and choose one of the optimal points, which will require justifying the possibility of introducing a selection criterion. Second, the form of trajectories can be simplified by using a complex scaling factor for absorbing potential, like in the technique [28]. In this case, the resonance position can be estimated by averaging over a fragment of the trajectory [5]. However, both results will contain a systematic error, the value of which depends on the length of the solution expansion. In the calculations of H − and Li resonances, it was found that CAP trajectories have a curl that is closed around the reference resonance position (Figures 6 and 7). For the H − anion, the energy and width of the lowest 1 S 0 resonance, calculated by averaging over the closed section of the trajectory (A), are in better agreement with the reference values than the energy and width of the resonance estimated from the position of the optimal point (27). In addition, there is a dependence of the energy and resonance width on the choice of the value of the AO scaling factor (26) ( Table 4), which is consistent with the results of calculations by the CAP and stabilization methods [46,84,86].     [51,53]. Sign + in panel (b) indicates the position of the calculated resonance [54]. Data are presented in a.u.
For 2 P 0 resonances of Li, methods A and B provide energy estimates that are close in accuracy. Both methods predict the value of the width of the lowest resonance, which is in poor agreement with the experimental value, but estimates of the width by method (A) are consistent with the independent theoretical results ( Table 5). Estimates of the second resonance width by method A agree better with theoretical values than estimates by method B. Calculations of 2 P 0 resonances by the CAP and stabilization methods give close values of energies and widths for both resonances ( Table 6). Comparison of the results of Hand Li calculations by the CAP and stabilization methods allows us to conclude that the accuracy of the ab initio calculation is primarily determined by choice of the AO basis and the method for the approximate solution of the stationary equation.  [51,53]. Sign + in panel (b) indicates the position of the calculated resonance [54]. Data are presented in a.u. For 2 P 0 resonances of Li, methods A and B provide energy estimates that are close in accuracy. Both methods predict the value of the width of the lowest resonance, which is in poor agreement with the experimental value, but estimates of the width by method (A) are consistent with the independent theoretical results ( Table 5). Estimates of the second resonance width by method A agree better with theoretical values than estimates by method B. Calculations of 2 P 0 resonances by the CAP and stabilization methods give close values of energies and widths for both resonances ( Table 6). Comparison of the results of H − and Li calculations by the CAP and stabilization methods allows us to conclude that the accuracy of the ab initio calculation is primarily determined by choice of the AO basis and the method for the approximate solution of the stationary equation.  [51][52][53]. c theory, Feshbach projection formalism [54]. d theory, saddle-point technique [79]. Table 6. Energy (a.u.) and width (×10 −5 a.u.) of the first and the second 2 P 0 resonant states of Li atom calculated by stabilization method (1) and CAP (2). For CAP, the results were obtained by method A. λ 2 P 0 [1s(2s2p) 3

Conclusions
In this paper, we analyzed the features of the CAP method for a single-channel problem with an explicit parameterization of the scattering matrix. It was found that there are three types of trajectories, the of which is determined by the choice of a branch of the stabilization graph. The first type of trajectory is characterized by a quasilinear segment close to the resonance pole, and a closed segment around the resonance characterizes the second type of trajectory. There are also trajectories of an intermediate type, with several closed sections near the resonant pole. When passing to a model with a finite number of terms in the expansion of the solution of the stationary equation, the number of optimal points is preserved for the trajectories, but it becomes difficult to classify the type of trajectory.
On the contrary, CAP trajectories for H − and Li correlate well with model trajectories of the second type. Accordingly, the estimates of the resonance parameters by averaging over the closed section of the trajectory are in better agreement with the results of independent calculations.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Based on the results from the book [68], it is possible to propose several methods for transforming the series where a and b are complex numbers. First, we can use the table integrals (7.511) and (7.512.10) and replace the series (A1) with the Gaussian hypergeometric function, and then replace the latter with the ratio of the gamma functions [71]. Secondly, the same result can be obtained using formulas for the sums of the series obtained using the Mellin transformation ( §1.4 and §2.4 [68]). However, implicit restrictions on the values of parameters a and b are assumed in both cases, which requires additional study. Alternatively, it is convenient to use the method from §2.1.3 of the book [68]. Replacing the expression b + n = Γ(b + n + 1)/Γ(b + n) and entering the Pochhammer symbols (a) n = Γ(a + n)/Γ(a), series (A1) can be represented as .