Finite Element Method of High-Order Accuracy for Solving Two Dimensional Elliptic Boundary-Value Problems of Two and Three Identical Atoms in a Line

We considered models of three identical atoms in a line with molecular pair interactions and diatomic molecule scattered by an atom or tunneling through potential barriers. The models are formulated as 2D elliptic boundary-value problems (BVPs) in the Jacobi and polar coordinates. The BVP in Jacobi coordinates solved by finite element method of high-order accuracy for discrete spectrums of models under consideration. To solve the scattering problems the BVP in polar coordinates are reduced by means of Kantorovich method to a system of second-order ordinary differential equations with respect to the radial variable using the expansion of the desired solutions in the set of angular basis functions that depend on the radial variable as a parameter. The efficiency of the elaborated method, algorithms and programs is demonstrated by benchmark calculations of the resonance scattering, metastable and bound states of the considered models and also by a comparison of results for bound states of the three atomic system in the framework of direct solving the BVP by FEM and Kantorovich reduction.


Introduction
With the development of modern computing power, there are more possibilities for setting and numerically solving multidimensional boundary-value problems with high accuracy.For this, new numerical methods of high accuracy order are being developed.When reducing the boundary value problem to an algebraic one in the finite element method (FEM) of order , one of the problems is the calculation of integrals on a finite element (we consider only simplicial finite elements) containing the products of two basis functions of Lagrange or Hermite interpolation polynomials of order  by the coefficients for the unknown functions [1][2][3].The FEM schemes and algorithms for solving the d-dimensional BVP are presented in [4][5][6].The algorithms for constructing 2-dimensional fully symmetric Gaussian quadratures are presented in [7] while d-dimensional ones will be published elsewhere.
The aim of this paper is to demonstrate the efficiency of the elaborated algorithms and program complexes KANTBP 4M, KANTBP 3, ODPEVP [8][9][10] by benchmark calculations of the resonance scattering below the dissociation threshold, metastable and bound states of the considered models and also by a comparison of results for bound states of the three atomic system in the framework of direct solving BVP by 2D FEM [3,7] and Kantorovich method (KM).
We elaborated also algorithms for calculating the asymptotic parametric angular functions, the effective potentials and the fundamental solutions of the SODEs in the form of expansions by inverse powers of radial variable [11] and apply them to the construction of the asymptotic states of the triatomic scattering problem, because in three-atomic systems at large values of the hyperradial variable the effective potentials of SODEs have the form of expansions in inverse powers of hyperradius even for short-range potentials of pair interactions.
The paper is organized as follows.In Section 2 we formulate the 2D BVP for the dimer and trimer models.In Section 3 the 2D BVP in polar coordinates is reduced to the SODEs.In Sections 4, 5 and 6 we demonstrate the efficiency of the elaborated technique by test calculations of the resonance scattering of a diatomic molecule on an atom, the metastable and bound states of three atomic system, and the resonance tunneling of a diatomic molecule through a Gaussian barrier and the metastable states, respectively.In Conclusion the results and perspectives are discussed.

Setting of the Problem
Consider a model of three identical particles (trimer) on a straight line with the masses   =  and the coordinates x ∈ R 1 ,  = 1, 2, 3, of one-dimensional Euclidian space R 1 coupled via the pair short-range potential Ṽ (|x  − x |), ,  = 1, 2, 3. Performing the change of variables corresponding to the cyclic permutation (, , ) = (1, 2, 3) of the permutation group  3 [12]: we get three pairs (︀  () ,  () )︀  of the scaled Jacobi variables.In the center-of-mass reference frame of the configuration space are connected by the relations: where the angles With the variable parameter  in T  () they and simply change the sign of pairs of Jacobi coordinates in Eq. ( 1), In Fig. 1 (a) the hyperplane (, )  ∈ Ω  is separated by three axes  (12)3 ,  (31)2 and  (23)1 into six sectors that together with three orthogonal axes  12 ,  31 ,  23 describe six channels of the scattering problem for three identical particles.Our choice is determined by the parameterization of the pairs (, )  = ( cos ,  sin )  , where the angle  ∈ [0, 2) is counted counterclockwise from the axis  (12)3 , for which  = 0.
As a result, we arrive at the Schrödinger equation for the wave function Ψ(, ) corresponding to the total energy  of the triatomic system in Jacobi coordinates (, )  = (  ,   )  ∈ Ω  : ✭❛✮ ✭❜✮ Here the potential function for the trimer with the pair potentials (below this case is referred to as Task 2, for example, see Fig. 1 or the potential function for a dimer in the field of a barrier potential (below this case is referred to as Task 3, see Fig. 1 is an even function with respect to the straight line  = 0 (i.e., x1 = x2 ), which allows one to consider the solutions of the problem in the half-plane  0. The use of the Dirichlet or Neumann boundary condition at  = 0 allows one to obtain the solutions, symmetric and antisymmetric with respect to the permutation of two particles.If the pair potential possesses a high peak at the pair collision point, then the solution of the problem in the vicinity of  = 0 is exponentially small and can be considered in the half-plane   min .
In this case the Neumann or Dirichlet boundary condition imposed at  min yields only a minor contribution to the solution.The equation, describing the diatomic molecular subsystem (dimer) (below this case is referred to as Task 1 ), has the form We assume the energy spectrum of the dimer to consist of the discrete part with a finite number  0 1 of bound states with the eigenfunctions   (),  = 1,  0 and the eigenvalues ε = −|ε  |< 0, and the continuous part with the eigenvalues ε > 0 and the corresponding eigenfunctions  ε().In certain cases the eigenfunctions of the continuous spectrum are approximated by the eigenfunctions of pseudostates of the discrete spectrum ε > 0,  = 1 +  0 , . . .calculated in the sufficiently large but finite interval  ∈ [ min ,  max ].
Using the above pairs of the basis functions (15), we rewrite the expansion (11) in the -representation where  = ,  or  = ,  for the unknown functions     () and     () in the ()representation, related to the functions     () and     () in the ()-representation as The averaging of Eq. ( 8) with the basis functions in -representation (17) yields the system of coupled ODEs with determined by the hyperspherical parameterization of the two-dimensional configuration space Here the potential curves (terms)   () are eigenvalues of the BVP (12) and the effective potentials (EPs)   () = −  (),   () =   () are expressed as integrals calculated in the reduced intervals  ∈ [0, 2  ] using the above ,  symmetry: For Task 3 the effective potentials Ŵ () =   () +    () are sums of   (), calculated using the potential curves and the parametric basis functions of Task 1, and the matrix elements    () of the barrier potentials As an example, we calculated with the required accuracy the parametric basis functions of BVP (12) and the effective potentials (21) for the models of Be 2 dimer and Be 3 trimer in collinear configuration using the FEM implemented in the program ODPEVP [10].The results of calculation on the grid Ω  [1.8/,   ] = {1.8/(24)3/(10)4/(5)5/(10) } for   = /2 for Be 2 dimer and   = /6 for Be 3 trimer.Using the obtained result in the uncoupled ()-representation and the transformation matrix  from Eq. ( 18), we can rewrite the system of ODEs in the coupled ()-representation with the effective potentials.
In Section 4 one can see that the ()-representation provides the required compatibility of the solutions of Eqs.(19) with the asymptotic boundary conditions of the scattering problem on the full axis and its half-axis.

Asymptotic Expressions of Scattering, Metastable and Bound States
The general solution   of the system of ODEs in the open channels   = 1, . . .,   is determined by a linear combination of the fundamental solutions  −  ′  (   ) and  +  ′  (   ) calculated using Eqs.( 19)- (21) following from [11] with the leading terms of the Hankel functions of the first and the second kind [17] 1/2 (   ) and (1) 1/2 (   ) below the dissociation threshold at  < 0 in the form of incoming and outgoing waves The components of the radial asymptotic solutions     () of the scattering problem in the open channels   = 1, . . .,   have the form while in the closed channel   =   + 1, . . .,  the asymptotic solutions    () are determined by the fundamental solutions  +  ′  (  ) with the leading term of Kelvin functions [17]  1/2 (   ) for the decaying waves These asymptotic solutions  () = {   ()}    =1 = {{   ()}  =1 }    =1 are used to have the conventional asymptotic boundary conditions for the components of the numerical solution    () of the system of ODEs (19) The scattering problem ( 3)-( 6) with the asymptotic boundary conditions ( 25) and ( 26) is reduced to a boundary-value problem for the set of close-coupled equations (19) with the boundary conditions at  =  min and  =  max : where is the required  ×  numerical matrix solution.These matrices and the
For metastable states the even and odd eigenfunctions obey the boundary conditions of the third kind (28), where the matrices ℛ( max ) = diag(ℛ( max )) depend on the complex energy eigenvalue  ≡   = Re   +  Im   , Im   < 0 sought for, and are expressed as [19] since the asymptotic solutions of this problem contain only outgoing waves in the open channels   = 1, . . .,   and closed ones   =   +1, . . .,  .In this case the eigenfunctions obey the orthogonality and normalization conditions For bound states the even () and odd () eigenfunctions obey the boundary conditions (28), with the matrices ℛ( max ) = diag(ℛ( max )) = 0.In this case the eigenfunctions obey the orthogonality and normalization conditions Taking the property ( 15) and ( 16) of the quasiangular parametric basis functions and the effective potentials (23) into account, we express the S-matrix in the () representation on the full axis ℛ 1 via the matrix   (24), (25) calculated on the half-axis ℛ 1 + .The matrix S is a unitary and symmetric scattering matrix consisting of the matrices S  , S  and S  = S   of the dimension   ×   determined by the relations S  = S  = (S +1 + S −1 )/2, S  = S   = (S +1 − S −1 )/2, where S +1 ≡ S  and S −1 ≡ S  are the matrices from (25).Here I is the unit matrix with the dimension 2  × 2  , S  and S  corresponds to the elastic scattering processes (in the considered case of 1D scattering it means reflection) of the dimer ()(or ()) on the atom  (or ):  + () →  + (), or () +  → () + , and S  and the matrices S  correspond to the inelastic rearrangement scattering processes (in the case of 1D scattering it means transmission)  + () → () +  or () +  →  + (), for which the conventional relations between inelastic and elastic scattering below breakup threshold at  < 0 follow from (32) that provide conservation of the Wronskian, , where I is the unit   ×   matrix.For the scattering of the dimer () by the potential barriers, considered on the full axis, the matrix S is the 2  × 2  scattering matrix (32) read as similar to [19] where

Bound, Metastable and Scattering States of the trimer
For the considered models, the eigenvalues and the hyperradial components of 2D eigenfunctions of the BVP for the set of ODEs (19) with Dirichlet boundary conditions were calculated with the predetermined accuracy using the FEM implemented in the KANTBP 4M program [8].
The comparison of sets of total and binding energies of g and u bound states of the trimer Be 3 calculated by KM and 2D FEM is presented in Table 1 and the corresponding eigenfunctions (11) are shown in Fig. 2. One can see that results obtained by KM and 2D FEM [3,7] are in agreement with an accuracy of the order 0.1 Å−2 .These solutions, i.e. the real eigenvalues and the corresponding eigenfunctions, were used as an initial approximation in the continuous analogue of Newton's method [21] with additional condition (F  |F  ′ ) = 0 to calculate the metastable states of the trimer Be 3 on the same finite-element grid.The corresponding problem of the dimer scattering on an atom was solved on the same grid and  = 15.
The calculated complex energies of the metastable states    ≡   = Re   +  Im   for the trimer are presented in Table 2.
Table 2 The sets of the first resonance energy values  at which the minimum of the transmission coefficient is achieved, the number  of the threshold   , the real and imaginary part of the complex energy eigenvalues   = Re    +  Im    in Å−2 of the even  and odd  metastable states of Be 3 numbered by the index  calculated with  = 15 equations ( 19)

Metastable and Scattering States of the Dimer Tunneling
The metastable states of the dimer Be 2 tunneling through Gaussian barriers of Task 3 were calculated for BVP calculated with  = 15 equations in the system (19) with matrix elements of potential barrier on the finite element grid Ω  = { min =1.81(12)4.21(15) max = 7.51} with the fifth-order Hermitian elements ( = 5).The corresponding problem of a dimer tunneling through the barriers was solved on the same grid.
The corresponding algebraic eigenvalue problem for metastable states was solved using the above mentioned continuous analogue of Newton's method.As the initial approximation the real eigenvalues and the corresponding orthonormalized eigenfunctions (31) were used.They were found as a result of solving the bound state problem with ℛ(  ) = 0 on the grid Ω  = { min =1.81(12) max = 4.21}.The complex values of energy of the metastable states    ≡   = Re   +  Im   for the dimer Be 2 tunneling through the Gaussian barriers, are presented in Fig. 6.
These metastable states are responsible for the resonance values of energy, corresponding to the maximal transmission coefficient, i.e., the quantum transparency of the potential barriers (see Figs. 1 (b) and 3 (b)), i.e., the shape resonances.The position of peaks presented in Fig. 6. is seen to be in quantitative agreement with the real part Re(), and the geometric halfwidth of the | | 2 11 () peaks agrees by the order of magnitude with the imaginary part Γ = −2 Im() of the complex energy eigenvalues  = Re() +  Im() of the metastable states.The obtained complex energy values corresponding to the resonance values of energy in the first open channel are in good agreement with the ones calculated analytically in the model of a rigid diatomic molecule with Morse potential tunneling trough the Gaussian barrier at the same values of parameters [13].From Fig. 6 one can see that as the energy of the initial excited state increases, the transmission peaks demonstrate a shift towards higher energies, the set of peak positions keeping approximately the same as for the transitions from the ground state and the peaks just replace each other.For example, the left epure shows that the positions of the 13th and 14th peaks for transitions from the first state coincide with the positions of the 1st and 2nd peaks for the transitions from the second state, while the right epure shows that the positions of the 25th and 26th peaks for transitions from the first state coincide with the positions of the 13th and 14th peaks for transitions from the second state and with the positions of the 1st and 2nd peaks for the transitions from the third state.

Conclusion
The model for three atomic beryllium system in a straight line was formulated as a 2D boundary-value problem for the Schrödinger equation in Jacobi and polar coordinates.Using the Kantorovich expansions this problem has been reduced to the boundary-value problem for a set of second-order ordinary differential equations.
The efficiency of the elaborated method, algorithms and programs has been demonstrated by benchmark calculations of the resonance scattering, metastable and bound states of the considered models and also by a comparison of results for bound states of the three atomic system in the framework of direct solving BVP by FEM and Kantorovich reduction.
The effects of resonant quantum transmission of diatomic molecule through the potential barrier and reflection from the effective potential well of a three atomic system (see Figs. 3 (b) and 3 (a)), arising in the scattering process were revealed, that are generated by metastable states of the composite system (cluster + barrier or well) with complex energy eigenvalues below the dissociation threshold of dimer, corresponding to the shape and Feshbach resonances, respectively.
The elaborated method, algorithms and programs [8][9][10] for solving the three-atomic scattering problem as well as diatomic molecule tunneling through potential barrier can be applied to the further analysis of quantum transparency and reflection effects [13], quantum diffusion [22] and the resonance scattering in triatomic systems using modern theoretical and experimental results [23,24].

Figure 1 .
Figure 1.(a) The profiles of 2D potential functions of Be 3 in Jacobi coordinates (1) in the center-of-mass plane and the relative arrangement of particles in accordance with the region of the center-of-mass plane.The numbers of sectors are given in boxes.(b) The profiles of 2D potential functions of Be 2 with barrier and the relative arrangement of particles and barrier.Here coordinates and potential functions are given in Å and Å−2 , respectively

Figure 2 .
Figure 2. The density plots of the eigenfunctions Ψ ,  (, ) displayed in sector 1 of (, )-plane (in Å) of the gerade () and ungerade () bound states with energies  ,  of the Be 3 trimer presented in Table 1.The negative, positive and near-zero values of the eigenfunctions are displayed by black, white and gray, respectively

Figure 5 .
Figure 5.The components   and isolines of the absolute values |Ψ(, )| of the solution Ψ(, ) for the trimer Be 3 in 1g metastable state with the real part of energy eigenvalues Re  = −188.94Å−2 , localized near the minimal of the trimer potential

Figure 6 .
Figure 6.The total probability |T| 2  () = ∑︀   =1  *    (lines) of penetration of the dimer for the initial states  through the repulsive Gaussian potential barriers versus the total energy  = Re  counted from the main threshold  = 0.The values of the threshold energies  =   ,  = 1, . . ., 5 corresponding to the energies of ground and excited initial states are shown by arrows.The real Re  and imaginary (−1) Im  part (with negative sign) of the metastable states energy (circles) the open channels   = max (), where  = diag{   }     =1 , is a diagonal matrix.In open channels it is defined as the matrix transforming the amplitudes of the incoming waves  −