Variational Solutions for Resonances by a Finite-Difference Grid Method

We demonstrate that the finite difference grid method (FDM) can be simply modified to satisfy the variational principle and enable calculations of both real and complex poles of the scattering matrix. These complex poles are known as resonances and provide the energies and inverse lifetimes of the system under study (e.g., molecules) in metastable states. This approach allows incorporating finite grid methods in the study of resonance phenomena in chemistry. Possible applications include the calculation of electronic autoionization resonances which occur when ionization takes place as the bond lengths of the molecule are varied. Alternatively, the method can be applied to calculate nuclear predissociation resonances which are associated with activated complexes with finite lifetimes.


Introduction
In 1978 Frank Weinhold together with Phil Certain and Nimrod Moiseyev, in their work on the complex variational principle (a stationary point rather than an upper bound as in Hermitian QM), paved the way for the use of electronic structure computational algorithms to metastable (resonance) states [1]. In this framework, the energies and inverse lifetimes of atoms and molecules are associated with the real and imaginary parts of the complex eigenvalues of non-Hermitian Hamiltonians, respectively.
Recently, there is increasing interest among chemists concerning the use of non-Hermitian QM for the calculations of molecular resonances. For example, the most recent developments of the QCHEM quantum chemistry package enable calculations of shape type and Feshbach molecular resonances [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. In addition, such methods have been employed in the study of RNA bases [19,20]. It is clear by now that in order to increase the accuracy of ab-initio calculations, associated with resonances, grid methods should be incorporated into the Gaussian basis set. See for example the hybrid Gaussian and b-spline method recently developed by Fernando Martin and his group [21,22]. Specifically, this will enable to improve the description of the outgoing ionized electrons. However, the standard finite difference method (FDM), that enables calculation of the bound states of a quantum system, is not applicable for the calculations of the complex poles associated with metastable (resonance) states. The origin of this failure can be traced back to the fact that the standard FDM does not satisfy the variational principle.
Here we show how a simple change in the selection of the grid points in FDM leads to a variational principle and enables calculation of both real and complex poles of the scattering matrix. This approach opens the gate to evaluate the resonances by FDM for atoms and molecules as well as mesoscopic systems. Illustrative numerical examples will be given for a "toy" model, which was first presented in a paper published together with Frank Weinhold many years ago and is used until now as a model for testing new algorithms for calculating resonances [1].
Numerical approaches for the analysis of physical systems can be classified into two prominent categories: grid and basis set approaches [23][24][25][26][27][28][29]. The basis set methods are equivalent to the use of an approximate representation of the identity operator. As a result, they provide upper bounds to the eigenvalues. Contrastly, in the grid based methods one represents the continuous space by a quantized finite number of grid points. These methods exhibit fast processing time, however, they generally do not provide an upper bound to the eigenvalues [30][31][32][33].
The standard grid method is the traditional finite difference method, which is abundantly used in the solution of second order partial differential equations, e.g., in the study of heat transfer problems [34], as well as in solving the Maxwell [35,36] and Schrödinger equations [37]. The crucial limitation of the standard FDM is that the convergence of the numerical results requires refining the grid spacing (mesh), which in turn increases the amount of storage and calculation. An important improvement of the accuracy and stability of the FDM has been recently described in Ref. [38] by combining two high-order exponential time differencing precise integration methods (PIMs) with a spatially global sixth-order compact finite difference scheme (CFDS). In addition, by modifying the representation of the Laplacian operator one can obtain a rigorous upper bound estimate of the true kinetic energy [39].
The first goal of this paper is to show that upper bounds to the spectrum of any given Hamiltonian can be obtained without modifying the representaion of the Laplacian and by using the same set of coupled equations as are used in the standard FDM. We refer to the proposed method as the "present" variational FDM, while the common approach is termed the "standard" non-variational method. The standard FDM typically converges to the exact spectrum from below (i.e., non-variational), this is attributed to the fact that the obtained spectrum of the kinetic energy operator in the standard FDM serves as a lower bound to the exact kinetic energies, see Figure 1 and also Ref. [39]. Note however that this characteristic behaviour is not true for any potential.
Building upon the Hylleraas Undheim and MacDonald theorem we prove that the proposed present FDM satisfies a variational principle with respect to the accurate solution within the finite box approximation. The variational principle guarantees the stability of the proposed scheme as the number of the grid points are increased. The stability of the present FDM calculations is obtained by holding the grid spacing to be as small as possible and constant, while increasing the number of grid points.
We first focus on the calculation of the bound discrete states of Hermitian Hamiltonians. Following, we show how the the present FDM can be utilized to calculate the energies and widths (inverse lifetimes) of mestasbale states states, embedded in the continuous part of the spectrum (so called resonances), which are associated with the poles of the scattering matrix.
We introduce the finite box quantization condition, assuming that this restriction does not serve as a limitation to obtain the bound state spectrum in the desired accuracy. That is, the exact result is considered as the result obtained by fixing the spatial range of the system and infinitely increasing the precision of the calculation. Physically, this is motivated by the fact that any realistic computation is conducted by using a finite number of grid points or basis states, i.e., finite size computers. Moreover, any realistic measurement has a corresponding fundamental uncertainty. Therefore, one can replace infinite space by a box of finite dimension, without practically effecting the physical description.
Within the finite-box approximation, the considered exact solution is determined by two parameters: the box size L max and the maximum number of grid points N max . In the present approach the grid spacing is defined by δx = L max /N max , where the box-size L is varied with the number of grid points N, i.e., for L = δxN. This should be distinguished from the standard method, where the grid spacing ∆x = L max /N varies with the number of grid points where N ≤ N max .
The paper is organized as follows. First, we describe the two FDMs procedures, leading to the spectrum of the Hamiltonian under study. We then provide a proof that the proposed FDM produces an upper bound to the spectrum of the Hamiltonian within the box quantization condition. Following, we present the numerical results for the calculation of the bound and metastable states (resonances), and compare to the exact results. Finally, we conclude by emphasizing the generality of our approach.

Methodology
When conducting a numerical calculation utilizing a grid based method, the Hamiltonian operatorĤ is represented by a N × N dimensional matrix, where N is the number of grid points and T and V are matrix representations of the kinetic and potential energies operators. The matrix T is calculated by utilizing m (m = 2l + 1; l = 1, 2 . . .) grid points to evaluate the Laplacian (second-order derivative). Specifically, we employ the central finite grid method, where m values of the wavefunction around the grid point x i (identical number on both sides) are used to approximate the second order derivative of the wave function (cf. Appendix A for further details). This leads to the following relation where w j = A −1 3,j+l+1 , and A is a m × m matrix, with elements A j+l+1,k+1 = j k /k! and j ∈ [−l, l]. This relation determines the matrix elements of the (2l + 1)-diagonal matrix T whose elements are given by while zero otherwise and i = 1, 2, . . . , N. Hereh is Planck's constant and µ is the particle's mass. Table 1 gives the coefficients for different values of m. We emphasize that the coefficients are symmetric w j = w −j ; this entails that the kinetic energy matrix is symmetric and real, which in turn implies that it is positive definite.  (2) for different values of m. Here we present only w j for j ≥ 0 as the weights remain symmetric for every j, e.g., The number of grid points included in the calculation of T (m) has a significant effect on the eigenvalues of the matrix. This can be witnessed by analyzing the convergence to the exact solution with increasing n, see Appendix B for a graphical representation.
Alternatively, the present (variational) FDM can also provide an upper bound to the spectrum of the Hamiltonian. For a properly defined grid, when the grid-difference is held fixed δx ≡ ∆x(N max ), and the box-size is increased with N: L(N) = δxN, the numerical result converges to the exact spectrum from above. In this case, as seen in Figure 1, the eigenvalues of T provide a curve which approaches the parabolic function y(n) = c(n/L max ) 2 , with L max = L(N max ), from above. As N approaches the maximal value N max , the obtained spectrum approaches the same converged numerical result of the kinetic energy, leading to the converged Hamiltonian spectrum. The kinetic energy matrix is calculated by using m = 3 (see Appendix A) and therefore is a tri-diagonal matrix. It is evident that the standard method provides lower bounds for the kinetic energy spectrum for every n ≥ 0 as N is increased (N = 201, 251, . . .), whereas the present (variational) FDM, first presented here, provides upper bounds to the exact values. In this case, the standard FDM converges faster than the present FDM but it is not a variational method.

The Variational Principle for Finite Difference Method
We consider a system confined in a finite box which is represented by Hamiltonian H. The box size, utilized in the calculation, is chosen so to not limit the accuracy of the calculated eigenstates. This is a common approximation, which is implicitly included in any numerical calculation, including all quantum chemistry packages used to obtain the electronic spectrum. In such calculations, the molecular HamiltonianĤ is replaced by a finite dimensional matrix H. Similarly, in the FDM we limit the 1D space to a < x < b. The exact Hamiltonian under study within the box quantization framework is obtained by Based on the Hylleraas, Undheim and MacDonald theorem [40][41][42], we prove that the eigenvalues of the N-grid point representation matrix of the Hamiltonian, H(N), serves as an upper bound to the exact spectrum. The N by N matrix H(N) satisfies the following eigenvalue equation where the columns of C(N) are the eigenstates of H(N) and E diag is a diagonal matrix containing the corresponding eigenvalues. Clearly, the N + 1 by N + 1 matrix H(N + 1) satisfies a similar eigenvalue equation. This matrix can be expressed in terms of H(N) matrix as follows where M = (H 1,N+1 , . . . , H N,N+1 ) T with H ij are the corresponding matrix elements of H(N + 1). Note that in the present case, where the system is represented by a grid, M never vanishes.
After some algebraic manipulations one obtains the relation where ε(N + 1) is one of the eigenvalues of H(N + 1). Equation (6) can be now solved by replacing ε(N + 1) by a parameter x and plotting both sides of the equation as a function of x. The intersection between the two curves are values of x = ε(N + 1) for which Equation (6) is satisfied.
Poles of Equation (6) are obtained whenever ε(N + 1) is one of the eigenvalues of H(N), i.e., one of the elements on the diagonal of E diag (N). In Figure 2 a schematic representation of the left hand side (l.h.s) and r.h.s of Equation (6) are plotted as function of the parameter x. By observing Figure 2 it is evident that This equation shows that the eigenvalues converge from above. Hence, the eigenvalues obtained using a finite number of grid points upper-bound the exact eigenvalues. The strict inequality between the eigenvalues of Equation (7)  To demonstrate how the two FDM schemes can be combined together to evaluate the system spectrum, we compare the FDM results to the analytical solution for two cases: the harmonic and Rosen-Morse potentials. The harmonic potential, V HO (x) = 1 2 µω 2 x 2 , includes an infinite number of bounded states with energies E n =hω n + 1 2 , where n = 0, 1, 2, . . ., µ is the particle mass and ω is the oscillator frequency. In contrast, the Rosen-Morse has finite number of bounded states n max with energies where n ≤ n max , and the potential is of the form V RM = −V 0 /cosh 2 (ax) [43]. The calculation is performed by the following procedure: First, we evaluate the maximum box size L max utilizing a semi-classical approximation. The semi-classical bound state function is well described when the box quantization condition is imposed on the quantum solution, such that |A( where the eigenenergies of interest lie in the range [min x (V(x)), E max ]. This evaluation is equivalent to employing the WKB method in order to recast the wavefunction in an exponential form [44][45][46]. Technically, this approximation is valid only for large action relative toh and smooth potentials, nevertheless, this condition is sufficient to evaluate L max even beyond this regime. In our calculations we take A ≈ 10 −7 . For higher dimensional space the box-size should be evaluated in a similar way, by choosing the spatial coordinates according to the classical turning points of the potential in the energy range under study.
We now compare the eigenenergies of the two FDM schemes for a varying number of grid points N. The standard procedure (L = const), typically, produces a lower bound, while keeping the grid density constant with increasing grid size gives an upper bound to the spectrum. This can be observed in Figures 3 and 4, which present the energy error, error(E n ) = E numerical n − E exact n /|E exact n |, as a function of quantum number n for the two potentials. The two cases demonstrate the varying convergence behaviour. In the case of the harmonic potential, the standard FDM shows a faster convergence from below relative to the present method. In contrast, the later method shows a rapid convergence for the Rosen-Morse potential from above. This demonstrates the utility of applying both methods, and combining them to evaluate the exact spectrum.

Illustrative Numerical Examples for the Calculations of Energies and Widths of Resonances
We now apply the present FDM to calculate the spectrum of a model potential Such a potential has been used to study new computational algorithms for calculating the energies and widths (inverse lifetimes) of shape type resonances [1]. The spectrum is characterized by a single bound state and metastable states with higher energies. In addition, in Ref. [47] this model was employed to calculate upper and lower bounds to the complex decay poles of the scattering matrix (resonances). In the following, the two finite different methods are applied to solve for the spectrum of H r = T + V r and the eigenenergies are plotted as a function of the number of grid points in Figure 5. The convergence of the standard (non-variational) method is characterized by non-intersecting lines, Figure 5a. As a result, this plot does not indicate which states are metastable. In contrast, the present variational FDM produces an informative picture of a typical stabilization graph, Figure 5b, allowing to distinguish the resonances from the other states in the quasi-discrete spectrum of the continuum. The possibility to isolate the resonances from the other states in the quasi-continuum can be utilized to calculate the resonance widths. See for example, the calculation of resonance widths (inverse lifetimes) and energies for L 2 methods [48], spherical-box quantization [49] and from the analytical continuation of real stabilization graphs, utilizing a Gaussian basis [19,20,. This is the primary result of our work, demonstrating that resonances can be calculated from stabilization graphs obtained by utilizing grid methods and not only by applying basis sets, as previously done. finite difference methods. The plots show the spectrum, E n ; n = 1, 2, . . ., as a function of the number of grid points N for a model potential V r (x), Equation (8). The full dark points on the right hand side indicate the exact values of the bound and the two nearest resonance energies as calculated by the uniform complex scaling approach. Both methods lead to accurate values as N increases, however, in the standard method, Panel (a), does not allow to distinguish between the resonances and the other quasi-discrete continuum states. In contrast, the variational principle of the present FDM leads to a typical stabilization graph. The stability of the metastable states in the continuum enables identifying them. Their positions (i.e., energies) is determined by the stable region. Model parameters: δx = L max /N max ≈ 8.5 × 10 −2 , L max = 60 and N max = 700. These numerical parameters lead to an error of 10 −9 in the value of bound state, as well as 10 −5 and 10 −2 in the values of the two first metastable states, respectively.
In order to obtain the resonance energies (in an improved accuracy) and the corresponding resonance widths we repeat on the FDM calculations with a uniformly rotated coordinate in the complex plane. Formally, the procedure maps the x coordinates to {x i → x i exp(iθ)} i=1,2,...,N , leading to complex eigenvalues. The real part of the eigenvalues that are invariant under the mapping (invariant under a change of θ) are the resonance energies (positions). While the resonance widths are associated with the imaginary parts of the complex eigenvalues multiplied by −2 (E n → ε n − iΓ/2).
For a formal justification for calculating the resonances by a rotation of the coordinated in the complex plane see the text book on non-Hermitian quantum mechanics and the references therein [71]. Specifically, the works of Aguilar and Combes [72], and Balslev and Combes [73], should be highlighted along with the work of Barry Simon that put the computations of resonances by the complex scaling transformations on a solid mathematical ground [74,75].
The results presented in Figure 6 were obtained by following the complex eigenvalues which their real parts are closest to the resonance values obtained in Figure 5 (θ = 0). This is a crucial property of the present method, as it enables following the convergence of the complex poles as the number of grid points N is increased. This cannot be achieved by the standard FDM. The stability which results from introducing absorbing boundary condition (here obtained by employing complex scaling) has been observed before for the helium resonances in Ref. [76], utilizing Hylleraas basis functions. Therefore, it is expected that similar stabilization graphs will be obtained by grid methods in combination with complex absorbing potentials.

Concluding Remarks
We show that by fixing the grid spacing in the finite difference method one obtains a variational principle within the finite box approximation. This property allows obtaining a stabilization graphs for the spectrum, which produces an accurate estimation for both the bound ground, excited states as well as the positions of narrow resonances. The stabilization graphs obtained by the presented FDM enable one to distinguish between the metastable (resonance) states that are localized in the interaction region and the other states in the continuum. The later states are not localized and have large amplitudes outside the interaction region. To obtain the resonance width (inverse lifetimes), we preformed a rotation of the grid points in the complex plane. This procedure also increases the accuracy of the the resonance positions. The calculation accuracy is determined by the grid spacing. As the grid spacing decreases the accuracy increases at the expense of increasing the number of grid points required for convergence.
Overall, the present study demonstrates how a simple change in a known numerical method (FDM in our case) might increase the broadness of its application. Our results suggest that grid methods should be added to the Gaussian basis set in the electronic structure quantum chemistry packages. This will allow introducing complex absorbing boundary conditions in the non-interacting region of the Hamiltonian. This approach paves the way for the calculations of the complex poles of the scattering matrix. Such poles can be associated with the peaks in measured reaction rates in cold molecular collisions.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

QM
Quantum mechanics FDM Finite difference method

Appendix A. Derivation of the Representation of the Laplacian by the Standard Finite Difference Method
The kinetic energy operator T is presented by a N × N matrix, which employs m (where m = 3, 5, . . .) grid points to approximate the representation of the Laplacian operator at each point. This approximation results in a discretized spectrum of the kinetic energy operator, given by E n = c(n/L) 2 , where n ∈ N is the quantum number, L = x N − x 1 is the size of the box which discretizes the kinetic energy spectrum and the proportionality constant c is problem dependent. In the solution of the time-independent Schrödinger equation the proportionality constant equals to π 2h2 /(2µ). Commonly, in a grid representation, the matrix representing the potential energy is diagonal with values V(x i ), where x i denotes the coordinate of the i th nodal point.
For the sake of clarity we give below a short description of the derivation of the approximate matrix representation of the Laplacian operator. Consider a 1-dimensional evenly spaced grid made up of N nodal points with a total length L. The spacing between adjacent nodal points is given by ∆x = L N−1 . We wish to approximate the second order derivative of the wave function ψ(x) at x = x i , compactly denoted as ψ i . For this purpose, we write the truncated Taylor series expansion around x i using j nodal steps, explicitly written as Notice that the first and last elements of ψ are ψ i−l and ψ i+l , respectively. The nodal values and its derivatives are related through the matrix A m×m , with elements A j+l+1,k+1 = j k /k!: By inverting Equation (A3), we isolate the second order derivative, which is proportionate to the third element of ψ (D) . This leads to a linear combination of the nodal values, with weights w j = A −1 3,j+l+1 . The derivative is then explicitly written as (∆x) 2 d 2 ψ dx 2

Appendix B. Convergence of the Kinetic Energy Representation by the Standard Finite Difference Method
The eigenvalues of T converge towards the analytical result with an increase in the number of grid points included in the calculation m. The different matrices T(m) are diagonalized for a constant number of grid points N, to obtain the eigenenergies E n for different values of m. Figure A1 presents the eigenvalues of the matrices T(m) in increasing order. As the m increases the representation of the kinetic energy operator becomes more accurate and the corresponding eigenenergy curves converge to the exact result. Therefore, the value of m has a major influence on the derivation of the upper bound within a specific potential.