1. 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 [
10,
11,
12,
13,
14,
15,
16,
17,
18,
19,
20,
21,
22,
23,
24,
25,
26,
27,
28,
29,
30,
31,
32,
33,
34,
35,
36,
37,
38,
39,
40,
41,
42,
43,
44,
45,
46,
47].
The CAP method’s development was aimed at multiconfiguration quantum chemical methods that provide the correct account of electron correlation effects. At present, computational procedures have been tested using the configuration interaction (CI) formalism [
7,
8,
11,
15,
17,
18,
21,
22], various variants of cluster expansion (EOM-CC, SAC-CI, FSMRCC) [
26,
27,
29,
32,
33,
34,
35,
36,
38,
39,
40], algebraic diagrammatic approach [
24,
47] and quasi-degenerate perturbation theory (XMCQDPT2) [
41]. The accuracy of the proposed approaches can be compared only for the lowest resonances
(
2Π
g) and
(
2Π), although simple systems (
,
, Li, etc.) could be used as test systems, for which the positions of the resonances are known with sufficient accuracy [
48,
49,
50,
51,
52,
53,
54].
For the lowest resonance
(
2Π) at an equilibrium internuclear distance, the energy estimates are in the range 1.90–2.08 eV [
31,
32,
33,
43,
47,
55] and exceed the experimental value of 1.50 eV [
56]. The resonance widths are estimated in the range 0.38–0.63 eV [
31,
32,
33,
43,
47,
55] at the experimental value of 0.4 eV. The maximum theoretical value of the resonance width (1.3 eV) was obtained in [
43]. The sensitivity can explain the error of theoretical results to a quality of the basis of atomic orbitals (AO) and parameters of the absorbing potential [
32,
33]. For the lowest resonance
(
2Π
g), the smallest deviations from the reference values of energy (2.316 eV) and width (0.414 eV) [
57] are less than 0.1 eV [
27,
49] and 0.12–0.2 eV [
24,
26,
27,
32,
43,
47,
55], respectively. In other calculations, deviations in the energy and width of the resonance have larger values [
7,
11,
18,
24,
26,
30,
32,
33,
41,
43,
55]. The calculation accuracy depends on the same factors as in
, i.e., the quality of the AO basis and the parameters of the absorbing potential [
18,
30,
32].
Simultaneously with the search for methods for solving the electronic problem, attempts to improve the accuracy of calculations by changing the methods of building trajectories or estimating the position of the resonance were made. It has been shown that replacing imaginary absorbing potential with a complex one makes it possible to eliminate nonphysical optimal points of the CAP trajectory [
28]. Moreover, a new technique for estimating the position of the resonance was proposed, based on the construction of Padé approximants [
55,
58]. The application of this technique for
(
2Π
g) and
(
2Π) resonances has improved the accuracy of the calculation [
55].
A feature of the implementation of the CAP is the use of a basis of atomic orbitals (AO), square-integrable over the entire three-dimensional space. This method of constructing the AO basis preserves the continuity of the wave function at the boundary of the overlapping region of the absorbing potential but does not allow to separate the errors arising from the electron correlation, as well as the errors of the CAP parameterization. This feature can be removed by stitching the solutions at the boundary of the overlapping region of the absorbing potential using the R-matrix formalism [
59,
60]. Similar approaches were previously proposed to calculate the parameters of scattering resonances by the complex coordinate method [
61,
62] and solving problems of reaction scattering [
25,
63,
64,
65]. The combination of the R-matrix with the CAP method was also used in the construction of the theory of the CAP [
7].
The approach proposed in [
7] makes it possible to explicitly parametrize the R-matrix and, therefore, to obtain the CAP-trajectories free from errors of the approximate solution of the Schrödinger Equation in the inner region. By comparing the behavior of the trajectories obtained for an explicitly parameterized R-matrix and real systems, it is possible to formulate the requirements necessary for the successful application of the CAP [
46].
This work aims to assess the possibilities of the combined approach [
7]. For comparison, we chose the lowest
1S resonance of the H
− anion and
2P
0 resonances of the lithium atom. The theoretical formalism is briefly presented in the second part of the article. The calculation details are given in the third part, followed by a discussion of the main results and a conclusion.
2. Theory
Let there be a radial equation
where
is the cyclic Planck constant,
m is the particle mass,
is an independent variable,
is the potential,
E is the eigenvalue, and
is the solution satisfying the initial condition
. Here we will assume that for the function
and its derivative, the continuity conditions are satisfied at the point
, i.e.,
The int and ext indexes here are introduced for the inner () and outer () domains.
Let us assume that there are functions
in the outer domain and a system
biorthogonally conjugated to it. Then the elements of these systems satisfy the condition
where
is the Kronecker symbol and the function
can be represented as series
with the coefficients
and
.
Multiplying Equation (1) by the functions
(
n = 1, 2, …) on the left and integrating over the range
, we obtain the equation
where
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
is the Dirac delta function, and
.
Expand the function
by functions
, and substitute the result into Equation (6). Then let us combine the coefficients
, functions
, and their derivatives in the obtained equations into the row vectors
,
and
whereas the integrals into the matrix
with elements
, which leads to the matrix equation
Multiplying (9) on the right by the column vector
, where
and
, we get the equation
which after reducing similar terms transforms to the form
where
is the element of
R-matrix.
Thus, solving Equation (11) with the given boundary condition
, 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
, 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 square-integrated over the range
. 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
, there is an absorbing potential
with a complex scaling factor
s (
) [
7,
25,
55]. To expand
, let us choose the Hermite functions
where
is the
n-th degree Hermite polynomial
,
is a parameter and
is a normalization factor. Functions (13) are eigenfunctions for the harmonic oscillator problem with mass
m, potential
,
(
), and eigenvalues
, where
.
In an interval
, functions
form a complete and closed orthonormal system for
[
66,
67]. Using the orthonormality of the functions
and representing the integral over the domain
as the sum of the integrals
it can be noted that for an even integrand at
the condition
is satisfied, i.e., orthogonality is preserved. For an odd integrand, expressing the Hermite polynomials in terms of the generalized Laguerre polynomials
and
, where
[
68], we have
Thus, it follows from expressions (15) and (16) that for real r and Hermite functions of even and odd degrees form two mutually nonorthogonal systems of orthonormal functions in the range .
When passing to complex values of the parameter
β properties (15) and (16) are preserved, which allows defining the biorthogonal conjugate functions by the condition
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].
Equation (11) assumes that
and
(or
and
). Hence it follows that for the expansion of
that only Hermite functions of even degrees can be used. Then
, and consequently
, which leads to the equation
particular cases of which were introduced in [
7,
60,
62].
Building the trajectory
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
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
, we obtain
where
Using the recursion relation,
it can be shown that
is a tridiagonal matrix with nonzero entries
,
and
. 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
the matrix
becomes diagonal (
), which allows replacing the right-hand side of (17) with the sum
which includes terms with
.
Substituting the values of
and
in (19), taking into account the change in the normalization (
), and performing the necessary transformations, we obtain
where
is a gamma function,
.
is a wave number, and
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
(
). Details of transformation (21) are given in [
68,
71] (see the
Appendix A). The substitution of (21) into (20) followed by multiplication by
gives the final expression
3. Calculation Method
In the numerical solution of Equations (17), (20) and (22), the function
at
was taken in the form
where
,
,
is the target charge,
k is the wavenumber, and
is the scattering matrix. For a system with isolated Breit–Wigner resonances, the scattering matrix can be represented as the sum
where
is the nonresonant component of the scattering matrix,
δ0 is the background part of the phase shift,
is the position of the
n-th pole, and
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
a.u., which corresponds to the energy of the lowest resonance in the potential [
73], and the nonresonant component of the phase shift
. The particle mass was set equal to
a.u.,
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].
The CAP method was used to calculate the two lowest
2P
0 [1s(2s2p)
3P] and
2P
0 [1s(2s2p)
1P] resonances of the lithium atom and
1S resonance of H
−. For the lowest
2P
0 [1s(2s2p)
3P] resonance, there are experimental estimates of energy and width [
51,
52,
53,
76], for the next
2P
0 [1s(2s2p)
1P] resonance where only the energy is known [
77,
78]. The experimental results are in satisfactory agreement with the theoretical ones [
54,
79,
80,
81,
82].
In the calculations of the Li atom, we used a technique previously tested on anions
and
[
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 (
) and
is the normalizing factor [
83,
84,
85,
86]. The AO basis for the correlation component included the functions
with
,
, and
, whereas the asymptotic one
with
,
, and
(after this referred to as the basis of AO 5s4p3d2f/22s22p). The calculations used the absorbing potential (12) with an imaginary scaling factor (
).
To select the point of superposition of the absorbing potential calibration, calculations were performed with
, 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
a.u. revealed that for the chosen basis, AO closed sections appear on trajectories at
a.u. Subsequent calculations were carried out for
R = 15 a.u. at
. 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
) (B) [
7,
8]. The minimum of the function determined the position of the optimal point
where
For an independent assessment of energy and width of
2P
0 resonances, we used the method of stabilization by a modified Coulomb potential [
83,
84,
85,
86,
87]. CAP-trajectories for the lowest
1S resonance of H
− were calculated using an 9s8p7d/26s26p26d AO basis set. The details of calculations can be found elsewhere [
46].
4. Results and Discussion
Regarding Hermite functions of even degrees, the expansion of Equation (17) solutions leads to three model problems. The first of them is a solution to (22) and is expanded in a series in an infinite number of scalable () functions (I). The second is reduced to solving Equation (20) and uses an expansion of in a series in a finite number of scalable functions (II). The latter corresponds to the solution (17) with an expansion in a finite number of functions with a constant scaling factor () (III).
For positive values of the parameter
s, the behavior of the branches of the function
depends on the definition of the parameter λ. In cases I and II, the branches
converge to the point
and finite values at
. In case III at
the branches
have nonzero limits and at
converge to the same finite values as in cases I–II. At the parameter value
, the branches
, which have been built for the problems I-III, practically coincide (
Figure 1).
Since the
is a multivalued function at
, it follows that for the absorbing potential with a complex scaling factor
(
), with the appropriate values of
, it is possible to construct several trajectories converging to the same real value
at
. It was found that in the case of problem I for an uncharged target (
), in the range of values
there are three types of trajectories passing in the vicinity of the resonance (
Figure 2). Trajectories of the first type (
) have a quasilinear section near the resonance pole, whereas the second type (
) bypasses the resonance, forming a section that is closed around it. The third type of trajectory (
) 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
the resonance position cannot be estimated unambiguously since the trajectories have more than one optimal point. At
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 (
Figure 3,
Figure 4 and
Figure 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
, there are only trajectories of the first type (
Figure 3,
Figure 4 and
Figure 5), while for
, the type of trajectory depends on the chosen value
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 (
) gives practically the same results as the above case, we give only a representative part of the results (
Table 1,
Table 2 and
Table 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 (
Figure 6 and
Figure 7). For the H
− anion, the energy and width of the lowest
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].
For
2P
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
2P
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.