Resonances are perhaps one of the most striking phenomena in chemistry [1, 2]. Molecules in metastable states (so called resonances) have enough energy to ionize or dissociate but do not do it right away. They have finite lifetimes and decay to the products which can be electrons, ions and radicals. The decay rates can vary from case to case and can be different by many order of magnitudes and there may be several open decay channels.

Let’s consider the following illustrative triatomic (ABC) molecular reaction that occurs on a ground electronic potential energy surface, A + B C → A B C # → A B + C Or → A C + B Or → A + B C ∗ . (1) Where [ABC] ^# represents an activated complex that has enough energy to dissociate to several different products and, BC* is the diatomic molecule BC in an excited ro-vibrational state. This reaction takes place for a specific collision energy (within a given uncertainty) at which the activated complex is created in a well defined metatsable state. This metastable state is known as a resonance state, as time passes it decays into the reaction products. The energy of the activated complex above the lowest energy threshold (i.e., above the minimal energy in which it can dissociates) is the resonance position, E _res . The decay rate or width of the activated complex, Γ _res , is inverse proportional to the lifetime of the activated complex, where 2πℏ is the proportional parameter. The decay rate, Γ _res , is a sum over the partial decay rates into all the possible products. That is, Γ res = Γ A B + C + Γ A C + B + Γ A + B C ∗ .

The difference between the energy of the activated complex in a resonance (metastable) state and the energies of the different bound-state energies of the reaction products provide the relative kinetic energies of the AB + C, AC + B and A + BC* products. The energy and width of the activated complex and the partial decay rates can be calculated from a single eigenfunction of the time independent nuclear Schrödinger equation when outgoing boundary conditions are imposed [1–3]. Notice that by imposing outgoing boundary conditions on the non-equilibrium reaction presented in Eq. 1 will turn the real physical molecular Hamiltonian into a non-Hermitian Hamiltonian, as will be explained in details below. The resonance via Padé (RVP) method, which is the focus of this review, enables the calculations of the energy and decay rate of such an activated complex. In principle, also the partial widths can be calculated by RVP, but the way to do it is out of the scope of this review.

A second illustration examines the autoionization process in the helium atom. Helium has an infinite number of discrete bound states. The bound states of helium are associated with the ground and singly excited electronic states. Contrary, the doubly-excited states must ionize after a finite period of time, i.e., these states are resonance states, He ∗ ∗ → He + + e − . (2)

Let us prove it. The starting point in our proof is to neglect the electronic repulsion and approximate the energy of a doubly-excited helium (He**) as, E n 1 > 1 , n 2 > 1 no - repulsion = − Z 2 1 n 1 2 + 1 n 2 2 13.6  eV  , where Z = 2. This value must be lower than the exact energy of a doubly-excited helium state since the electronic repulsion was neglected. Therefore, whenever − Z 2 1 n 1 2 + 1 n 2 2 13.6  eV  > − Z 2 1 n ion 2 13.6  eV  ionization takes place, where n ₁, n ₂ are the neutral He electronic levels and n _ion is the level of He⁺. This inequality can be reformulated as, n 1 2 n 2 2 n 1 2 + n 2 2 > n ion 2 .

The lowest doubly-excited states occur when n ₁ = n ₂ = 2, for which n 1 2 n 2 2 n 1 2 + n 2 2 = 2 , i.e., the energies of these He** states are higher than the ionization threshold (corresponding to n _ion = 1). As long as n 1 2 n 2 2 n 1 2 + n 2 2 ≤ 4 (corresponding to n _ion = 2), there is only one open channel–autoionization decay into the helium cation in its ground electronic state. In case n 1 2 n 2 2 n 1 2 + n 2 2 > 4 , there will be several open decay channels involving also excited helium-cation states. Note that the spontaneous-emission time of He** is longer by several order of magnitude than the resonance lifetime. The spontaneous-emission time can be estimated using the Einstein coefficient for spontaneous emission (A, which is proportional to the calculated transition dipoles and the energy differences between the relevant states) [4], where the emission time is T = A ⁻¹. The resonance lifetimes are inverse proportional to the calculated widths. Therefore, auto-ionization would dominate the dynamics of He**.

These kind of resonances are referred to Feshbach type resonances [2] since their energies (for n ₁ = n ₂ = 2) are below their own ionization threshold, i.e., He*⁺(2s or 2p), but above the He⁺(1s) threshold. Thus, Feshbach resonances are associated with two-electron processes. Therefore, dynamical electronic correlation, in which one go beyond the mean field (Hartree Fock) approximation, must be considered. Imposing outgoing boundary conditions on the electronic solutions of the time-independent Scrödinger equation makes the electronic Hamiltonian non-Hermitian, as in the above molecular dissociation example. Thus, the resonance energy positions and decay rates of He**, and any chemical systems in autoionization states, can be immediately obtained from the eigenvalue spectrum of a non-Hermitian Hamiltonian (see the next section for more details). Moreover, below we show that RVP can be used to study the autoionization process in Eq. 2 with great accuracy.

Another type of metastable electronic states are the shape-type resonances [2], which lay above their own ionization threshold, therefore they represent a one-electron transition and they can be obtained within the Hartree-Fock approximation. This is a one-electron process, in which an electron tunnels through a potential barrier that results from a mean repulsion potential that corresponds to all the other electrons. Of course, for accurate calculations of energies and widths (inverse lifetimes) one should go beyond the mean field approximation. The simplest examples for electronic shape-type resonances are the ground electronic states of the molecular Hydrogen, Nitrogen and CO anions [5, 6]. Notice that in such anionic-resonance cases the decay process is referred to as autodetachment.

Uracil anion is an example for a biochemical system in which ionization and dissociation may occur simultaneously. Attachment of an electron to uracil leads to two types of electronic resonance anion states. Shape resonances appear as an electron is attached to one of the unoccupied π* orbitals of the neutral ground state of uracil. Alternatively, an electron can be attached to excited states of the neutral uracil, forming an electronica Feshbach resonances. In addition to the autoionization process associated electronic resonance states, also dissociative electron attachment (DEA) processes may follow upon the creation of a uracil anion. Uracil may undergo DEA into C₃H₃N₂O⁻ by eliminating CO and H [7]. These autoionization and DEA process may result in damage to the RNA strand. Below we show that RVP can be used in studying such processes, by calculating the complex energies of the lowest three shape-type states of the uracil anion as well as the transition dipoles between them.

All in all, electronic resonances refer to autoionization (e.g., Eq. 2) and nuclear resonances refer to a situation in which a molecule pre-dissociates (e.g., Eq. 1). Notice that autoionization can become a more complicated and “rich” phenomenon, as in the case of: Auger [8–11], ICD (interatomic Coulombic decay) [12–20], ETMD (electron-transfer mediated decay) [15, 20–23], etc. Nuclear shape-type and Feshbach-type resonances are described in detail in Chapter 2 in Ref. 1 that is dedicated to Resonances Phenomena in Nature. Of course, often the decay of the electron and nuclear resonance states may happen simultaneously. As for example, A + B C → A B C # → A B + C + + e − .

The RVP method, which is described below, is applicable to all these cases.

Let us first briefly explain the motivation to look for such a comparison. In scattering theory resonances are associated with wavefunctions or eigenstates of the time-independent Schrödinger equation (TISE), which only consist of outgoing waves at the asymptote. The physical reason is clear. The system is prepared in a metastable state and as time passes it breaks apart into sub-systems as described above. Solving the TISE with such outgoing boundary conditions (OBCs) results in a discrete spectrum with real and complex eigenvalues, which are associated with bound states and resonances, respectively [1–3]. Notice that such a spectrum characterizes the non-Hermitian quantum mechanics (NHQM) formalism, i.e., by imposing OBCs we turn the QM problem into the NH regime. The bound state eigenvalues are real and the corresponding eigenfunctions are exactly as usual (i.e, decay asymptotically to zero). The resonance eigenvalues are complex, E res − i 2 Γ res . The real part, E _res , corresponds to the resonance energy position, while the imaginary part, Γ _res , corresponds to the resonance decay rate (inversely proportional to the resonance lifetime) [1]. Alas, the asymptotes of the corresponding resonance eigenfunctions exponentially diverge. This asymptotic exponential divergent of the TISE is also obtained by wavepacket propagation calculations as illustrated in Figure 1A.

This seems to contradict the representation of resonances as the solutions of the time-dependent Schrödinger equation (TDSE) known as wavepackets. Since resonances are embedded in the continuous part of the spectrum they cannot be described using a single stationary eigenstate. Though, at any given time of the dynamics, the wavepackets are represented as superpositions of the (Hermitian) Hamiltonian eigenstates. Thus, wavepackets by definition are square integrable functions, i.e., their asymptotes decay exponentially and do not exponentially diverge.

The answer to this puzzle was given in Refs. [1, 3, 24]: before a wavepacket decays at the asymptote it actually diverges exponentially (at very large but finite distance from the interaction region). Figure 1B illustrates this point at different dynamic times. The conclusion is quite clear. It is appropriate to impose OBCs in order to describe a decaying system since the corresponding solution of the TISE with OBCs also display the same divergence as the wavepacket solution of the TDSE does.

So which method, out of the two, is recommended for calculating resonances? The answer is that neither of them can be recommended for realistic molecular systems. The reason is as follows: within Hermitian QM we cannot use the TISE and we are forced to solve the TDSE; this presents a major numerical difficulty when considering standard quantum chemistry packages for calculating resonance states. Within NHQM, solving the TISE with OBCs does not allow one to use square integrable basis sets, which transform the partial differential eigenvalue problem into a matrix eigenvalue problem, as in the Hermitian case. Therefore, a practical approach for describing resonances would be to transform the problem into the NH regime while retaining the square-integrability of the eigenstates. Such a procedure (the so-called complex scaling, which is described in the next section) can be used within many-body electronic structure formalism.

Let us briefly explain why the asymptote of the resonance wavefunction diverges exponentially in space but decays exponentially in time. The time-independent Schrödinger equation is given by, H ̂ | Ψ E 〉 = E | Ψ E 〉 , (3) where E = E _threshold + |ℏk|²/(2m) and E _threshold is the ionization/dissociation threshold energy. The asymptote of an eigenfunction associated with an above threshold energy and with only one open decay channel is given by, ⟨ r | Ψ E ⟩ → r → ∞ A out E e ikr r + A i n E e − i k r r . (4)

The scattering matrix is the ratio between the amplitudes of the incoming [A _in (E)] and the outgoing [A _out (E)] waves. The poles of the scattering matrix are complex and the associated wavevectors take discrete values, k n = k n R e − i k n I m = | k n | e − i ϕ n (5) when ϕ _n ≥ 0, for which A _in (E _n ) = 0 (see chapter 4 in Ref. [1]).

Therefore, since e i k n r = e i | k n | r ⁡ exp − i ϕ n → r → ∞ ∞ , (6) the resonance function is not part of the Hilbert space. That is, it is not a square integrable function since it exponentially diverges in space. Due to this spatial behavior the resonance wavefuctions can not be expanded with a basis set of square integrable functions, as in the calculations of bound electronic states. Therefore we need to find out how we can transform the electronic coordinates such that the resonance wavefunction will remain square integrable and can be described as a finite linear combination of square integrable basis functions. As for example Gaussians which are used in standard electronic structure calculations.

Now, let’s turn to the exponential decay of the resonance function in time. The time phase factor e − i E n t / ℏ = e − i Re E n t / ℏ e − I m E n t / ℏ → t → ∞ 0 , (7) where E _n = ReE _n + iImE _n (and for resonances ImE _n < 0) and the number of particles is conserved only when ℏ k n R e t = r (for explanation see Ref. [1] on the coupling of space and time even in the non-relativistic quantum mechanics framework). An interesting situation, which demonstrates the decay of the resonance function in time is presented in Figure 2. When the initial wavepacket mainly populates the two narrowest resonances of a model Hamiltonian, which is given in the caption of Figure 2. As one can see from Figure 2 first the shorter resonance decays and only then the resonance with longer lifetime (smaller width) decays. The plot of the log of the survival probability as function of time gives two straight lines that their slopes provide the decay rates of the two resonances initially populated.

It is straightforward to realize that Eq. 6, when r → re ^iθ and using Eq. 5, becomes, e + i k n r ⁡ exp i θ = e + i | k n | r ⁡ exp i θ − ϕ n → r → ∞ 0 , (8) when θ ≥ ϕ _n . And since (Eq. 5) k n = | k n | e − i ϕ n = 2 m E n − E threshold / ℏ , the condition θ ≥ ϕ _n yields, tan 2 θ ≥ 2 Γ n R e E n − E threshold . (9) Where {E _n } are the complex resonances energies, Γ _n = −2Im[E _n ], m is the mass of the emitted particle (e.g., an electron) and E _threshold is the ionization/dissociation threshold energy. For additional details see Chapter 5.2 in Ref. [1]. In such a case, the resonance state can be represented as a square integrable function, thus, it can be expanded in terms of localized basis functions (such as Gaussians), similarly to a bound state. Upon such complex coordinate rotation, the Hamiltonian becomes non-Hermitian and the complex resonance energies are obtained regardless of the value of the rotation angle θ inside the interval [θ _c , π/4]. (The upper limit on the θ interval is also explained in Chapter 5.2 in Ref. [1].)

This raises the question: is it possible to rotate the coordinate into the complex plane? If the physical potential is an analytical function, as in the case of atomic potentials, the coordinates of the system can be rotated into the complex half lower plane when θ _c < θ < π/4. However, the molecular electronic Hamiltonian within the Born-Oppenheimer approximation is singular at the nuclei positions and therefore a uniform analytical continuation into the complex plane by r _e → r _e e ^iθ is not allowed. This complication results in different methodologies for tackling this problem. The rigorous solution to this problem is to chose complex electronic coordinates that remain real inside an interaction volume, which include all the molecular nuclei. By that we avoid the singularities in the electron-nucleus attractive potential terms, and rotate the coordinate into the complex plane only out of the interaction volume by using θ. Section 4 presents several methods, which introduce such complex electronic coordinates that avoid these singularities. By using one of these complex electronic coordinates, the non-Hermitian (NH) Hamiltonian, H ̂ ( r e θ ) , is obtained. This NH operator can be represented using a finite matrix that is spanned by finite number of Gaussian basis functions. The non-Hermitian Hamiltonian matrix elements are given by H ζ , ζ ′ θ = ⟨ G ζ | H ̂ r e θ | G ζ ′ ⟩ , (10) where { G ζ ( r e ) } ζ = 1,2 , … N represents a set of Gaussian basis functions.

Here we are coming to a delicate point which is important for the RVP method. Rather than computing the complex matrix elements for a non-Hermitian Hamiltonian operator [ H ̂ ( r e θ ) ] , we can calculate the complex matrix elements by keeping the molecular Hamiltonian operator real [ H ̂ ( r e ) ] , as usual, and use complex transformed Gaussian basis functions instead. Moreover, since we want to move into the complex plane outside the interaction volume, we can complex transform only the diffuse Gaussians of the employed basis set [25]. To a good approximation, these diffuse Gaussians span the space outside the interaction volume. Therefore, transforming only these diffuse Gaussians will transform only the electronic coordinates outside the interaction volume. In this case, the non-Hermitian Hamiltonian matrix elements become H ζ , ζ ′ θ = ⟨ G ζ − θ ∗ | H ̂ r e | G ζ ′ − θ ⟩ , where the complex diffuse Gaussians are given by [26], G ζ r e − θ = x e − x N n y e − y N m z e − z N l ⁡ exp − ζ η − 2 r e − r N 2 ζ = 1,2 , … N , (11) and where η = αe ^iθ and the electronic vector position is r _e = (x _e , y _e , z _e ) centered on the nuclei R _N = (x _N , y _N , z _N ). Notice that one needs to avoid the complex conjugate in the matrix elements (the so called c-product, see Ref. [1]) and therefore we get, H ζ , ζ ′ θ = ⟨ G ζ + θ | H ̂ r e | G ζ ′ − θ ⟩ . (12)

Examining the NH Hamiltonian matrix elements, one can see that going into the complex plane can be made even simpler. From Ref. [27] we know that the Hamiltonian matrix elements, unlike the operator itself, can be analytically dilated. Therefore, one does not even have to use complex diffuse Gaussians to obtain these matrix elements. Instead, one can calculate these elements as a function of the parameter η but with θ = 0, i.e., calculate the matrix elements as a function of α. Then, analytically dilate this parameter into the complex plane by substituting η = αe ^iθ . In other words, to simplify things further and avoid the use of complex diffuse Gaussian basis functions, one can use real diffuse Gaussian basis functions, which depend on the real parameter, α, and then make this parameter complex by taking θ ≠ 0. Doing so, one can obtain the NH matrix elements by using real functions and the real Hamiltonian.

Needless to say that even this kind of simplification is not straightforward to apply and requires the modifications of the standard (Hermitian) quantum chemistry packages, which include variety of ab-initio methods (based on MP (Møller Plesset perturbation theory), CI (configuration interaction), CC (coupled cluster) and more.

Here we are coming to the new approach we developed, in which we take one step further in the analytic continuation direction and replace the analytic dilation of the Hamiltonian matrix elements with that of a single eigenvalue. A proof of this point is given in the next two sections. Obviously, there is a great numerical advantage to analytical continuation of a single eigenvalue over an N × N matrix elements, where N is the number of basis functions employed. Since the Padé approximant (see Section 2.2) is used for the analytical continuation the method is know as Resonance via Padé (RVP). There are additional numerical advantages to RVP. Standard NHQM methods commonly work directly in the complex plane in order to calculate the resonance eigenvalue [6, 25–40]. RVP belongs to a subgroup of methods, which move into the NHQM regime via analytical continuation [41–44] It is based on the stabilization technique [45–49], where the real energy levels are plotted as a function of a parameter (α as denoted above) that controls the diffuseness of the Gaussian basis functions (see details in Section 2.1). The stabilization calculation, which is computationally the most demanding step, is followed by a very “cheap” analytical continuation step [50] (see Section 3). Therefore, analytical-continuation methods that are based on the stabilization technique hold two distinct advantages: First, Matsika and co-workers showed that the computation time required by the stabilization technique is an order of magnitude lower than the time required by a NHQM approach that works literally in the complex plane [43]. Second, the versatility in generating the stabilization graph opens the door for various applications. That is, stabilization graphs can be calculated using any standard quantum chemistry packages, with any existing Hermitian electronic structure method, see for example Ref. [44].

In this section we explain the concept behind the RVP method and its connection to the analytical continuation of the real Hamiltonian matrix elements (HMEs); we stated above that in order to calculate resonances in the complex plane one can analytically dilate the real HMEs into the complex plane. This means that in order to calculate resonances one needs to obtain the matrix elements as a function of a real parameter, α. Therefore, the characteristic polynomial for this matrix will also depend on α, and in turn the solution of this polynomial will also depend on α. The solution of this polynomial is in fact the eigenvalue (energy level) that we are looking for. We conclude that if the real Hamiltonian matrix elements can be analytically dilated into the complex plane, also the eigenvalues can be dilated into the complex plane. Further details are given below.

Thus, one can avoid the numerical diagonalization of a non-Hermitian complex matrix and replace it with Hermitian electronic-structure calculations that yield real eigenvalues, which can be dilated into to complex plane. Similarly, one can carry out dilation into the complex plane and obtain complex properties other than the energy, such as complex dipole transitions and other properties that are calculated as expectation values. This is because the eigenvectors can be expressed as a linear combination of the HMEs, therefore we claim that they also depend on α, and can be analytically dilated into the complex plane. Thus, one can avoid the numerical calculations of the eigenvectors of the non-Hermitian complex matrix, and obtain complex energies and expectation values by analytical continuation.

Unfortunately, closed form expressions of the eigenvalues and eigenvectors as function of the HMEs are known only for small Hamiltonian matrices, i.e., less than 5 × 5. Namely, such closed form expressions are known only for Hamiltonian matrices constructed from two, three or four basis functions. Nevertheless, we can use the Padé approximant in order to get closed form expressions for the eigenvalues and expectation values for relatively large matrices, which are typically used in electronic structure calculations. However, before going into the procedure for large and finite analytical expressions (that are required for actual chemical problems) we want to establish the proof for the small dimensional matrices.

Let us show it for the simple 2 × 2 Hamiltonian matrix. The HMEs are function of a real scaling parameter η, where we chose η = αe ^iθ with θ = 0. The eigenvalues are given by E ± η = H 11 η + H 22 η 2 ± 1 2 H 11 η + H 22 η 2 + 4 H 12 η H 21 η .

It is clear that if the HMEs are known analytical functions of η, one can dilate them into the complex plane by substituting complex values for η → αe ^iθ with θ ≠ 0. Then, analytically calculate the stationary points in the complex plane (to satisfy the complex variational principle as described in Ref. [1]). The stationary points are associated with the resonance complex eigenvalues, where Re[E _±] are the energy positions and − 2Im[E _±] = Γ are the widths. Similarly, the real dipole transitions can be expressed as function of a real scaling parameter, D + , − η = Ψ + † η T d η Ψ − η .

For the 2 × 2 real Hamiltonian matrix case the vectors [Ψ_±(η)] are analytically obtained as function of the real HMEs. The components of the eigenvectors [ ψ ± ( η ) ] 1 [ ψ ± ( η ) ] 2 are given as usual by, ψ ± η 1 = H 12 η E ± − H 11 η ψ ± η 2 where ψ ± η 1 2 + ψ ± η 2 2 = 1 .

Thus the complex dipole transitions can be associated with stationary solutions in the complex plane obtained by analytic dilation of D _+,−(η) into the complex plane.

Similar expressions can be derived for the 3 × 3 and 4 × 4 cases. However, for larger matrices we shell use the Padé approximant in order to get closed form expressions for the eigenvalues and dipole transitions as a function of (a real) η. And then, carry out the analytical continuation into the complex plane and look for stationary solutions, which are associated with the resonance positions and widths or the complex dipole transitions, as detailed in Section 2 and Section 3.

In the above section we saw that complex resonance eigenvalues and dipole transitions can be obtained by analytical continuation from real eigenvalues and dipole transitions, which are obtained from standard (Hermitian and real) calculations. However, this analysis is limited to small number of basis functions (<5) for which we have closed form analytical expressions. The extension of this approach to large matrices, i.e., to large number of basis functions, requires the use of a numerical scheme to describe either { E j ( η ) } j = 1,2 , . . , N or { D j , j ′ ( η ) } j , j ′ = 1,2 , . . , N by a closed form analytical expression. One such scheme is based on the Padé approximant, in which the energy and dipoles are expressed as a rational fraction of polynomials, E j η ≈ Σ p = 0 N p a j , p η p Σ q = 0 N q b j , q η q (13) and, D j , j ′ η ≈ Σ p = 0 M p c j , j ′ , p η p Σ q = 0 M q d j , j ′ , q η q . (14)

The main question is from where do we get the data from which we will fit the coefficients in the above expression. The answer to this question is that we take it from stabilization graphs, see Section 2.1 for more details.

Other questions we need to answer are:

(Q1) How to select the initial real data points from the stabilization graphs?

These questions are answered in Section 2 (specifically Section 2.2), which discusses the RVP formalism in details.

As mentioned above, our new approach to introduce non-Hermiticity and obtain resonances is by analytic dilation of eigenvalues from the stable part of a stabilization graph into the complex plane. Since we use the Padé approximant as our analytical continuation method, we call our new approach Resonance Via Padé, RVP. This approach, unlike uniform complex scaling, uses standard, Hermitian, calculations to obtain the resonance position and width, and does not modify the Hamiltonian.

In the RVP method, as a first step, the energy spectrum is calculated using Hermitian codes as a function of a generalized box quantization parameter, E(α) where α = η (with θ = 0 and η = αe ^iθ ). For example, the eigenvalues can be calculated as a function of the number of basis functions (BFs) [45] or when finite given BFs are scaled by a real factor [47, 48, 51]. In practice, to scale the BFs by a real factor, one can use Gaussian base functions and divide the exponents of the Gaussians by the real factor α (see Eq. 11). Calculating the energy spectrum with these BFs will produce E(α). Note that in this case, α < 1 will cause the spatial distribution of the Gaussians to compress, and α > 1 will cause it to expand. Therefore, we typically employ the range: 0.6 < α < 2.0.

In such calculations, when continuously increasing α, the discrete energy levels of the quasi-continuum spectra are highly affected (i.e., lowered). Resonance states, unlike the delocalized quasi-continuum states are much less affected by small variations of α [1, 52], since resonance states are typically much more localized in the interaction region, see Figures 3.5 and 3.6 in Ref. [1]. Therefore, while quasi-continuum states significantly change as α is varied, resonance states remain relatively stable. This is why a graph portraying E(α) as function of α around the resonance energy is named a stabilization graph, as illustrated in Figure 3A.

In such a graph, an energy level crossing is expected between the highly affected delocalized quasi-continuum states and the stable resonance energy. However, since states with the same symmetry cannot cross each other in the adiabatic representation, avoided crossings are obtained in the graph. In these avoided crossings, a transition from a localized state to a delocalized quasi-continuum state occurs. Therefore, the avoided crossings are associated with branch points (BPs) in the complex energy plane (see chapter 9 in Ref. [1]). Consequently, E(α) is not an analytic function of α.

Nevertheless, while the avoided crossings correspond to a mixing between two functions: a localized function and a delocalized quasi-continuum function, the stable part of the stabilization graph, in between two avoided crossings, corresponds to a single function, localized in the interaction region. Therefore, the stable part of the stabilization is expected to be locally analytic, while the avoided crossings are expected to be non-analytic. Thus, the stable regions contain all the relevant information for analytical continuation into the complex plane. A numerical proof for this point is given in Refs. [52, 53] and the figures within.

This is how the RVP approach avoids the non-analytic parts in the stabilization graph. In this approach, a single energy level, obtained from the stable part of the stabilization graph, is analytically dilated using the Padé approximant. Namely, the stable part of the stabilization graph, which has a smaller slope than the unbound energies [46], is fitted as a function of a real scaling parameter to a ratio between two polynomials (like Eq. 13): E α = P α Q α . (15) where P(α) and Q(α) are polynomial functions of a real scaling parameter, α. As the focus of the analytic continuation scheme is not on the avoided crossing regions, but rather on the stable part of the stabilization graph, excellent results are obtained [50, 53].

As previously explained, the stable part of the stabilization graph is expected to be locally analytic. Yet this is not all, this stable region is also known to contain information on the resonance lifetime. It has already been shown by Hazi and Taylor. [46], that the slope of the stable region is related to the width of the resonance. That is, as the resonance width increases the slope of the stable part increases. Furthermore, as shown in Section 3.4 in Ref. [1] (and Figures 3.4–3.8 wherein), one can even estimate the resonance width by analysing the localized function that is associated with the stable region’s eigenvalues. In this case, the width is proportional to the square of the ratio between the normalized amplitude of the function out of the interacting region and in the interaction region. So, the stable region in a stabilization graph contains enough information on the resonance lifetime and all the relevant information for analytic dilation into the complex plane.

Therefore, it is clear why this region should be used if one wants to carry analytical continuation to the complex plane and gain insight on the resonance lifetime. In other words, the answer to Q1 in Section 1.5 is clear - the initial real data points from the stabilization graphs one needs to use is the data from the stable part of the stabilization graph.

Indeed, as a second step in the RVP method, data from the stable region is analytically dilated into the complex plane using the Padé approximant (Eq. 15) in order to locate stationary points (SPs), resonances. In Ref. [53], an analytical path from the stable region towards a complex stationary point is shown. This path goes between the BPs and bypass them. This is an additional proof that using Padé, one can always remain on an analytic path in the complex plane that goes towards a stationary point. Notice that the existence of such a path results from the use of a finite basis set, which are always used in any electronic-structure calculation [53].

In practice, within the RVP method, an analytic Padé function is fitted using the Schlessinger point method [54] to data from the stable region. The Schlessinger point method requires a set of M data points (α _i ) and their corresponding eigenvalues (E _i ), and then the Schlessinger truncated continued fraction assumes the following form: C M α = E η 1 1 + z 1 α − α 1 1 + z 2 α − α 2 ⋮ z M − 1 α − α M − 1 , (16) where the z _i coefficients need to be determined recursively such that C M α i = E α i , i = 1,2 , … , M . (17)

This truncated continued fraction can be transformed to a Padé like form (Eq. 15). Thus, by choosing to use the Schlessinger point method for the Padé approximant, one obtains automatically from the data points, all the information on the Padé function, namely the answers to Q2 and Q3 in Section 1.5 above. Q2 deals with the degree of the polynomials in the Padé function, and Q3 deals with their coefficients. Clearly, both are determined by the data set itself satisfying Eq. 17.

Once the z _i coefficients are determined, and Eq. 16 is completed, an analytic continuation into the complex plane is performed. This is done by substituting a complex η, instead of α, i.e. η = αe ^iθ with θ ≠ 0. Then, SPs, resonances, can be identified by generating α- and θ-trajectories and looking for cusps in the complex plane [52, 55]. Alternatively, SPs can be identified by solving the algebraic equation d E d η = 0 [52]. While solving an algebraic equation is easier than looking for cusps in the complex plane, the SPs found by solving the algebraic equation are not necessarily associated with resonances [52]. This means that if we solve the algebraic equation some SPs are unphysical. Therefore, to identify the resonance energy and lifetime, we use a clusterization technique described in Ref. [50] as the final step of the RVP method. In practice, we generate Eq. 16 for different input sets, where all sets are taken from the stable region. This way we get a large number of SPs by solving the algebraic equation for each set. The SPs depend on the input points chosen for the analytical dilation, where as mentioned, some SPs are also unphysical. The physical SPs should not depend strongly on the variation of the input data, unlike the unphysical ones. Therefore, we examined the SPs by statistical distribution and look for clusters of SPs obtained from different input data. The mean value of the cluster is reported as the resonance complex energy. In this way, we finally answer Q4 in Section 1.5 above: We understand that we locate the stationary solution in the complex plane, by solving the algebraic equation d E d η = 0 for many input sets, and by using a statistical clusterization technique.

To sum up, by using the RVP method, analytic continuation of a single eigenvalue level into the complex plane can be employed, provided that the input data is taken from the stable region of a stabilization graph since it is a local analytic region. Therefore, we can divide the RVP method to three steps. In the first step the stable region data is fitted into a Padé/Schlessinger function and dilated into the complex plane. An analytic path from the stable region to the complex SPs, which avoids any of the BPs, exist when a finite basis set is employed. In the second step, the SPs are located using the algebraic equation d E d η = 0 , for different input data sets. Then, as a final step, the clusterization technique locates the physical complex resonance values out of the total SPs.

The calculated RVP complex potential energy surfaces (CPESs) that are presented below play a crucial role in the interpretation and analysis of the reaction rates (RRs) measured in cold molecular collision.