Algorithms for Solving the Boundary-Value Problems for Atomic Trimers in Collinear Configuration using the Kantorovich Method

Full Text

1. Introduction At the present time, the processes of resonance scattering of diatomic molecules by atoms via three-particle clustering and metastable states, as well as the processes of dissociation of the molecule induced by collisions with atoms are a subject of intense theoretical and experimental studies [1-3]. To analyse such processes they use triatomic model systems with the atoms bound by pair realistic molecular and van-der-Waals potentials, which possess bound and metastable states in the vicinity of the dissociation threshold of the diatomic molecule. The study of such models stimulates the development of new methods and symbolic-numerical algorithms for solving multidimensional boundary value problems with non-separable variables and the construction of asymptotic states of the scattering problem in a system of three atoms, based on the Kantorovich method and the finite element method (FEM), implemented in the problem-oriented software packages [4-9]. The aim of this paper is to formulate the 2D BVPs in the Jacobi and polar coordinates for the model of an atomic trimer in the collinear configuration and to develop the algorithms for calculating the asymptotic parametric angular functions and the effective potentials arising in the Kantorovich method and needed for constructing the asymptotic states of the triatomic scattering problem. The developed algorithms will be tested by the example of beryllium atoms trimer model. The paper is organised as follows. In Section 2 we set the 2D BVP. In Section 3 the 2D BVP is reduced to the 1D BVP for a set of second-order ODEs using the Kantorovich Received 16th September, 2016. The authors thank Profs. V.P. Gerdt, F.M. Penkov and Drs. P.M. Krassovitskiy, L.L. Hai for collaboration. The work was supported by the Russian Foundation for Basic Research (grants No. 14-01-00420, 17- 01-00298). The reported study was funded within the Agreement N 02.a03.21.0008 dated 24.04.2016 between the Ministry of Education and Science of the Russian Federation and RUDN University. method. As an example, the eigenvalues and the hyperradial eigenfunction components are calculated by means of the Kantorovich method and FEM for the model of beryllium trimer in the collinear configuration. In Section 4 we present the algorithms for calculating the asymptotes of the parametric basis functions in polar coordinates at large values of the parameter (hyperradial variable) and the effective potentials of the system of ODEs. In Conclusion the results and perspectives are discussed. 2. Setting of the problem Consider a 2D model of three identical particles with the mass and the coordinates ∈ | - | ∈ | - | R1, = 1, 2, 3, coupled via the pair potential ˜ ( ) , = 1, 2, 3. Performing the change of variables at cyclic permutation (, , ) = (1, 2, 3): + - 2 √2 ≡ () = - , ≡ () = √3 , 0 = √3 (1 + 2 + 3), we arrive at the Schro¨dinger equation for the wave function in the center-of-mass system { ∈ R1|1 + 2 + 3 = 0} - ∂2 - ∂2 + 2 ( (, ) - ) - ∂2 - ∂2 + 2 ( (, ) - ) Ψ(, ) = 0. (1) Ψ(, ) = 0. (1) (︂ ∂2 ∂2 ˜ ˜ )︂ In the case of a diatomic molecule with identical nuclei coupled via the pair potential | - | | - | | - | | - | ˜ ( 1 2 ) and moving in the external potential field ˜ ( 3 ), = 2, 1 of the third atom having the infinite mass, the same equation (1) is valid for the variables = 1 - 2, = 1 + 2, the origin of the coordinate frame being placed on the infinite-mass atom, 3 = 0. Here the potential function for a trimer with the pair potentials (below this case is referred to as Task 2), √ ⃒)︁ ˜ (︁⃒ √ ⃒)︁ √ ⃒)︁ ˜ (︁⃒ √ ⃒)︁ 2 2 2 2 ˜ (, )=˜ (||)+˜ (︁⃒⃒⃒ - 3 ⃒⃒ + ⃒⃒ + 3 ⃒⃒ , (2) or the potential function for a dimer in the field of barrier potentials (below this case is referred to as Task 3) 2 2 2 2 ˜ (, )=˜ (||)+˜ (︁⃒⃒⃒ - ⃒⃒⃒)︁+˜ (︁⃒⃒⃒ + ⃒⃒⃒)︁, (3) is symmetric 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. Using 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 maximum in the vicinity of 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 setting the Neumann or Dirichlet boundary condition at min gives only a minor contribution to the solution. The equation, describing the molecular subsystem (dimer), has the form -d2 + 2 ( () - ˜) -d2 + 2 ( () - ˜) () = 0. (4) () = 0. (4) (︂ d2 ˜ )︂ -| | -| | We assume that the dimer has the discrete spectrum, consisting of a finite number 0;;? 1 of bound states with the eigenfunctions (), = 1, 0 and eigenvalues ˜ = ˜ < 0, and the continuous spectrum of eigenvalues ˜ > 0 with the corresponding eigenfunctions ˜(). As a rule, the solutions of the discrete and continuum spectra of the BVP for Eq. (4) can be found only numerically, except simplified models having exact solutions and used for computer modelling of bimolecular chemical reactions. In certain cases the eigenfunctions of the continuous spectrum are approximated by the eigenfunctions of pseudostates of discrete spectrum ˜ > 0, = 1 + 0, . . . calculated in sufficiently large but finite interval [10]. The proposed algorithm is illustrated by the example of the molecular interaction approximated by the Morse potential of Be2 with the reduced mass /2 = 4.506 Da of the nuclei [10, 11] ˜ ˜ () = 2 (), ˜ () = {exp[-2( - ˆ )��] - 2 exp[-( - ˆ )]}. (5) Here = 2.96812 ˚A-1 is the potential well width, ˆ = 2.47 ˚A is the average distance between the nuclei, and = 1280 K (1 K=0.184766 ˚A-2,1 ˚A-2=5.412262 K) is the potential well depth. This potential supports five bound states [12] having the - - energies = (/ 2)˜, = 1, . . . , 0 = 5 presented in Table 1. The parameter values are determined from the condition (˜2 ˜1)/(2 ) = 277.124cm-1, 1K/(2 c)=0.69503476 cm-1. Table 1 The discrete spectrum energies of dimer Be2 and the binding energy - - - - = ( ) of gerade () and ungerade () states of trimer Be3 counted of - - - - = ˜1 = 193.06˚A-2 = 1044K dimer energy Be2 calculated in grid Ω = 4.1(20)7(10)10 at max = 10 -˜ , , -˜1=1044.879 649K -˜2=646.157 093K -˜3=342.791 979K -˜4=134.784305K -˜5=22.134 073K = 196.02˚A-2 = 1060.86K 1, = 142.37˚A-2 = 770.51K 2, = 93.95˚A-2 = 508.50K 3, = 52.77˚A-2 = 285.63K 4, = 32.32˚A-2 = 174.95K 5, = 22.31˚A-2 = 120.75K = 107.52˚A-2 = 581.90K 1, = 67.41˚A-2 = 364.84K 2, = 34.60˚A-2 = 187.28K 3, = 11.79˚A-2 = 63.83K 4, = 0.8˚A-2 = 4.4K 5, = 5.14˚A-2 = 27.87K To solve the discrete spectrum problem we applied the FEM of the seventh or- der using the Hermitian interpolation polynomials with double nodes [9]. The grid { } { } 0, . . . , , . . . , was used to calculate the values of both the function and its deriva- tives. 3. Reduction of the BVP using the Kantorovich method Using the change of variables = sin , = cos , we rewrite Eq. (1) in polar coordinates (, ), Ω,�� = ( ∈ (0, ∞), ∈ (0, 2)) - ∂ ∂ + 2 Λ(, ) - - ∂ ∂ + 2 Λ(, ) - Ψ(, ) = 0, Λ(, ) = -d2 + (, ), (6) Ψ(, ) = 0, Λ(, ) = -d2 + (, ), (6) (︂ 1 ∂ ∂ 1 )︂ d2 2 where for a trimer with pair potentials (, )= (| sin |) + (| sin( - 2/3)|) + (| sin( - 4/3)|), (7) and for a dimer with pair potential in the external field of barrier potentials: (, )= (| sin |) + (| sin( - /4)|) + (| sin( + /4)|). (8) The solution of Eq. (6) is sought in the form of the Kantorovich expansion ∑︁ ∑︁ max Ψ (, ) = (; ) (). (9) =1 ∈ ∈ ∈ ∈ Here () are unknown functions and the orthogonal normalised angular basis functions (; ) 2(Ω) in the interval Ω = (0, 2) are defined as eigenfunctions, corresponding to the eigenvalues of the Sturm-Liouville problem for the equation ∫︁2 (Λ(, )- ()) (; ) = 0, d(; ) (; )= , (10) 0 ∈ ∞ ∈ ∞ where (), = 1, . . . is the set of the real-valued eigenvalues forming a purely discrete spectrum at each value of the parameter (0, + ). For the problems under consid- eration 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) (=/3, /4 or /2) (, )= ( sin ). Task 2. The case of three pair potentials, Eq. (7), in the interval ∈ (0, 2=/3). Task 3. The case of one pair potential and two penetrable or almost impenetrable barrier potentials, Eq. (8), 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. If the pair potential possesses a high peak in the vicinity of the pair collision point, then the solution of the problem (6) will be considered in the half-plane Ω, = ( (min, ), [min(), 2 min()]) with the Neumann or Dirichlet boundary condition. Since the potential of the boundary-value problem (10) is symmetric with respect to = , the gerade (; ) = (2 ; ) and ungerade (; ) = (2 ; ) solutions, satisfying the Neumann or the Dirichlet boundary condition respectively, will be considered separately in the interval ∈ ∈ [min(), ]. The parametric angular basis functions with successive numbers = 1, . . . , 0 are referred to as cluster states with () < 0 and those with ;;? 0 + 1 as pseudostates with () > 0 corresponding to discrete and continuous spectrum of (4) at large values of the parameter , respectively. The system of coupled ODEs in the Kantorovich form has the form [︂ 1 d d () ]︂ max - d d + 2 - () + ∑︁ ∑︁ =1 () () = 0, (11) 1 d d () = () + d () + () d . (12) Here the potential curves (terms) () are eigenvalues of the BVP (10) and the effective potentials () = -��(), () = () are given by the integrals calculated using the above symmetry on reduced intervals ∈ [0, 2]: () = - 2 d (; ) ∫︁ ∫︁ d (; ) , () = 2 d��(; ) d (; ) d�� . (13) ∫︁ ∫︁ 0 d 0 d For Task 3 the effective potentials ˆ () = () + () are sums of (), calculated using the potential curves and the parametric basis functions of Task 1, and the integrals of barrier potentials () multiplied by the basis functions () = ∫︁ ∫︁ d��(; )( ( sin( - /4)) + ( sin( + /4))) (; ). 0 { } { } As an example, we calculated with a required accuracy the parametric basis functions of BVP (10) and the effective potentials (13) for the models of Be2 dimer and Be3 trimer in collinear configuration by means the FEM using the programme ODPEVP [4]. 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 are shown in Figs. 1, 2, and 3. 1. (b) Figure 1. The potential curves of Be2 (in K, 1K=0.18˚A-2), i.e., the energy eigenvalues depending upon the parameter (in ˚A): (a) () and (b) ˜ = ()/2. Here = 1, . . . , 10 In the FEM the eigenfunctions ��(; ) of the problem (10) are approximated by a finite sum of local functions () with coefficients ( , ) specified in the nodes of the finite element grid [4] , , (; ) = ∑︁ ( (; ) = ∑︁ ( , , ) (). (14) =0 Gusev A. A. et al. Algorithms for Solving the Boundary-Value Problems for . . . Gusev A. A. et al. Algorithms for Solving the Boundary-Value Problems for . . . (a) (b) (c) (d) Figure 2. Be+Be2: The potential curves of Be3 (in K), i.e., the energy eigenvalues depending upon the parameter (in ˚A): (a) () and (b) ˜ = ()/2, (c) the isolines of 2D potentials of Be3 trimer, and (d) the diagonal effective potentials (). Here = 1, . . . , 10 (a) (b) (c) (d) 61 61 Figure 3. The effective potentials (13) (a) -1(), (b) -1(), (c) 1(), (d) 1(). Here = 2, . . . , 10 =0 =0 The functions { ()}�� , = , form a basis in the space of polynomials of the -th order. After substituting the expansion (14) into the variational functional corresponding to BVP (10) and minimizing it [13, 14] we obtain the generalized eigenvalue problem A = B. (15) Here A is the stiffness matrix; B is the positive definite mass matrix; is the vector approximating the solution on the finite-element grid; and is the corresponding eigenvalue [4]. The required accuracy of the calculated eigenvalues, eigenfunctions and their derivatives with respect to parameter, and effective potentials is provided by the estimations given below. Let (), (; ) be the exact solution of (10) and , be the numerical solution of (15). Then for a bounded positively defined operator Λ(; ) the following estimates are valid [13] 2 ⃦⃦ ⃦⃦ 2 ⃦⃦ ⃦⃦ | () - | :( 1�� , ⃦ (; ) - ⃦ :( 2+1, 1 > 0, 2 > 0, (16) where ||(; )||0 = where ||(; )||0 = (; )d; is the maximal grid step, is the order of finite (; )d; is the maximal grid step, is the order of finite 2 ∫︀max 2 min min 0 min min elements, is the number of the corresponding eigensolution, and the constants 1 and 2 do not depend on the step . It is necessary to mention that the second estimate of Eq. (16) is valid also for the solution ∂ (; )/∂ of the problem [4] ∂ (; ) ∂Λ(, ) ∂ () (︂ )︂ (︂ )︂ (Λ(, )- ()) = - - (; ), ∂ ∫︁ ∫︁ 2 ∂ ∂ (17) ∂ (; ) d (; ) =0. ∂ 0 This fact guarantees the same accuracy for eigenfunctions and their derivatives within the present method. Theorem 1. Let Λ(; ) be a bounded positively defined operator on the finite interval ∈ ∈ ∈ ∈ (min, max). Let ∂ (, )/∂ be also bounded for each value of the parameter . Then for the exact solutions, ∂ ()/∂, ∂ (; )/∂ 2, from (17) and the potential matrix elements, (), (), from (13), and the corresponding numerical values, ∂/∂, ∂/∂ ∈ 1, and �� , , the following estimates are valid [4] ⃒⃒⃒ ∂ () ∂ ⃒⃒⃒ 2 ⃦⃦⃦ ∂ (; ) ∂ ⃦⃦⃦ +1 ∂ ⃒ ∂ ⃒ ⃒⃒ ∂ - :( 3 , , ⃦ ⃦ ⃒ ⃦ ∂ - ∂ ⃦⃦0 :( 4 , (18) ⃒⃒ () - ⃒⃒ :( 52, ⃒⃒ () - ⃒⃒ :( 62, where is the maximal grid step, is the order of finite elements, , are the numbers of the corresponding solutions, and the constants 3, 4, 5 and 6 do not depend on the step . For this model the eigenvalues and the hyperradial components of 2D eigenfunctions of the BVP for the set of ODEs (11) with Neumann boundary conditions were calcu- lated with the required accuracy by means the FEM using the program KANTBP [9]. The discrete energy spectrum of the dimer Be2 calculated of the grid Ω��(1.8, 10) = { } { } 1.8(24)3(10)4(5)5(10)10 and the set of binding energies of the trimer Be3 calculated of the grid Ω = 4.1(20)7(10)10 with the eighth-order Lagrange elements ( = 8) are shown in Table 1, and the components of the trimer eigenfunctions (9) are shown in Fig. 4. () () Figure 4. Components ,��=,(, ) ≡ () of gerade () and ungerade () bound states of the trimer Be3 with total energy in ˚A-2 4. Algorithms for calculating the asymptotes of parametric angular basis functions and effective potentials Algorithm for calculating the cluster parametric angular basis functions and effective potentials of the lower part of the discrete spectrum Let us calculate the solution of the Sturm-Liouville problem (10) at large (Λ(; )- ()) (; )≡ (Λ(; )- ()) (; )≡ - ∂2 + ( sin )- () - ∂2 + ( sin )- () (; )=0. (19) (; )=0. (19) (︂ ∂2 2 )︂ Using the new variable ′ defined as = ′/, ′ = arcsin(/) we get - ∂′2 + ( sin( /)) - - ∂′2 + ( sin( /)) - (︂ ∂2 ′ ())︂ ′ 2 2 2 2 ( ; ) = 0. (20) ( ; ) = 0. (20) + ∆ )) = ( ) + + ∆ )) = ( ) + , (21) ! , (21) ! In the argument of the potential we add and subtract ′ and expand the potential in Taylor series in the vicinity of ′, ( sin( /)) = ( ( sin( /)) = ( ′ ′ ′ ′ ∑︁ d (′) (∆′) =1 d′ =1 d′ =1 d′ =1 d′ - ≪ - ≪ where the small correction ∆′ = sin(′/) ′ 1 is presented in the form of Taylor series ∆′ = ∑︁=1 (-1) (2 + 1)! ′2+1 2 . Then the Sturm-Liouville problem (19) is reduced to (︃ ∂2 ′ ∑︁ ()(′) ())︃ ′ - ∂′2 + ( ) + =1 ∫︁ ∫︁ ′max 2 - 2 ( ; ) = 0, (22) ⟨()| ()⟩ ≡ ′0 d′((′; ) (′; ) = . Here first terms ()(′) of the asymptotic expansion of ( sin(′/)) read as ′3 d (′) - - (1)(′)= (︂ (︂ 6 , )︂ )︂ d′ (2)(′)= ′5 360 5′ d2 (′) d′2 d (′) +3 , d′ )︂ )︂ (23) (3)(′)=- ′7 45360 (︂35′ 2 d3 (′) d′3 +63′ d2 (′) d′2 d (′) +9 , )︂ )︂ d′ (4)(′)= ′9 5443200 (︂175′ 3 d4 (′) d′4 +630′ 2 d3 (′) d′3 +369′ d2 (′) d′2 +15 d (′) . d′ Note that !�� ()(′) are the derivatives of the potential of the BVP (22) with respect to the parameter -2. So, we apply the modified version of the program ODPEVP for calculating the parameter derivatives of the solution up to the given order to determine the asymptotic expansion of the cluster eigenfunctions and eigenvalues. ∈ ∈ In the framework of FEM using similar expansions (14) we reduce the eigenvalue problem on the given finite-element grid Ω′ (′min, ′max) in the finite interval ′ (′min, ′max) to a generalized eigenvalue algebraic problem with respect to the eigenvalues ∈ { }=1 and the corresponding eigenvectors () = {(()) }=1 of the order . In { } { } the considered example Ω′ (1.8, 10) = 1.8(24)3(10)4(5)5(10)10 with Lagrange elements ≈ ≈ of the eighth order NPOL=p=8 (the numbers of elements between the nodes given in parentheses), providing the accuracy (2) for the eigenvalues and (+1) 10-8 for the eigenfunctions, being the maximal grid step. equals the number of grid nodes, which in the considered example was = 393. ()() - ()() = 0, (24) (()) () = 1. (25) Here the matrix ��() is presented in the form of inverse power series max max () = (0) + () = (0) + (). (26) (). (26) ∑︁ -2 =1 ! =1 ! =1 ! =1 ! We choose as unperturbed the matrix operator ��(0) corresponding to the differential ∂ ∂ 2 one, - ∂′2 + (′), and numerically solve the corresponding algebraic problem (0)(0) - (0)(0) = 0, (0) (0) = 1. (27) The solution ��(), () is sought in the form of inverse power series () = (0)() + max ∑︁ -2 ∑︁ -2 ! =1 (), () = (0) + max ∑︁ -2 ∑︁ -2 ! =1 (). (28) The substitution of Eq. (28) into Eq. (22) leads to the system of inhomogeneous algebraic equations for the corrections ��() and (): (��(0)() - (0)()) + (��()(0) - ()(0)) = (), -1 ≡ ∑︁ ! ((-)() - ()(-)), ≡ 0, () =1 !( - )! () ∂(0) () ∂(0) (1) (29) ≡ ∂ , ≡ ∂ . Differentiating the normalization condition (25), we obtain additional conditions for the corrections ��(): -1 -1 2 2 !( - )! !( - )! (0) () ≡ () = - 1 ∑︁ ! (-) (). =1 =1 =1 =1 Multiplying (29) by ��(0) and taking the zero values of the first term and the nor- malization condition into account, we get the formula for the corrections () of the eigenvalues: () () () = (0) ()(0) - (0) . (30) The vector ��() is calculated by solving the system of algebraic equations () ≡ (0)() - (0)() = () () ≡ - ∑︁ ! (()(-) - ()(-)), - - !( )! =1 ∑︁ ! ∑︁ ! - - (-) () = 0, !( )! =0 (31) where the latter equation is a result of differentiation of the normalization condition. Since ��(0) is an eigenvalue of (27), the matrix in Eq. (31) is degenerate. In this case the algorithm for solving Eq. (31) can be written in three steps as follows: Step k1. Calculate the solutions ��() and of the auxiliary inhomogeneous systems of algebraic equations ¯ () = ¯(), ¯ = (32) with the non-degenerate matrix ¯ and the right-hand sides ¯() and ¯ = �� , ( - )( - ) ̸= 0, {︂ {︂ = = 0, = , 0, = , = = , ( - )( - ) = 0, ¯() {︂ (), ̸= , {︂ , = , 0, = , 0, = , where is the number of vector (0) element having the maximal absolute value. Step k2. Evaluate the coefficient ��() () () () = - 1 () = - 1 1 1 2 2 2 2 - ((0) - ) , () = () (0), = (0). (33) Step k3. Evaluate the vector ��() () = () () {︃ - , ̸= , {︃ - , ̸= , (), = . (34) The execution of the above procedure yields the required asymptotic expansion for the eigenvalues () = ��(0)() + max -2() = (0)() + ∑︁ -2 () (35) 2 ∑︁ ∑︁ =1 max max =1 ! and the corresponding expansion of the eigenfunctions in the nodes ′ of the given grid max max (′; ) = ((0)) + (′; ) = ((0)) + (()) . (()) . ∑︁ -2 =1 ! =1 ! =1 ! =1 ! Following Eq. (14) we present it in the form of a piecewise-polynomial function (+1) (+1) (′ ∈ [′, ′ ]; ) = ∑︁ (′ ′=0 ′=0 ; ) ∏︁ ′ - ′′′ . (36) ′=0 ′′=0,′′̸=′ ′=0 ′′=0,′′̸=′ +′ +′ ′′ ′′ - ′′′ - ′′′ We calculate the discrete spectrum solutions (; ) of the problem (19) on the grid = ′/, related to the solutions (′; ) of the problem (20) as ( = ′/; ) = √ (′; ). The corresponding piecewise-polynomial functions are ( ∈ [, (+1)]; )= (=′/ ∈ [′/, ′(+1)/]; ) = -′ = ∑︁ ( ��+ ′ ; ) ∏︁ -′′ ′ - ′′ = ∑︁ ( + ′ ; ) ∏︁ ′ ′′ = -′ ′=0 ′′=0 ′′̸=′ ′=0 ′′=0 ′′̸=′ ′ ′′ √ -′ = ∑︁ (′+′ ; ) ∏︁ ′=0 ′=0 ′′ = ′ ′ -′ ′′ ′′=0 ′′=0 ′′ ′ =√ ∑︁ (︃((0)) ∑︁ -2 (()) )︃ ∏︁ -′′′ . (37) ′=0 + + =1 ! ′′=0 ′′ -′′′ The nodes ′′ in Eq. (37) are the same as in Eq. (36). ′′ ′ { } { } ∂ ∂ The calculation of the matrix elements is similar to the derivation of Newton-Cotes formulae [15]. In each of the subintervals of the nonuniform finite element mesh the calculated values of functions and their derivatives are approximated, in the considered case of Ω[1.8/, 10/] = 1.8/(24)3/(10)4/(5)5/(10)10/ by a Lagrange polynomial of the order = 8, which provides the same relative accuracy of (+1). The obtained approximations of functions are polynomials of , explicitly depending on . By analytical differentiation of Eq. (37) with respect to the parameter we arrive at the explicit expansion of the derivative ∂ (;��) . Integrating the products of functions (; ) and/or - - their parametric derivatives in each subinterval of the length = (+1) of the nonuniform finite element mesh Ω and summing the obtained results, we get the required asymptotic expansion of matrix elements (13) for cluster states, , = 1, . . . , 0: () = ∑︁ ∑︁ =1 (2 1) - - 2-1 , () = ∑︁ ∑︁ =1 (2) 2 . (38) For the beryllium trimer in collinear configuration, the coefficients ��() of the expansion (35) and the coefficients () of the expansion (38) at = are presented in Table 2, and the first coefficients of the expansions (38) are ⎜ ⎜ (1) ⎛ ⎜ 0 10.3201 -3.81790 1.83170 -.903311 ⎟ ⎟ ⎞ ⎞ -10.3201 0 12.2419 -5.54402 2.68986 ⎟ = ⎜⎜⎝ ⎛ ⎜ ⎜ ⎜ 3.81790 -12.2419 0 11.5096 -5.30745 -1.83170 5.54402 -11.5096 0 7.99437 0.903311 -2.68986 5.30745 -7.99437 0 0 19.6114 -3.63427 1.58999 -6.31271 -19.6114 0 32.2355 -12.5531 22.0625 ⎟⎟⎠ , ⎞ ⎟ ⎟ ⎟ = ⎜ = ⎜ (3) ⎜⎝ ⎛ ⎜ ⎜ ⎜ 3.63427 -32.2355 0 52.0431 -53.8134 -1.58999 12.5531 -52.0431 0 105.275 6.31271 -22.0625 53.8134 -105.275 0 0 68.2366 -6.68774 -14.8770 128.505 -68.2366 0 91.6918 27.6600 -356.227 ⎟⎟⎠ , ⎞ ⎟ ⎟ ⎟ = ⎜ = ⎜ (5) ⎜⎝ 6.68774 -91.6918 0 68.9065 595.104 14.8770 -27.6600 -68.9065 0 -523.035 -128.505 356.227 -595.104 523.035 0 ⎟⎟⎠ , 68 68 Coefficients ��() and () of the expansions (35) and (38) Table 2 1 2 3 4 5 (0) -193.06601252 -119.39267226 -63.338854932 -24.904560537 -4.0897890782 (1) -127.73059638 -317.30568182 -408.25519644 -385.66879485 -223.26092504 (2) -215.85831875 -672.02680456 -1198.7991074 -1892.8857386 -3079.0060449 (3) -667.46086524 -2093.0917304 -3656.2842793 -5140.5360777 -3803.8470973 (4) -2590.0140097 -8317.9479477 -15062.850569 -23930.063302 -63419.196580 (5) -11287.527768 -37281.527942 -69636.528276 -107755.44601 164277.292590 (2) 127.98059638 317.55568182 408.50519643 385.91879484 223.51081140 (4) 414.47883932 1204.0042402 2009.0880648 3041.7739524 5091.1634806 (6) 1912.4702415 5431.4025830 8482.5056376 10047.915436 -664.09784640 (8) 9895.5812579 28061.019351 44415.159972 62195.582167 189759.82059 (10) 54127.043138 154284.86093 246059.16425 311539.17441 152663.97768 (1)+ (2) 0.2499999984 0.2499999944 0.2499999948 0.2499999953 0.2498863674 Bulletin of PFUR. Series Mathematics. Information Sciences. Physics. No 4, 2016. Pp. 56-76 Bulletin of PFUR. Series Mathematics. Information Sciences. Physics. No 4, 2016. Pp. 56-76 Table 3 Convergence of the expansion (35) to the numerical values (NUM) at = 20 1 2 3 4 5 (0)/2 -193.06601252 -119.39267226 -63.338854932 -24.904560537 -4.0897890782 +(1)/4 -193.38533901 -120.18593646 -64.359492923 -25.868732525 -4.6479413908 +(2)/6 -193.38668813 -120.19013663 -64.366985417 -25.880563060 -4.6671851786 +(3)/8 -193.38669856 -120.19016933 -64.367042547 -25.880643381 -4.6672446137 +(4)/10 -193.38669866 -120.19016966 -64.367043135 -25.880644316 -4.6672470910 +(5)/12 -193.38669866 -120.19016966 -64.367043142 -25.880644327 -4.6672470749 (NUM) (NUM) = /2 = /6 -193.38669866 -193.38669866 -120.19016966 -120.19016966 -64.367043142 -64.367043142 -25.880644327 -25.880644327 -4.6672471622 -4.6672471630 ⎜ ⎜ (2) ⎛ ⎜ 127.981 -67.2820 -85.3811 73.3913 -44.4889 ⎟ ⎟ ⎞ ⎞ -67.2820 317.556 -161.487 -40.5195 46.8714 ⎟ = ⎜⎜⎝ ⎛ ⎜ ⎜ ⎜ -85.3811 -161.487 408.505 -231.091 45.2520 73.3913 -40.5195 -231.091 385.919 -215.798 -44.4889 46.8714 45.2520 -215.798 223.511 414.479 -142.333 -526.641 464.604 -624.169 -142.333 1204.00 -420.171 -812.366 1119.09 ⎟⎟⎠ , ⎞ ⎟ ⎟ ⎟ = ⎜ = ⎜ (4) ⎜⎝ ⎛ ⎜ ⎜ ⎜ -526.641 -420.171 2009.09 -947.285 -464.551 464.604 -812.366 -947.285 3041.77 -2605.05 -624.169 1119.09 -464.551 -2605.05 5091.16 1912.47 -670.689 -2112.77 815.251 4153.61 -670.689 5431.40 -1689.40 -3295.46 -646.435 ⎟⎟⎠ , ⎞ ⎟ ⎟ ⎟ = ⎜ = ⎜ (6) ⎜⎝ -2112.77 -1689.40 8482.51 -2016.27 -12077.1 815.251 -3295.46 -2016.27 10047.9 12431.9 4153.61 -646.435 -12077.1 12431.9 -664.074 ⎟⎟⎠ . The convergence of the expansion (35) to the numerical values (NUM) is presented in Table 3. One can see that the first six terms of the asymptotic expansion provide the accuracy of 11 significant digits. Remark 1. The effective potentials ��(1)+ (0)=(1/4)-2 (presented, e.g., in the 1/2 √︀- 0 1/2 √︀- 0 last line of Table 2) lead to the asymptotic cluster fundamental solutions of the ODEs (11), determined by the Bessel functions ( ��+ ), =1, . . . , , while for pseu- dostates the asymptotic fundamental solutions of the ODEs are (√) with the integer =(-0)=1, 2, . . .. Algorithm for calculating the parametric angular functions of pseudostates and the effective potentials - - - - The eigenfunctions of pseudostates () ;;? 0, ( 0) = 1, 2, . . ., are localised beyond the potential well. Then the (0 1)-th node is located at the boundary of the potential well. Here and below we consider the case of = /2 for the atomic dimer illustrated by Fig. 5. Figure 5. Eigenfunctions (; ) corresponding to the eigenvalues () ;;? 0 of the gerade () pseudostates = + 1 = 6, . . . , 10 at = 100 ≈ - ≈ - From this fact the estimate of the eigenvalues for pseudostates () ( 0)2 follows, namely, the eigenvalues of the corresponding BVPs (22) in new variable ′, = ()/2, will be a small quantity (see Fig. 1b). So, the solution ��() ( = ()/2) of the derived equation is sought in the form of a power series at / < 1 ∑︁ ∑︁ max () () = 2+ �� . (39) =1 Then the numerical values of the function ��(; ) = (′) and its derivative ′ ′ { ′1 ′0 ′ ′ ′ } ′ ′ { ′1 ′0 ′ ′ ′ } { } { } ′(; ) = ′(′) on the specified grid Ω = 1 = 0, . . . , = /, . . . , = in the polar system of coordinates are determined via the values of the function () and its derivative ( ) on the grid Ω = , . . . , , . . . , = , found with the same accuracy accepted in the FEM scheme chosen above in the form of the power series of the small parameter : max () (1) () max () (1) () ( )=(0)+∑︁ ( , . . . , ) , ′ ( )=(0)+∑︁ ( , . . . , ) , (40) =1 =1 =1 =1 =1 =1 =1 =1 using the Runge-Kutta method, in which the all terms contain the power max + 1 of 1/ and the higher ones are neglected. The expansion coefficients ��() ≡ ()(′ ) and () ≡ ()(′ ), calculated at the grid nodes ′ for the BVP (22) with the potential (′) defined in Eq. (5) are presented in Fig. 6. Figure 6. Expansion coefficients ()((1), . . . , ()) and ()((1), . . . , ()), =0, 2, 3, 4, 5, = 1 with () given by Eq. (44) calculated at the nodes ′ of the grid Ω ′ One can see that in the vicinity of the potential well the corrections to the eigenfunctions are small, and at large ′ they become essential. The coefficient ��(0), the derivative of the wave function with = 0, exponentially tends to a constant for > 5.5. From these observations the condition for choosing follows. However, to avoid analytical calculations of the exponential terms in the effective potentials (13) between the weakly bound cluster states and pseudostates, it is sufficient to choose ′ = 10. The interval 0 :( :( /2 is divided into two subintervals by the point = /: | | | | 0 < :( and /2 > > . In the calculations the point was chosen from the condition ( > ) < , where is a preassigned number, and the left-hand boundary of the interval 0 = 0. In the case of a high barrier, at the pair collision point, when the eigenfunctions in its vicinity are close to zero, the left boundary of the interval changes, 0 = 0/ > 0. The eigenfunctions (; ) are calculated in the form (41) (41) ∫︁ ∫︁ {︃ ()�� (; ), 0 :( :( , sin sin (; ) = ()√︁ 2 {︁cos}︁ (√︀ ()( - /2)), < :( /2, /2 2 0 d((; ))2 = 1. Here () and () are the normalisation factors, and ��(; ) is determined from the numerical solution () in Cartesian coordinates using the transformation = /. From the continuity of the eigenfunctions and their derivatives, ( - 0; ) = ( + 0; ), d ( d d - 0; ) = d ( (42) + 0; ), we get the equation for the eigenvalue (): {︃ tan(√︀()(- )) even }︃ 2 2 √︀() -cot(√︀ ()( - )) odd =0, (43) (43) - - 2 2 ′ (; ) ′ () = = . (; ) () Substitute (40) into (43), and then substitute (39) into the resulting equation. Ex- panding both sides of the equation in inverse powers of , we arrive at the system of linear equations, from which the expansion coefficients (), and then the coefficients () and () are determined. - - Since the values of the function ��(; ) and its derivative ′ (; ) on the grid Ω�� are known, for the calculation of the first integral we use the quadrature formula of the Newton-Cotes type. The second integral is calculated analytically using the expansion (39). We have the analytical expression in the interval ��() < :( /2, and the explicit dependence of its values upon the parameter on the grid Ω. For the considered potential (5) we get the asymptotes of the and potential curves at = 0 for / < 1: ∑︁ ∑︁ max () ()=2+ �� ; for -parity at odd , or -parity at even . (44) =1 The first terms of these expansions are expressed as (1)=4.48622572, (2)=15.0946662, (3)=-88.1643242+0.29884084, (4)=-770.503442+3.35166834, (5)=-13803.8542+162.148504-0.29227346, (6)=-68926.7632+1584.07824-4.43293086. The calculated eigenvalues in comparison with the numerical solution obtained by means of the program ODPEVP [4] are presented in Table 4. Table 4 Convergence of the expansion (44) at =50 and the numerical results (NUM) gerade (-parity) 1 3 5 7 2 1.000000 9.00000 25.00000 49.00000 +(1)/ 1.089724 9.80752 27.24311 53.39650 +(2)/2 1.095762 9.86186 27.39405 53.69235 +(3)/3 1.095059 9.85570 27.37792 53.66353 +(4)/4 1.094936 9.85464 27.37517 53.65878 +(5)/5 1.094893 9.85428 27.37437 53.65775 +(6)/6 1.094888 9.85425 27.37432 53.65774 (NUM) 1.094887 9.85424 27.37431 53.65776 ungerade (-parity) 2 4 6 8 2 4.000000 16.00000 36.000000 64.0000 +(1)/ 4.358898 17.43559 39.230082 69.7423 +(2)/2 4.383049 17.53219 39.447445 70.1287 +(3)/3 4.380266 17.52152 39.425152 70.0934 +(4)/4 4.379781 17.51968 39.421409 70.0877 +(5)/5 4.379613 17.51911 39.420448 70.0868 +(6)/6 4.379597 17.51906 39.420407 70.0868 (NUM) 4.379592 17.51905 39.420408 70.0869 - - - - The described algorithm is implemented in the Maple-Fortran environment. The asymptotic expansions, obtained using it at = 50, coincide with the numerical solution given by the finite element method to 5-6 significant digits for the eigenvalues and to 4-5 significant digits for the eigenfunctions. The asymptotes of the effective potentials (13) between the states 1 = 0 and 2 = 0 of the same ( or ) parity at 0 = 5, , = 0 + 1, . . . for / < 1 have the form: ∑︁ ∑︁ max (+2) max (+4) 12 ()= 12 +2 =0 , 12 ()= 12 +4 =0 . (45) ∑︁ ∑︁ First terms of these expansions read as (2) 12 =2.2431128 , (3) 12 =5.0315554 , 12 (2-2) 12 (2-2) 1 2 12 1 2 12(2-2) 1 1 2 2 (4) 12 =188.67822 (2-2) +0.22413 1 2 , 1 1 2 2 (2-2) (4) 12(2+2) (5) 12(2+2) 1 2 1 2 =45.14538 +69.9350912, =10.06311 1 2 , =10.06311 1 2 , 12 (2-2)2 12 (2-2)2 1 2 1 2 1 2 1 2 1 2 12(2+2) 1 2 33 (6) =-1642.273 +8.044 , +476.76212 1 1 2 2 1 1 2 2 12 (2-2)2 (2-2)2 (4) =0.6289444+2.06914422, (5) =2.8215867+79.2177462, 2 4 2 4 11 1 11 1 11 11 (6) = - 102.642068 + 138.328874 + 0.826991 . 1 1 1 1 Using Eq. (35) and Eq. (41), we get the asymptotic expansions for () and () between the cluster states =1, . . . , 0 and pseudostates -0=1, 2, . . . , for /<1 ∑︁ ∑︁ max (+5/2) max (+7/2) ()= �� +5/2 =0 , ()= �� +7/2 =0 , (46) ∑︁ ∑︁ The first terms of these expansions read as 1 1 1 1 1 1 (5/2)=-0.428911, (7/2)=-1.443145, (9/2)=-25.01600+0.1837913, 2 2 2 2 2 2 (5/2)= 1.273900, (7/2)=+4.286254, (9/2)= 75.66328-0.5466573, 3 3 3 3 3 3 (5/2)=-2.497511, (7/2)=-8.403301, (9/2)=-152.3919+1.0754213, 4 4 4 4 4 4 (5/2)= 3.668834, (7/2)=+12.34441, (9/2)= 231.9007-1.5990053, 5 5 5 5 5 5 (5/2)=-4.098659, (7/2)=-13.79066, (9/2)=-247.9253+2.0188733, 1 1 1 1 1 1 (7/2)= 22.53006, (9/2)=+77.24936, (11/2)=1551.8678-9.8358923, 2 2 2 2 2 2 (7/2)=-26.57543, (9/2)=-93.70381, (11/2)=-2203.032+11.896573, 3 3 3 3 3 3 (7/2)=-12.19890, (9/2)=-32.64199, (11/2)=262.11512+4.4004933, 4 4 4 4 4 4 (7/2)= 85.06821, (9/2)=+273.8820, (11/2)=4387.8346-35.919253, 5 5 5 5 5 5 (7/2)=-120.4823, (9/2)=-391.5948, (11/2)=-8001.428+54.831113. 5. Conclusion The model for beryllium trimer in collinear configuration is formulated as a 2D boundary-value problem for the Schro¨dinger equation in polar coordinates. This problem is reduced using the Kantorovich expansions to the boundary-value problem for a set of second-order ordinary differential equations. The symbolic-numeric algorithms are proposed and implemented in Maple to evaluate the asymptotic expansions (35), (38), (45) and (46) of the parametric BVP eigensolutions and the effective potentials () in inverse powers of . It can be used for calculation of the asymptotes of the fundamental solutions of the system of second-order ODEs at large values of [10] and construction of asymptotic states of three atomic scatterring problem. The proposed approach can be applied to the further analysis of quantum transparency effect [10, 11], quantum diffusion [16-18] and the resonance scattering in triatomic systems using modern theoretical and experimental results [19-21] and algorithms and programs [4-9]. References 1. V. Efimov, Few-Body Physics: Giant Trimers True to Scale, Nature Phys. 5 (2009) 533-534. 2. M. Zaccanti, B. Deissler, C. D’Errico, M. Fattori, Jona-Lasinio, S. M. Muller, G. Roati, M. Inguscio, M. G., Observation of an Efimov Spectrum in an Atomic System, Nature Phys. 5 (2009) 586-591. 3. J. Voigtsberger, et al., Imaging the Structure of the Trimer Systems 4He3 and 3He4He2, Nature Commun. 5 (2014) 5765. 4. 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. 5. A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, A. G. Abrashkevich, POTHEA: A Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Parametric Self-Adjoined 2D Elliptic Partial Differential Equation, Comput. Phys. Commun. 185 (2014) 2636-2654. 6. A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, A. G. Abrashkevich, Description of a Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Coupled Parametric Self-Adjoined Elliptic Differential Equations, Bulletin of Peoples’ Friendship University of Russia. Series “Mathematics. Information Sciences. Physics” (2) (2014) 336-341. 7. O. Chuluunbaatar, A. A. Gusev, A. G. Abrashkevich, A. Amaya-Tapia, M. S. Kaschiev, S. Y. Larsen, S. I. Vinitsky, KANTBP: A Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adiabatic Approach, Comput. Phys. Commun. 177 (2007) 649-675. 8. O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich, KANTBP 2.0: New Version of a Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adiabatic Approach, Comput. Phys. Commun. 179 (2008) 685-693. 9. 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 10. S. I. Vinitsky, A. A. Gusev, O. Chuluunbaatar, L. Hai, V. Derbov, P. Krassovitskiy, 1. Go´´zd´z, Symbolic Numerical Algorithm for Solving Quantum Tunneling Problem of a Diatomic Molecule Through Repulsive Barriers, Lect. Notes Comp. Sci. 8660 (2014) 472-490. 11. P. M. Krassovitskiy, F. M. Pen’kov, Contribution of resonance tunneling of molecule to physical observables, J. Phys. B 47 (2014) 225210. 12. J. Wang, G. Wang, J. Zhao, Density Functional Study of Beryllium Clusters with Gradient Correction, J. Phys. Cond. Matt. 13 (2001) L753-L758. 13. G. Strang, G. J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, New York, 1973. 14. K. J. Bathe, Finite Element Procedures in Engineering Analysis, Englewood Cliffs, Prentice Hall, New York, 1982. 15. I. S. Berezin, N. P. Zhidkov, Computing Methods, Vol. 1, Pergamon Press, Oxford, 1965. 16. A. A. Gusev, L. L. Hai, Algorithm for Solving the Two-Dimensional Boundary Value Problem for Model of Quantum Tunneling of a Diatomic Molecule Through Repulsive Barriers, Bulletin of Peoples’ Friendship University of Russia. Series “Mathematics. Information Sciences. Physics” (1) (2015) 15-36. 17. E. Pijper, A. Fasolino, Quantum Surface Diffusion of Vibrationally Excited Molecular Dimers, J. Chem. Phys. 126 (2007) 014708. 18. A. A. Gusev, A. G´o´zd´z, V. L. Derbov, S. I. Vinitsky, O. Chuluunbaatar, P. M. Krassovitskiy, Models of Resonant Tunneling of Composite Systems Through Repulsive Barriers, JINR News (1) (2014) 22-26. 19. A. V. Mitin, Ab Initio Calculations of Weakly Bonded He2 and Be2 Molecules by MRCI Method with Pseudo-Natural Molecular Orbitals, Int. J. Quantum Chem. 111 (2011) 2560-2567. 20. J. M. Merritt, V. E. Bondybey, M. C. Heaven, Beryllium dimer-caught in the act of bonding, Science 324 (2009) 1548-1551. 21. K. Patkowski, V. Sˇpirko, K. Szalewicz, On the Elusive Twelfth Vibrational State of Beryllium Dimer, Science 326 (2009) 1382-1384.

About the authors

A A Gusev

Joint Institute for Nuclear Research

Email: gooseff@jinr.ru
Dubna, Russia

Ochbadrah Chuluunbaatar

Joint Institute for Nuclear Research

Email: chuka@jinr.ru
Institute of Mathematics, National University of Mongolia, Ulan-Bator, Mongolia Dubna, Russia

S I Vinitsky

Joint Institute for Nuclear Research

Email: vinitsky@theor.jinr.ru
RUDN University (Peoples’ Friendship University of Russia), Moscow, Russia Dubna, Russia

V L Derbov

Saratov State University

Email: derbov@sgu.ru
Saratov, Russia

References