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

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