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

**Authors:**Gusev AA^{1}-
**Affiliations:**- Laboratory of Information Technologies Joint Institute for Nuclear Research

**Issue:**Vol 26, No 3 (2018)**Pages:**226-243**Section:**Mathematical Modeling**URL:**https://journals.rudn.ru/miph/article/view/18988**DOI:**https://doi.org/10.22363/2312-9735-2018-26-3-226-243

Cite item

## Full Text

## Abstract

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.

## Full Text

1. 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-3]. The FEM schemes and algorithms for solving the d-dimensional BVP are presented in [4-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-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. 2. Setting of the Problem Consider a model of three identical particles (trimer) on a straight line with the masses �� = � and the coordinates �¯� ∈ R1, � = 1, 2, 3, of one-dimensional Euclidian space R1 coupled via the pair short-range potential �˜ (|�¯� - �¯� |), �, � = 1, 2, 3. Performing the change of variables corresponding to the cyclic permutation (�, �, �) = (1, 2, 3) of the permutation group �3 [12]: �¯� + �¯� - 2�¯� √2 �� ≡ �(��) = �¯� - �¯�, �� ≡ �(��)� = √3 , �0 = √3 (�¯1 + �¯2 + �¯3), � we get three pairs (︀�(��), �(��)� )︀ of the scaled Jacobi variables. In the center-of-mass reference frame of the configuration space x¯ = (�¯1, �¯2, �¯3) ∈ �3, i.e., on the hyperplane Ω = {x¯ = (�¯1, �¯2, �¯3) |�¯1 + �¯2 + �¯3 = 0} ⊂ �3 these Jacobi maps (�, �)� ∈ Ω�� ∼ �2 are connected by the relations: (︂�� )︂ �� ≡ (︂ �(��) )︂ �(��)� = T�� (��) (︂��)︂ �� , (1) T�� (��) = (︂ cos �� sin �� )︂ , - sin �� cos �� where the angles �� = �� = �� = -�� = -�� = -�� = (2�)/3, i.e. cos �� = cos �� = cos �� = -1/2, sin �� = sin �� = sin �� = √3/2, are such that �� + �� + �� = 2�, �� + �� + �� = -2�. With the variable parameter in T�� () they and simply change the sign of pairs of Jacobi coordinates in Eq. (1), (︂ -�(��) )︂ = T ( (︂ �(��) )︂ �) . (2) -�(��)� �� �� ± �(��)� 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¨odinger equation for the wave function Ψ(�, �) corresponding to the total energy � of the triatomic system in Jacobi coordinates (�, �)� = (�� , �� )� ∈ Ω�� : (︂ ∂2 ∂2 � ˜ ˜ )︂ -∂�2 - ∂�2 + n2 (� (�, �) - �) Ψ(�, �) = 0. (3) 228 RUDN Journal of MIPh. Vol. 26, No 3, 2018. Pp. 226-243 ✭❛✮ ✭❜✮ Figure 1. (a) The profiles of 2D potential functions of Be3 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 Be2 with barrier and the relative arrangement of particles and barrier. Here coordinates and potential functions are given in ˚A and ˚A-2, respectively In the case of a diatomic molecule (dimer) with identical nuclei coupled via the pair potential, �˜ (|�¯1 - �¯2|), moving in the external potential field �˜ �(|�¯� - �¯3|), � = 2, 1, of the third atom having the infinite mass, the same equation (3) is valid for the variables � ≡ �3 = �¯2 - �¯1, � ≡ �3 = �¯1 + �¯2, with respect to the origin of the coordinate frame placed at the infinite-mass atom, �¯3 = 0. 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 (a)), �˜ (�, �) = �˜ (|�|) + �˜ (︃⃒ - ⃒� ⃒ ⃒ √ ⃒)︃ 3� ⃒ ⃒ ⃒ + �˜ (︃⃒ ⃒� + ⃒ ⃒ √ ⃒)︃ 3� ⃒ ⃒ ⃒ , (4) ⃒ 2 ⃒ ⃒ 2 ⃒ 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 (b)) (︂⃒� - � ⃒)︂ (︂⃒� + � ⃒)︂ ⃒ ⃒ �˜ (�, �) = �˜ (|�|) + �˜ � ⃒ ⃒ ⃒ ⃒ ⃒ + �˜ � ⃒ 2 ⃒ ⃒ ⃒ , (5) 2 ⃒ is an even function with respect to the straight line � = 0 (i.e., �¯1 = �¯2), 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 (︂ d2 � ˜ )︂ -d�2 + n2 (� (|�|) - ˜) �(�) = 0. (6) Gusev A. A. Finite Element Method of High-Order Accuracy for Solving . . . 229 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]. The approach proposed below is illustrated by the example of Be2 with the reduced mass �/2 = 4.506 Da of the nuclei [13] and the molecular interaction approximated by the Morse potential � � (�) = �˜ (�), n2 �˜ (�) = �˜ {exp[-2(� - �ˆ�� )�] - 2 exp[-(� - �ˆ�� )�]}. (7) Here � = 2.96812 ˚A-1 is the potential well width, �ˆ�� = 2.47 ˚A is the average distance between the nuclei, and �˜ = 1280 K, � = (�/n2)�˜ = 236.51 ˚A-2 (1 K=0.184766 ˚A-2, 1 ˚A-2=5.412262 K) is the potential well depth. This potential supports only five bound states corresponding to the covalent bounding of the Be2 dimer [14] having the energies � = (︀�/n2)︀ ˜�, � = 1, . . . , �0 = 5. The parameter values are determined from the condition (˜2 - ˜1) / (2�n�) = 277.124 cm-1, 1 K/(2�n c)=0.69503476 cm-1. To solve the discrete spectrum problem we applied the seventh-order FEM using the Hermitian interpolation polynomials with double nodes [8]: � = {-193.066, -119.392, -63.338, -24.904, -4.089}˚A-2, � = 1, . . . , �0 = 5. As an exam- � ple, following [13], we use below the Gaussian barrier potentials � �(��) = � exp (︀-�2/�)︀ with � = 236.51 ˚A-2 and � = 0.0523 ˚A2. The values of parameters of the repulsive Gaussian barrier potential were estimated following the experimental observation of the quantum diffusion of hydrogen atoms on the copper surface [15]. 3. Boundary-Value Problems Using the change of variables � = � sin , � = � cos , we rewrite Eq. (3) in polar coordinates (�, ), Ω�, = (� ∈ (0, ∞), ∈ (0, 2�)) (︂ 1 ∂ ∂ 1 )︂ d2 2 -� ∂��∂� + �2 Λ(, �) - � Ψ(, �) = 0, Λ(, �) = -d2 + � � (, �), (8) where for a trimer with pair potentials � (, �) = � (�| sin |) + � (�| sin( - 2�/3)|) + � (�| sin( - 4�/3)|), (9) and for a dimer with pair potential in the external field of barrier potential: � (, �) = � (�| sin |) + � �(�| sin( - �/4)|) + � �(�| sin( + �/4)|). (10) where � �(�| sin( ± �/4)|) = � exp (︀-�2 sin2( ± �/4)/�)︀. The solution of Eq. (8) is sought in the form of the Kantorovich expansion [16] Ψ�� (, �) = �max ∑︁ �=1 �� (; �)���� (�). (11) Here ���� (�) are unknown matrix functions, � = 1, . . . , �max = 2�. The angular basis functions �� (; �) ∈ �� ∼ �2(Ω) in the interval Ω = ∈ [0, 2�), which is divided into �max 230 RUDN Journal of MIPh. Vol. 26, No 3, 2018. Pp. 226-243 �=1 subintervals Ω� = ∈ (2�(� - 1)/�max, 2��/�max): Ω�; = ⋃︀�max Ω�, are determined at each value of the parameter � ∈ (0, +∞) as the eigenfunctions corresponding to the real discrete eigenvalues 1(�) < 2(�) < . . . < � (�) < . . . of the Sturm-Liouville problem for the equation (︀Λ(, �) - � (�))︀�� (; �) = 0. (12) � The functions �(�)(; �) ≡ �� (; �) have the parity (-1)�, � = 0, 1 with respect to the inversion of Jacobi coordinates (2), i.e., �(�)(; �) = (-1)��(�)( ± �; �); here and below � � � = 0. These functions satisfy the orthogonality and completeness conditions �max ⟨��|�� ⟩ = ∑︁ 2��/�max ∫︁ d��(; �)�� (; �) = ��� , (13) �=12�(�-1)/�max ∑︁ ∑︁ |��⟩⟨��| = � ��(; �)��(0; �) = �( - 0). (14) � For the three problems under consideration the potential function � (, �) depending on the parameter � can be defined as follows. Task 1. The case of one pair potential in the intervals ∈ (0, 2�) (� = �/6, �/4 or �/2) � (, �) = � (� sin ). Task 2. The case of three pair potentials, Eq. (9), in the interval ∈ (0, 2� = �/3); the potential has the symmetry of the �3� dihedral point group [12]. Task 3. The case of one pair potential and two penetrable or almost impenetrable barriers, Eq. (10), in the interval ∈ (0, � = �/2) or in the intervals ∈ (0, � = �/4-) and ∈ (� = �/4 - , �/2), 0 < ≪ �/4. The solutions symmetric with respect to the permutation of two particles satisfy the Neumann boundary condition at = 0 and = 2�, while the antisymmetric ones satisfy the Dirichlet boundary condition at these points, i.e., at the aforementioned points � = ��� = 0. If the pair potential possesses a high peak in the vicinity of the pair collision point, then the solution of the problem (8) will be considered in the half-plane Ω,� = (� ∈ (�min, ∞), ∈ [min(�), 2� - min(�)]) with the Neumann or Dirichlet boundary condition. The potential function � (, �) of the boundary-value problem (12) is symmetric under the reflection �ˆ�: → (4� - 2)� - with respect to the lines = (2� - 1)� in each sector of the cycle, numbered by � = 1, . . . , �max: �ˆ�� (, �) = � (, �). Therefore, the set of eigenfunctions �� (; �) is separated into two subsets, namely, the even and odd ��=±1(; �) ones: �ˆ��� (; �) = �� ((4� - 2)� - ; �) = ±��=±1(; �). � � � � This fact allows separate calculation of gerade �� (; �) = �� ((4� - 2)� - ; �) or � � ungerade ��(; �) = -�� ((4� - 2)� - ; �) eigenfunctions in the reduced interval � � � ∈ [min(�), �], subjecting them to Neumann or Dirichlet boundary condition at the boundary point = � of the interval, respectively. Below the parametric angular basis functions ��=±1(; �) with the numbers � = 1, . . . , �0 are referred to as cluster states with � (�) < 0, and those with � � �0 + 1 as pseudostates with � (�) > 0, corresponding to the discrete and continuous spectrum of BVP for Eq. (6) at large values of the parameter �, respectively. To reveal the above structurization property, we introduce the linear Gusev A. A. Finite Element Method of High-Order Accuracy for Solving . . . 231 combinations of these functions ��,� (; �) for the trimer and �←,→(; �) for the dimer: ��,� � �=-1 � �=1 √ � (; �) = (︀±�� (; �) + �� (; �))︀ / 2, �←,→ �=-1 �=1 √ (15) � (; �) = (︀±�� (; �) + �� (; �))︀ / 2. The action of the parity operator �ˆ� on these functions �ˆ��� (; �) = �� (; �), �ˆ��� (; �) = �� (; �), � � � � (16) �ˆ��←(; �) = �→(; �), �ˆ��→(; �) = �←(; �), � � � � consists in the permutation of states |�⟩ ↔ |�⟩ or | ←⟩ ↔ | →⟩ with respect to the lines = (2� - 1)�. Indeed, the parametric cluster functions at � = 1, . . . , �0 and at large � have maximal in the vicinity of = (2� - 2)� and = 2��, respectively, that correspond to the eigenfunctions of cluster states of the BVP for Eq. (6). In the particular � case of sector 1, the dimer functions �←,→(; �) have maximal in the vicinity of = 0 � and = �, i.e. for � > 0, � =←, or for � < 0, � =→ at large �, respectively while the trimer functions ��,� (; �) have maximal in the vicinity of = 0 and = �/3, i.e. for �12/�(12)3 > 0, � = �, or for �31/�(31)2 < 0 � = � at large �, respectively. Using the above pairs of the basis functions (15), we rewrite the expansion (11) in the �-representation � ∑︁ ∑︁ � � Ψ�� (, �) = �=1 �� (; �)���� (�), (17) � �� where � = �, � or � = �, � for the unknown functions �� � �� (�) and �� � (�) in the (��)- �� representation, related to the functions �� � ��� (�) and �� (�) in the (��)-representation as ���� (�) = �� (︃�� )︃ ��� (�) ��� (�) ��� (︂�� (�))︂ �� = � ��� (�) 1 , � = √ 2 (︂ 1 1)︂ -1 1 . (18) 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 [︂ 1 d d �(�) ]︂ �max ∑︁ -� d��d� + �2 - � ���� (�) + �=1 ��� (�)���� (�) = 0, (19) 1 d d ��� (�) = ���(�) + � d�����(�) + ���(�) d�. (20) 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: ��� (�) = -⟨��|∂��� ⟩, ��� (�) = ⟨∂���|∂��� ⟩. (21) �� For Task 3 the effective potentials �ˆ �� (�) = ��� (�) + � � (�) are sums of ��� (�), calculated using the potential curves and the parametric basis functions of Task 1, and 232 RUDN Journal of MIPh. Vol. 26, No 3, 2018. Pp. 226-243 �� the matrix elements � � (�) of the barrier potentials � � �� (�) = ⟨��|� �(�| sin( - �/4)|) + � �(�| sin( + �/4)|)|�� ⟩. (22) 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 Be2 dimer and Be3 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 Be2 dimer and � = �/6 for Be3 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. (︂����� (�) ����� (�))︂ (︂�����(�) 0 )︂ 1. (23) ��� (�) = ����� (�) � ���� (�) = � 0 ����� (�) �- 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. 4. 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] �(2) (�� �) and �(1) (�� �) below the 1/2 � 1/2 � dissociation threshold at � < 0 in the form of incoming and outgoing waves � �� �� ∑︁ [︁ - - + + ]︁ � (�) = ���′ (��� �)��′ + ���′ (��� �)��′ , � � � � �′ � =1 � at ��� = √︀� - �� > 0 in the open channels �� = max � ::: �. The scattering matrix �� (�) or �� } (�), where � = diag{�� �� � � , is a diagonal matrix. In open channels � ���′ � ���′ � �� =1 �′ it is defined as the matrix transforming the amplitudes of the incoming waves �- into � �′ those of the outgoing waves �+ � [18] �� �+ ∑︁ � - � �′ �� ′ = � �′ � =1 �� �� (�)� . (24) ��� The components of the radial asymptotic solutions � �� (�) of the scattering problem in the open channels �� = 1, . . . , �� have the form � �� ��� (�) = �� ∑︁ [︁ ��′ � �- (�� �)���′ + �+ � � �′ � (�� �)�� (�)]︁ , (25) � ′ �� =1 � � ��′ � Gusev A. A. Finite Element Method of High-Order Accuracy for Solving . . . 233 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 � �� + + ��� (�) = ��′ ���′ (��� �). (26) � � � � � 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) at large � = �max �� d d �� ���� (�) = ���� (�), d� d� ���� (�) = ���� (�). (27) 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: d� (�) ⃒ ⃒ � (�min) = 0, ⃒ d� ⃒�=�max = �(�max)� (�max), (28) � }�� =1 where �(�max) is a � × � symmetric matrix function of �, � (�) = {�� (�) � = � � {{���� (�)}�=1}�� =1 is the required � × � numerical matrix solution. ��� ±1 �′ These matrices and the �� × �� matrices S = {��,� (�) �� ′ }��,�� =1 sought for in the open channels �� = max � ::: � are calculated directly from (27) using the program � � KANTBP 4M [8]. 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] � {︂ ��, Re �� � � }︂ √︁ � √︁ � ����� (�max) = � ��, Re �� < � , �� = �� - �� , �� = �� - ��, (29) 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 ⎡ �max ⎤ ∫︁ (F�|F�′ ) = (�� + ��′ ) ⎣ �min � F� (�)F�′ (�)d� - ���′ ⎦ + ���′ = 0, (30) � ���′ = -F� (�max)F�′ (�max). 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 ⟨F�|F�′ ⟩ = �max ∫︁ � F� (�)F�′ (�)d� = ���′ . (31) �min 234 RUDN Journal of MIPh. Vol. 26, No 3, 2018. Pp. 226-243 + 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 (︂S�� S�� )︂ S = S�� S�� , S S = SS = I, (32) † † �� 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, |S�� |2 = I - |S�� |2 �� and S�� = S� , 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] (︂R← T→)︂ ← → S = T R , S†S = SS† = I, (33) where I is the unit matrix with the dimension 2�� ×2�� consisting of the amplitudes of the reflected and transmitted waves R� = R� (�) and T� = T� (�), where � =←, → indicates the direction of the incident wave propagation with respect to the �-axis, i.e., � =← and � =→ for � > 0 and � < 0, respectively, and the �� × �� matrices R� = R� (�) and � → T� = T� (�) are expressed as R← = R→ = (S+1 + S-1)/2, T← = T = (S+1 - S-1)/2. For the scattering of the dimer (��) on potential barriers similar relations determine � | � the reflection R← = S�� and R→ = S�� , and transmission T← = S�� and T→ = S�� amplitudes. For the reflection coefficient |R� |2 = R† R� and the transmission coefficient T� |2 = T† T� the conventional relation below breakup threshold at � < 0 following from (32) and constant Wronskian, |T� |2 = I-|R� |2 is valid, where I is the unit �� ×�� matrix. 5. 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 set of even (g) and odd (u) bound states of the trimer Be3 (Task 2 ) were calculated on the grid Ω� = {�min = 4.24(1)4.33(10)6.13(1)6.33(23)�max = 11.39}, where in parentheses the number of fifth-order Hermitian elements [20] is indicated, for the number of equations � = 15 in the system (19). The comparison of sets of total and binding energies of g and u bound states of the trimer Be3 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˚A-2. Gusev A. A. Finite Element Method of High-Order Accuracy for Solving . . . 235 � Table 1 The total ��(*) and binding energies -��(*)=��(*)-�trsh of g and u bound states of relative to the exact threshold �trsh= - 193.06 (in ˚A-2): � (�) and the trimer ��3 � -�� � (�) calculated by KM in polar coordinates using � =12 basis functions (g or � u); ��(� ) and -��(� ) calculated by the 2D FEM using fifth order interpolation Lagrange polynomials on triangular finite elements in sector 1 of Jacobi coordinates � 1g 2g 3u 4g 5u 6g -��(�) 389.08 335.43 300.58 287.02 260.47 245.84 -��(�) � 196.02 142.37 107.52 93.96 67.41 52.78 -��(� ) 389.09 335.45 300.60 287.05 260.50 245.88 -��(� ) � 196.03 142.39 107.54 93.99 67.44 52.82 � 7u 8g 9u 10g 11u 12g -��(�) 227.66 225.39 215.37 204.85 198.21 193.86 -��(�) � 34.60 32.33 22.31 11.79 5.15 0.80 -��(� ) 227.70 225.42 215.41 204.89 198.24 193.91 -��(� ) � 34.64 32.36 22.35 11.83 5.18 0.85 � Figure 2. The density plots of the eigenfunctions Ψ�,�(, �) displayed in sector 1 of ��,� (�, �)-plane (in ˚A) of the gerade (�) and ungerade (�) bound states with energies � of the Be3 trimer presented in Table 1. The negative, positive and near-zero values of the eigenfunctions are displayed by black, white and gray, respectively 236 RUDN Journal of MIPh. Vol. 26, No 3, 2018. Pp. 226-243 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 Be3 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 |���|2 is achieved, the number � of the threshold �, the real and imaginary part of the complex energy eigenvalues �� � = Re �� � + Im �� in ˚A-2 of the even � and odd � metastable states of Be3 numbered by the index � calculated with � = 15 equations (19) � � |���|2 �� Re �� � Im �� � �type -193.066 1 thr -189.676 1 1 · 10-7 -188.94 -4 · 10-2 1g -164.654 1 4 · 10-6 -164.72 -1 · 10-2 1u -156.882 1 3 · 10-6 -157.04 -2 · 10-2 2g -140.545 1 1 · 10-4 -140.57 -5 · 10-3 2u -132.485 1 1 · 10-6 -132.47 -4 · 10-3 3g -124.256 1 1 · 10-3 -124.16 -6 · 10-3 3u -120.638 1 1 · 10-6 -120.75 -8 · 10-2 4g -119.392 2 thr -113.248 -89.319 2 2 0.10 9 · 10-4 -113.24 -89.16 -2 · 10-2 -3 · 10-2 5g 6g -77.271 2 0.77 -76.51 -4 · 10-6 4u -70.309 2 0.35 -70.30 -2 · 10-3 7g -63.385 2 0.41 -65.14 -3 · 10-4 5u -63.338 3 thr -42.858 3 0.06 -42.87 -6 · 10-3 8g -29.396 3 0.13 -29.19 -4 · 10-2 9g -24.899 3 0.19 -25.82 -1 · 10-3 6u -24.904 4 thr -6.799 4 0.40 -7.12 -1 · 10-3 10g -4.089 5 thr 0 thr These metastable states are responsible for resonance energies, corresponding to the minimal probability of inelastic scattering of the dimer by the atom, i.e., to the resonance quantum reflection from the potential well (Feshbach resonances, see Figs. 1 (a) and 3 (a). As an example, in Fig. 4 we display the eigenfunctions of the scattering problem for gerade and ungerade states corresponding to the minimum of the transmission coefficient |��� |2 = |��� |2 = ·10-7 at � = -189.676, as well as the metastable state 1g from Table 2. √ The isolines of the absolute values |Ψ�,� (�, �)| of the linear combinations Ψ�,� (�, �) = (Ψ� (�, �) ± Ψ�(�, �))/ 2 demonstrates the effect of resonance reflection from the effective Gusev A. A. Finite Element Method of High-Order Accuracy for Solving . . . 237 ✭❛✮ ✭❜✮ �� Figure 3. The effective diagonal potentials ��� (�) = � (�)�-2 + ��� (�) for the Be3 trimer (a) and effective diagonal potentials ��� (�) = � (�)�-2 + ��� (�) + � � (�) for the tunneling problem of the dimer Be2 through the Gaussian barrier (b) Figure 4. Upper panel: the isolines of the absolute values |Ψ�,�(, �)| of corresponding gerade (left-hand panel) and ungerade (right-hand panel) solutions in sector 1 of (�, �)-plane (in ˚A) for the scattering of Be atom with the energy -� = 189.676 ˚A-2 (relative to the three-body threshold) on the dimer Be2, corresponding to the 1g metastable state from Table 2. Lower panel: isolines of the absolute values |Ψ�,� (, �)| of the linear combinations Ψ�,� (, �) = (±Ψ�(, �) + Ψ� (, �))/√2 238 RUDN Journal of MIPh. Vol. 26, No 3, 2018. Pp. 226-243 potential well. It can be seen from the figures that the shape of the wave functions of the gerade scattering states (Fig. 4 (a)) and metastable states (Fig. 5) are similar and they are localized in the vicinity of the potential well (Fig. 3 (a)). At the same time, for the same energy value � = -189.6 ˚A-2, the wave function of ungerade scattering states (Fig. 4 b) is a typical nonresonant wave function. Figure 5. The components �� and isolines of the absolute values |Ψ(, �)| of the solution Ψ(, �) for the trimer Be3 in 1g metastable state with the real part of energy eigenvalues Re � = -188.94 ˚A-2, localized near the minimal of the trimer potential 6. Metastable and Scattering States of the Dimer Tunneling The metastable states of the dimer Be2 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 Be2 tunneling through the Gaussian barriers, are presented in Fig. 6. 11 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 (�) 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 Gusev A. A. Finite Element Method of High-Order Accuracy for Solving . . . 239 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. Figure 6. The total probability |T|2 (�) = ∑︀�� � * ��� (lines) of penetration of the �� �=1 �� 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) 7. Conclusion The model for three atomic beryllium system in a straight line was formulated as a 2D boundary-value problem for the Schr¨odinger 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-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].×

## About the authors

### A A Gusev

Laboratory of Information Technologies Joint Institute for Nuclear Research
**Author for correspondence.**

Email: gooseff@jinr.ru

Candidate of Physical and Mathematical Sciences, senior researcher of Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow region, Russia, 141980

## References

- P. Ciarlet, The Finite Element Method for Elliptic Problems, North-holland Publ. Comp, Amsterdam, 1978. doi: 10.1137/1.9780898719208.
- V. G. Korneev, Schemes of the Finite Element Method for High Orders of Accuracy, Leningrad State University, Leningrad, 1977, in Russian.
- A. A. Gusev, V. P. Gerdt, O. Chuluunbaatar, G. Chuluunbaatar, S. I. Vinitsky, V. L. Derbov, A. Go´´zd´z, Symbolic-Numerical Algorithm for Generating Interpolation Multivariate Hermite Polynomials of High-Accuracy Finite Element Method, Lecture Notes in Computer Science 10490 (2017) 134-150. doi: 10.1007/978-3-319-66320-3 11.
- A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar, G. Chuluunbaatar, V. P. Gerdt, V. L. Derbov, A. Go´´zd´z, P. M. Krassovitskiy, High-Accuracy Finite Element Method: Benchmark Calculations, European Physics Journal - Web of Conferences 173 (2018) 03009.
- A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar, G. Chuluunbaatar, V. P. Gerdt, V. L. Derbov, A. Go´´zd´z, P. M. Krassovitskiy, Interpolation Hermite Polynomials For Finite Element Method, European Physics Journal - Web of Conferences 173 (2018) 03010.
- A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar, G. Chuluunbaatar, V. P. Gerdt, V. L. Derbov, A. Go´´zd´z, P. M. Krassovitskiy, Algorithm for calculating interpolation hermite polynomials for high-accuracy finite element method, in: Computer Algebra: International Conference Materials, Plekhanov Russian University of Economics, 2017, pp. 89-95.
- A. A. Gusev, V. P. Gerdt, O. Chuluunbaatar, G. Chuluunbaatar, S. I. Vinitsky, V. L. Derbov, A. G´o´zd´z, Symbolic-Numerical Algorithms for Solving the Parametric SelfAdjoint 2D Elliptic Boundary-Value Problem Using High-Accuracy Finite Element Method, Lecture Notes in Computer Science 10490 (2017) 151-166. doi: 10.1007/9783-319-66320-3 12.
- A. A. Gusev, L. L. Hai, O. Chuluunbaatar, S. I. Vinitsky, Program KANTBP 4M for Solving Boundary-Value Problems for Systems of Ordinary Differential Equations of the Second Order. URL http://wwwinfo.jinr.ru/programs/jinrlib/kantbp4m
- A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, A. G. Abrashkevich, KANTBP 3.0: New version of a Program for Computing Energy Levels, Reflection and Transmission Matrices, and Corresponding Wave Functions in the Coupled-Channel Adiabatic Approach, Computer Physics Communications 185 (2014) 3341-3343. doi: 10.1016/j.cpc.2014.08.002.
- O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich, ODPEVP: A Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Parametric Self-Adjoined Sturm- Liouville Problem, Computer Physics Communications 180 (2009) 1358-1375. doi: 10.1016/j.cpc.2009.04.017.
- A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, V. L. Derbov, Algorithms for Solving the Boundary-Value Problems for Atomic Trimers in Collinear Configuration using the Kantorovich Method, Bulletin of Peoples’ Friendship University of Russia. Series: Mathematics. Information Sciences. Physics (4) (2016) 56-76.
- J. F. Cornwell, Group Theory in Physics, Academic Press, New York, 1984.
- P. M. Krassovitskiy, F. M. Pen’kov, Contribution of Resonance Tunneling of Molecule to Physical Observables, Journal of Physics B: Atomic, Molecular and Optical Physics 47 (22) (2014) 225210. doi: 10.1088/0953-4075/47/22/225210.
- J. Wang, G. Wang, J. Zhao, Density Functional Study of Beryllium Clusters, with Gradient Correction, Journal of Physics: Condensed Matter 13 (33) (2001) L753-L758. doi: 10.1088/0953-8984/13/33/101.
- L. J. Lauhon, W. Ho, Direct Observation of the Quantum Tunneling of Single Hydrogen Atoms with a Scanning Tunneling Microscope, Physical Review Letters 85 (2000) 4566-4569. doi: 10.1103/PhysRevLett.85.4566.
- L. V. Kantorovich, V. I. Krylov, Approximate Methods of Higher Analysis, Wiley, New York, 1964.
- M. Abramovits, I. A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972.
- R. G. Newton, Analytic Properties of Radial Wave Functions, Journal of Mathematical Physics 1 (1960) 319-348. doi: 10.1063/1.1703680.
- A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar, V. L. Derbov, A. G´o´zd´z, P. M. Krassovitskiy, Metastable States of a Composite System Tunneling through Repulsive Barriers, Theoretical and Mathematical Physics 186 (2016) 21-40. doi: 10.1134/S0040577916010037.
- A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, V. L. Derbov, A. G´o´zd´z, L. L. Hai, V. A. Rostovtsev, Symbolic-Numerical Solution of Boundary-Value Problems with Self-Adjoint Second-Order Differential Equation using the Finite Element Method with Interpolation Hermite Polynomials, Lecture Notes in Computer Science 8660 (2014) 138-154. doi: 10.1007/978-3-319-10515-4 11.
- I. V. Puzynin, T. L. Boyadjiev, S. I. Vinitsky, E. V. Zemlyanaya, T. P. Puzynina, O. Chuluunbaatar, Methods of Computational Physics for Investigation of Models of Complex Physical Systems, Physics of Particles and Nuclei 38 (2007) 70-116. doi: 10.1134/S1063779607010030.
- E. Pijper, A. Fasolino, Quantum Surface Diffusion of Vibrationally Excited Molecular Dimers, Journal of Chemical Physics 126 (2007) 014708. doi: 10.1063/1.2424699.
- A. V. Mitin, Unusual Chemical Bonding in the Beryllium Dimer and its Twelve Vibrational Levels, Chemical Physics Letters 682 (2017) 30-33. doi: 10.1016/j.cplett.2017.05.071.
- J. M. Merritt, V. E. Bondybey, M. C. Heaven, Beryllium Dimer-Caught in the Act of Bonding, Science 324 (2009) 1548-1551. doi: 10.1126/science.1174326.