Algorithms for solving the boundary-value problems for atomic trimers in collinear configuration using the Kantorovich method Текст научной статьи по специальности «Физика»

UDC 519.632.4

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

A. A. Gusev*, O. Chuluunbaatar*^, S. I. Vinitsky**, V. L. DerbovS

* Joint Institute for Nuclear Research, Dubna, Russia ^ Institute of Mathematics, National University of Mongolia, Ulan-Bator, Mongolia * RUDN University (Peoples' Friendship University of Russia), Moscow, Russia S Saratov State University, Saratov, Russia

The model of atomic trimers with molecular pair interactions for collinear configuration is formulated as a 2D boundary-value problem (BVP) in the Jacobi and polar coordinates. The latter is reduced to a 1D BVP for a system of second-order ordinary differential equations (ODEs) by means of the Kantorovich method using the expansion of the desired solutions over a set of angular basis functions, parametrically dependent on the (hyper)radial variable. The algorithms for solving the 1D parametric BVP by means of the finite element method (FEM) and calculating the asymptotes of the parametric angular functions and effective potentials of the system of ODEs at large values of the parameter are presented. The efficiency of the algorithms is confirmed by comparing the calculated asymptotic solutions and effective potentials with those of the parametric eigenvalue problem obtained by applying the FEM at large values of the parameter. The applicability of the algorithms is demonstrated by calculating the asymptotic expansions of the parametric BVP solution, effective potentials and sets of binding energies for the beryllium trimer in the collinear configuration.

Key words and phrases: boundary-value problems, the Kantorovich method systems of second-order ordinary differential equations, finite element method

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, 1701-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 M and the coordinates Xi G R1, i = 1, 2,3, coupled via the pair potential -(|Xi — Xj|) i,j = 1,2, 3. Performing the change of variables at cyclic permutation (a,@,^) = (1, 2, 3):

____Xa + Xg — 2x^ V2. ,

X = X(ap) = xa - Xf3, y = =--L, Xo = + x2 + X3),

we arrive at the Schrodinger equation for the wave function in the center-of-mass system {Xi G R1 Ixi + X2 + X3 = 0}

( d2 d2 M ~ ~ \

(—B? - W2 + F— *))*('■<■*) = 0

(1)

In the case of a diatomic molecule with identical nuclei coupled via the pair potential y(|xi — x2\) and moving in the external potential field y6(|Xi — X3I), i = 2,1 of the third atom having the infinite mass, the same equation (1) is valid for the variables

x = x1 — x2, y = x1 + x2,

the origin of the coordinate frame being placed on the infinite-mass atom, x3 = 0.

Here the potential function for a trimer with the pair potentials (below this case is referred to as Task 2),

yOr,y)=y(M)+y

"-V3y

2

+ V

+ V3y

2

(2)

or the potential function for a dimer in the field of barrier potentials (below this case is referred to as Task 3)

x-y

x+y

2

(3)

is symmetric with respect to the straight line x = 0 (i.e., x1 = x2), which allows one to consider the solutions of the problem in the half-plane x ^ 0. Using the Dirichlet or Neumann boundary condition at x = 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 x = 0 is exponentially small and can be considered in the half-plane x ^ xmin. In this case setting the Neumann or Dirichlet boundary condition at ¡Emm gives only a minor contribution to the solution. The equation, describing the molecular subsystem (dimer), has the form

( d2 M -, . V

(x) = 0.

(4)

We assume that the dimer has the discrete spectrum, consisting of a finite number n0> 1 of bound states with the eigenfunctions <(x), j = 1,n0 and eigenvalues ij = —|ij |< 0, and the continuous spectrum of eigenvalues i > 0 with the corresponding eigenfunctions <e(x). 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 ij > 0, j = 1 + n0,... 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 M/2 = 4.506 Da of the nuclei [10,11]

M

V (x) = ^y(x), V(x) = D{exp[—2(x — Xeq )a] — 2exp[-(x — xeq)a]}. (5)

Here a = 2.96812 A-1 is the potential well width, xeq = 2.47 A is the average distance between the nuclei, and D = 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 £i = (M/h2)ii, i = 1,... ,n0 = 5 presented in Table 1. The parameter values are

determined from the condition (i2 — i1)/(2nHc) = 277.124cm-1, 1K/(2^ft c)=0.69503476 cm-1.

Table 1

The discrete spectrum energies of dimer Be2 and the binding energy Eb = —(E — Ea) of gerade (g) and ungerade (w) states of trimer Be3 counted of Ea = £i = —193.06A-2 = —1044K dimer energy Be2 calculated in grid

np = 4.1(20)7(10)10 at jmaX = 10

- £i E1, Ei,u

-èi=1044.879 649K -£2=646.157 093K - £3=342.791 979K - £4=134.784305K -£5=22.134 073K Elg = 196.02À-2 = 1060.86K E^g = 142.37À-2 = 770.51K EÏg = 93.95À-2 = 508.50K EÏg = 52.77À-2 = 285.63K = 32.32À-2 = 174.95K E^>g = 22.31À-2 = 120.75K E^q = 5.14À-2 = 27.87K E i,u = 107.52À-2 = 581.90K = 67.41À-2 = 364.84K E^u = 34.60À-2 = 187.28K E^tU = 11.79À-2 = 63.83K Ei,a = 0.8À-2 = 4.4K

To solve the discrete spectrum problem we applied the FEM of the seventh order using the Hermitian interpolation polynomials with double nodes [9]. The grid {x0,... ,xi,..., xn} was used to calculate the values of both the function and its derivatives.

3. Reduction of the BVP using the Kantorovich method

Using the change of variables x = p sin^, y = p cos^, we rewrite Eq. (1) in polar coordinates (p,<p), = (p G (0, <x>),<£ G (0,2^))

/18 8 1 \ d2

(—PdpPdp+P^P) — E) = 0 A(^,P) = — V +P2V(^ (6)

where for a trimer with pair potentials

V(p,p) = V{p\ sinpj) + V(p\ sin(p - 2tt/3)|) + V(p\ sin(p - 4tt/3)\), (7)

and for a dimer with pair potential in the external field of barrier potentials:

V(p,p) = V(p\ sinpj) + Vb(p\ sin(p - n/4)\) + Vb(p\ sin(p + n/4)\). (8)

The solution of Eq. (6) is sought in the form of the Kantorovich expansion

Jmax

( p,p) = Y1 <J (V;p)Xji0 (p). (9)

i=i

Here Xji0 (p) are unknown functions and the orthogonal normalised angular basis functions <j(p; p) G L2(Q) in the interval Q = p G (0, 2n) are defined as eigenfunctions, corresponding to the eigenvalues of the Sturm-Liouville problem for the equation

(A(p,p)-ej (p))0j (p;p) = 0, Jdp0i(p;p)0j (p;p)=%, (10)

0

where £j (p), j = 1,... is the set of the real-valued eigenvalues forming a purely discrete spectrum at each value of the parameter p G (0, +ro). For the problems under consideration the potential function V(p,p) depending on the parameter p can be defined as follows.

Task 1. The case of one pair potential in the intervals p G (0, 2pa) (pa=n/3, n/4 or n/2) V(p,p)=V(p sinp).

Task 2. The case of three pair potentials, Eq. (7), in the interval p G (0, 2pa=n/3).

Task 3. The case of one pair potential and two penetrable or almost impenetrable barrier potentials, Eq. (8), in the interval p G (0,pa = n/2) or in the intervals p G (0, pa = n/4 - e) and p G (pa = n/4 - e,n/2), 0 < e < n/4.

The solutions symmetric with respect to the permutation of two particles satisfy the Neumann boundary condition at p = 0 and p = 2pa. 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 Qv,p = (p G (pmin, ^),p G [pmin(p), 2pa - pmin(p)j) with the Neumann or Dirichlet boundary condition. Since the potential of the boundary-value problem (10) is symmetric with respect to p = pa, the gerade <j(p; p) = <(2pa - p; p) and ungerade (p;p) = - < (2pa - p; p) solutions, satisfying the Neumann or the Dirichlet boundary condition respectively, will be considered separately in the interval p G [pmin(p), paj. The parametric angular basis functions with successive numbers j = 1,... ,n0 are referred to as cluster states with £j(p) < 0 and those with j > n0 + 1 as pseudostates with £j (p) > 0 corresponding to discrete and continuous spectrum of (4) at large values of the parameter p, respectively.

The system of coupled ODEs in the Kantorovich form has the form

1 d p d + £*(p) F

— I~pT~ +--2--F

pdp dp p2

xn0 (p) + (p)^ (p) = 0, (11)

3 = 1

Wi, (p) = Hji(p) + pdpp^(p) + Q^(p)jp. (12)

Here the potential curves (terms) £j (p) are eigenvalues of the BVP (10) and the effective potentials Qij(p) = —Qji(p), H^-(p) = Hji(p) are given by the integrals calculated using the above symmetry on reduced intervals p G [0, 2pa]:

2ip0

r> ( \ ¡a ai u r \ (a (p;p)

Qn(p) = -1 d,p&(p;p)—dp—,Hij(p) = J d'p—dp--dp—. (13)

For Task 3 the effective potentials Wij (p) = Wij (p) + V-j (p) are sums of Wij (p), calculated using the potential curves and the parametric basis functions of Task 1, and the integrals of barrier potentials V^ (p) multiplied by the basis functions

fa

V* (p) = j d(p<i((p; p)(Vb(p sin(p — n/4)) + V6(p sin(p + n/4)))< (p; p).

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 Q^[1.8/p,a] = {1.8/p(24)3/p(10)4/p(5)5/p(10)a} for a = n/2 for Be2 dimer and a = n/6 for Be3 trimer are shown in Figs. 1,2, and 3.

S 0

-4000

-8000

-12000

gvv........... ------ — ---------

-even

----odd

4 6 P

(a)

10

(b)

Figure 1. The potential curves of Be2 (in K, 1K=0.18A-2), i.e., the energy eigenvalues depending upon the parameter p (in A): (a) £j (p) and (b) £j = £j (p)/p2.

Here j = 1,..., 10

In the FEM the eigenfunctions 0(p; p) of the problem (10) are approximated by a finite sum of local functions Nf(p) with coefficients 4>l(pSS r,p) specified in the nodes pps r of the finite element grid [4]

L

#p;p) = £4>\plr ,p)N?(p). (14)

1=0

10о 12 14

р,А

Р> А

\ N f //7\

У(12)з1

tv 4 5

J A ., \ |к/ \ \V/i23)

(а)

(Ь)

(с)

Figure

А): (а)

2. Ве+Вег: The potential curves of Вез (in К), i.e., the energy eigenvalues depending upon the parameter p (in £j(p) and (b) ij=£j(p)/p2, (c) the isolines of 2D potentials of Вез trimer, and (d) the diagonal effective potentials

Hrj(p). Here j = l,...,10

nK A"2

-6 4

-11 о* з

-17

-22 2

-28

-33 1

-ЗЯ

0

iq, 12 14

P, A

iq, 12 14

P, A

P,A

9 10 tV 12 13 14 15 P,A

(a) (b) (c) (d)

Figure 3. The effective potentials (13) (a) Hjj-i(p), (b) Qjj-i(p), (c) Нц(р), (d) Qji(p). Here j = 2,..., 10

О с

w

(U <

> >

(U

r+ >

n о

in

sr

3'

era

c+

tr

(U

И

о

С p

a

9>

$

С

a>

Hi

о cr

The functions { Nf (p)}f= 0, L — np, form a basis in the space of polynomials of the p-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

Ap 4>

— e

(15)

Here Ap is the stiffness matrix; Bp is the positive definite mass matrix; $>h is the vector approximating the solution on the finite-element grid; and e h 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 tj(p), fy(p;p) be the exact solution of (10) and eh, be the numerical solution of (15). Then for a bounded positively defined operator A(p;p) the following estimates are valid [13]

I tj (p) - < C!h2p, fa (p;p) - $ < C2hp+1, d > 0, C2 > 0,

(16)

where ||v(p; p)||0 = j v2(p; p)dp; h is the maximal grid step, p is the order of finite

^min

elements, j is the number of the corresponding eigensolution, and the constants c1 and c2 do not depend on the step h. It is necessary to mention that the second estimate of Eq. (16) is valid also for the solution d^j (p; p)/dp of the problem [4]

(A( p,p)-e j (p))

90 j (p;p)

dp

2 n

d p <M p;p)

(dA(p, p) de3 (p)

V dp 90 j (p;p)_

d p

d p d p 0.

fa (p;p),

(17)

This fact guarantees the same accuracy for eigenfunctions and their derivatives within the present method.

Theorem 1. Let A(p; p) be a bounded positively defined operator on the finite interval p G (pmin,pmax). Let dV(p,p)/dp be also bounded for each value of the parameter p. Then for the exact solutions, dej (p)/dp, d^j (p; p)/dp G V2, from (17) and the potential matrix elements, Qij(p), Hij(p), from (13), and the corresponding numerical values, deh/dp, /dp G V1, and Qh, Hh, the following estimates are valid [4]

d p d p

j(p)

< C3h2p,

dfa (p;p) 9 $

d p

d p

< c4hp+1,

(18)

\Qn(p) -Qij\

< c5h2p, \ Hij (p) -Hï\ < c6h2p

where h is the maximal grid .step, p is the order of finite elements, i, j are the numbers of the corresponding solutions, and the constants c3, c4, c5 and c6 do not depend on the step h.

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 calculated 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 Qa(1.8,10) =

0

{1.8(24)3(10)4(5)5(10)10} and the set of binding energies of the trimer Be3 calculated of the grid Qp = 4.1(20)7(10)10 with the eighth-order Lagrange elements (p = 8) are shown in Table 1, and the components of the trimer eigenfunctions (9) are shown in Fig. 4.

10

0,4

0,0 0,4

0,0 -0,4 0,4 0,0 -0,4

0,0

-0,4 0,2 0,0 -0,2

\ X]<1o,(p,E=-300.58)

^■o)(p,E=-260.48) "

/\ )f0)<P, E=-227.67)

xf4,<"(p£=-204.86)

\J"' xf'ol(P,E=-193.89) _

10

p,A

Figure 4. Components x)

(p,E) = Xj )(p) of gerade (g) and ungerade (w) bound

states of the trimer Be3 with total energy E 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 p

(A(p;p)-£j(p))<(p;p)= (-+p2V(psinp)-£j-(p)^ <(p;p)=0. (19) Using the new variable x' defined as p = x'/p, x' = p arcsin(x/p) we get

(--X2 + V(psin(xVp)) - <(X;p) = 0. (20)

In the argument of the potential we add and subtract xX and expand the potential in Taylor series in the vicinity of x',

V(psin(x'/p)) = V(x' + Ax')) = V(x') + £ ^dx^^¡xP"' (21)

k=1

where the small correction Ax' = p sin(x'/p) — x' ^ 1 is presented in the form of Taylor series

vk

Ax' = £ ( )

fc=i (2k + 1)! p2k

Then the Sturm-Liouville problem (19) is reduced to

i & A (fc)(X) (p)\ ± f , ^ n

+V (x) +g - h (x ; p) = °

^max

(<Mp)|h(p)) = J dx'(Mx'(xX;p) = hj.

xo

Here first terms V(k)(x; ) of the asymptotic expansion of V(p sin(x//p)) read as

(22)

V ("<x')=-T ^ ■

V"V)=x° ("^) /7 /

V (T)(x )=-(W 2 dTVixi+63xd2v (x )

V (x )= 4536° ^35x dx3 +63x dx2 +9 dx )'

T,(4), A x'9 / ,3d4V(x') nnn 2 d3V(x') nnn ,d2V(x') dV(x')\

V(4)(x')=—-( 175x 3-^+63°x 2-^+369x-^+15—.

v ; 54432°° V dx'4 dx3 dx2 dx

(23)

' 5443200 \1,5x '630x >369x

Note that fciy (k)( X) are the derivatives of the potential of the BVP (22) with respect to the parameter p-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 (x'min, x'max) in the finite interval X £ (x'min, x'max) to a generalized eigenvalue algebraic problem with respect to the eigenvalues A £ {Xj}j=1 and the corresponding eigenvectors c(p) = {(Ci(p))J}j=1 of the order L. 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 0(h2p) for the eigenvalues and 0(hp+1) « 10-8 for the eigenfunctions, h being the maximal grid step. L equals the number of grid nodes, which in the considered example was L = 393.

A(p)c(p) — c(p)BA(p) = 0, (24)

(c(p))TBc (p) = 1. (25)

Here the matrix A(p) is presented in the form of inverse power series

kmax « — 2k

A(p) = A(0) + V ^A(k). (26)

k=i k'

We choose as unperturbed the matrix operator A(^ corresponding to the differential one, — -^72 + V(x'), and numerically solve the corresponding algebraic problem

A(°)c(°) — X(°)B c(°) = 0, c(0)TB c(°) = 1. (27)

The solution A(p), c(p) is sought in the form of inverse power series

fcmax —2k fcmax —2k

A(p) = Af (p) + £ c(p) = c(°) + £ pF-c(fc). (28)

fc=i ' fc=i '

The substitution of Eq. (28) into Eq. (22) leads to the system of inhomogeneous algebraic equations for the corrections A(fc) and c(fc):

( A(°)c(k) — a(°)b c(k)) + (A(k)c(°) — X(k)B c(°)) = b{k)

hk) - £ (Bc(k—n)X(n) — ^c^ 6(i) - 0, (29)

n=1

A(fc) - —^, A(fc) - -2d.

Differentiating the normalization condition (25), we obtain additional conditions for the corrections c( fc):

C(°)TBC(fc) - F™ = — „ kl „c(k—n)TBc(n).

N 2 i n!(k — n)!

n=1 v '

Multiplying (29) by c(°)T and taking the zero values of the first term and the normalization condition into account, we get the formula for the corrections

\(k) of the

eigenvalues:

X(k) = c(°)TA(fc)c(°) — c(°)T b(k). (30)

The vector c(k) is calculated by solving the system of algebraic equations

Kc(k) -A(°)c(k) — A(°)Bc(k) = b(k) k j •

l)(k) = _ ___(A(n)c(k—n) - X(n)Bc(k—n))

n=1 n! (k — n)!v (31)

Y ,nkl „ c(k—n)TB c(n) =0, n!(k — n)!

n=°

where the latter equation is a result of differentiation of the normalization condition. Since A(°)

is an eigenvalue of (27), the matrix K 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 v(k) and w of the auxiliary inhomogeneous systems of algebraic equations

K v(k) = 6(fc), Kw = d (32)

with the non-degenerate matrix K and the right-hand sides b(k) and d

Kij, (i — s)(j — s) = 0,

Кц = {

% (i - s )0'- *0 = 0, Й-) = I "Г- :i = »■ d, = ( Кi =

J l 0, j = s, ' \ 0, J = s,

where s is the number of vector c(0) element having the maximal absolute value. Step k2. Evaluate the coefficient 7(k)

(k) P(k)

7(k) = — 7l(0)— , 7(k) = v(k)Tc(0), 72 = wTc(0). (33)

(Cs ) — 72)

( k )

Step k3. Evaluate the vector cj

c(k) = j — 7(k)-j, . = (34)

j I 7(k), 3 = s.

The execution of the above procedure yields the required asymptotic expansion for the eigenvalues

/ \ kmax kmax _2 k

j( p) (0) -2 k ( k) (0) p- ( k)

ef00+ ^f^f = Af(f) + £ (35)

2 3 \r> 1 ^ j 1 м j

И k=1 k=1

and the corresponding expansion of the eigenfunctions in the nodes xi of the given grid

kmax p-2k

Фз(x'i;p) = (c(0)b- + £ ^-(c(%j.

к=1 !

Following Eq. (14) we present it in the form of a piecewise-polynomial function

P VII

1 1 ГУ»' _ ГУ»'

iii v I "i—г p''

ф(x e [xpk,xp(k+i)];p) = ф(xpk+p;p) ц x , . (36)

P=0 p"=0,p"=p' P' P"

We calculate the discrete spectrum solutions фj(p; p) of the problem (19) on the grid pi = xi/p, related to the solutions фj (x'i; p) of the problem (20) as

ф (Pi = xi/p;p) = л/рфз (x'i ;p). The corresponding piecewise-polynomial functions are

ф(Pi e [Ppk^(k+i^p^^ (pi=xi/p e [x'pk/p,xp(k+i)/p];p) =

P P P P

" " A . -Г-Г tpp-x'w'

= (ppk+P' ;f) П (п,-Рп '„=(^pk+p' ;f) П X^v

P'=0 P'=0 ^ P=0 P'=0 XP' XP'

p"=p p"=p

; p) n

p=0

_ ry t

pi ^ pi '

p' ' = 0 p' '=p'

p ( -2k \ p

=VP E (C(0))«+£ V(C(k)M n

p=0 \ k=1 ' p '=[

p' '=0 p ' =p

^p' ^p '

V~ ■ (37)

The nodes xXp, 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 Qv[1.8/p, 10/p] = {1.8/p(24)3/p(10)4/p(5)5/p(10)10/p} by a Lagrange polynomial of the order p = 8, which provides the same relative accuracy of 0(hp+1). The obtained approximations of functions are polynomials of p, explicitly depending on p. 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 <j (<p; p) and/or

their parametric derivatives in each subinterval of the length h = pp(i+1) — pp of the nonuniform finite element mesh and summing the obtained results, we get the required asymptotic expansion of matrix elements (13) for cluster states, i,j = 1,...,n0:

qh(p) = £

k=1

Ç-1) p2k—1

Hij(p)

kmax TT k=1

(2 k)

p

2 k

(38)

( k )

For the beryllium trimer in collinear configuration, the coefficients ej of the expansion

(35) and the coefficients HJ(k) of the expansion (38) at i = j are presented in Table 2, and the first coefficients of the expansions (38) are

cS) =

/ 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

V 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.63427 -32.2355 0 52.0431 -53.8134

-1.58999 12.5531 -52.0431 0 105.275

V 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

6.68774 -91.6918 0 68.9065 595.104

14.8770 -27.6600 -68.9065 0 -523.035

V -128.505 356.227 -595.104 523.035 0 /

Coefficients e f] and ii(fe) of the expansions (35) and (38)

Table 2

3 1 2 3 4 5

e?> -193.06601252 -119.39267226 -63.338854932 -24.904560537 -4.0897890782

4X) -127.73059638 -317.30568182 -408.25519644 -385.66879485 -223.26092504

£f) -215.85831875 -672.02680456 -1198.7991074 -1892.8857386 -3079.0060449

(3) 4 -667.46086524 -2093.0917304 -3656.2842793 -5140.5360777 -3803.8470973

(4) -2590.0140097 -8317.9479477 -15062.850569 -23930.063302 -63419.196580

(6) 4 -11287.527768 -37281.527942 -69636.528276 -107755.44601 164277.292590

H(2) n 127.98059638 317.55568182 408.50519643 385.91879484 223.51081140

H{i) J] 414.47883932 1204.0042402 2009.0880648 3041.7739524 5091.1634806

i7(6) n 1912.4702415 5431.4025830 8482.5056376 10047.915436 -664.09784640

iï(8) n 9895.5812579 28061.019351 44415.159972 62195.582167 189759.82059

i7(10) n 54127.043138 154284.86093 246059.16425 311539.17441 152663.97768

eV+HW ] 1 n 0.2499999984 0.2499999944 0.2499999948 0.2499999953 0.2498863674

Table 3

Convergence of the expansion (35) to the numerical values (NUM) at p = 20

CO ¡=

3

o

a fo

s

I

-3

3

s

a

o S CO

3 1 2 3 4 5

4°V -193.06601252 -119.39267226 -63.338854932 -24.904560537 -4.0897890782

+4XV -193.38533901 -120.18593646 -64.359492923 -25.868732525 -4.6479413908

+e?)/p6 -193.38668813 -120.19013663 -64.366985417 -25.880563060 -4.6671851786

+sf/PS -193.38669856 -120.19016933 -64.367042547 -25.880643381 -4.6672446137

+sf/P10 -193.38669866 -120.19016966 -64.367043135 -25.880644316 -4.6672470910

+e?)/p12 -193.38669866 -120.19016966 -64.367043142 -25.880644327 -4.6672470749

(NUM) <pa = tt/2 (NUM) <pa = tt/6 -193.38669866 -193.38669866 -120.19016966 -120.19016966 -64.367043142 -64.367043142 -25.880644327 -25.880644327 -4.6672471622 -4.6672471630

«S

iZ o

13

H,

(2) _

H,

(4) _

H,

(6)

/ 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

V -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

-526.641 -420.171 2009.09 -947.285 -464.551

464.604 -812.366 -947.285 3041.77 -2605.05

V -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

-2112.77 -1689.40 8482.51 -2016.27 -12077.1

815.251 -3295.46 -2016.27 10047.9 12431.9

V 4153.61 -646.435 -12077.1 12431.9 -664.074 /

The convergence of the expansion (35) to the numerical values Table 3. One can see that the first six terms of the asymptotic accuracy of 11 significant digits.

(NUM) is presented in expansion provide the

Remark 1. The effective potentials £j<1)+#j<jl)=(1/4)p-2 (presented, e.g., in the last line of Table 2) lead to the asymptotic cluster fundamental solutions of the ODEs (11), determined by the Bessel functions J\/2(^/—E+£jp), j=1,..., no, while for pseu-

dostates the asymptotic fundamental solutions of the ODEs are Jm(\^Ep) with the integer m=(j-no)=1, 2,....

Algorithm for calculating the parametric angular functions of pseudostates and the effective potentials

The eigenfunctions of pseudostates Sj(p) > 0, (j — n0) = 1, 2,..., are localised beyond the potential well. Then the (n0 — 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 fa(ip; p) corresponding to the eigenvalues £j(p) ^ 0 of the gerade (g) pseudostates j = n +1 = 6,..., 10 at p = 100

From this fact the estimate of the eigenvalues for pseudostates Sj (p) ~ (j — n0)2 follows, namely, the eigenvalues of the corresponding BVPs (22) in new variable x', £j = £j (p)/p2,

will be a small quantity (see Fig. 1b). So, the solution en(p) (e n — £n (p)/p2) of the derived equation is sought in the form of a power series at n/p < 1

fcmax (fc)

£n(p) = n2+ V .

fc= 1

P

(39)

Then the numerical values of the function B(p^; p) — B(x'i) and its derivative B' (pi; p) — pB' (x'i) on the specified grid = {<p1 = p0,.. .,<pi — Xi/p,..., — Pe} in the polar system of coordinates are determined via the values of the function B(xi) and its derivative B'(xi) on the grid Qxi, {xi — XQ , ... , X^ , ... , X — Xg }, found with the same accuracy accepted in the FEM scheme chosen above in the form of the power series of the small parameter en:

B, (xi)=Br +

B?\e

(1) n 5

Jk)

, in

)

fc=1

b'0 (^(0)+i; ^^

b(k)(e

(1), n 5

Jk)

• . , in

)

(40)

fc=i

using the Runge-Kutta method, in which the all terms contain the power kmax + 1 of

1/p and the higher ones are neglected. The expansion coefficients B= B^^x'j) and

bf^ = b(-k^(xJi), calculated at the grid nodes x'i for the BVP (22) with the potential V(x') defined in Eq. (5) are presented in Fig. 6.

Figure 6. Expansion coefficients B(k) (e^,.. ., e^ ) and b( fe)(el1),.. . ,£in)), k=0, 2, 3, 4, 5,

1 with e^ given by Eq. (44) calculated at the nodes x\ of the grid Qx

fc

fc

P

P

n

One can see that in the vicinity of the potential well the corrections to the eigenfunctions are small, and at large x' they become essential. The coefficient 6(0), the derivative of the wave function with en — 0, exponentially tends to a constant for x > 5.5. From these observations the condition for choosing x£ 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 x'e — 10.

The interval p0 < p < n/2 is divided into two subintervals by the point — x£/p: (p0 < p < and n/2 > p > <ps. In the calculations the point xe was chosen from the condition |y(x > xe)| < e, where e is a preassigned number, and the left-hand boundary of the interval p0 — 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,

po = x0/p > 0. The eigenfunctions fa (<p; p) are calculated in the form

f Aj (P)Bj (v; P), Vo < <p < Pe,

^(V] p)=\ Cj(p)Jf {cs°ns} - ./2)), < k/2,

ff/2

2 dp(M<-p; p))2 = 1. Jf0

(41)

Here Aj (p) and Cj(p) are the normalisation factors, and B(<p; p) is determined from the numerical solution B(x) in Cartesian coordinates using the transformation <p = x/p. From the continuity of the eigenfunctions and their derivatives,

- 0; p) = + 0; p), d^t n ^ d^ (42)

d^ ^ - 0; P) = ^^ + 0; P),

we get the equation for the eigenvalue £n(p):

j tan(^en(p)(p£-f )) even _Vf"M=o

)

[ -COt^£n(p)(p£-f )) Odd n j R

R= B'n(Ve; P) = pB'n(X£) Bn(p£; p) Bn(x£) .

(43)

Substitute (40) into (43), and then substitute (39) into the resulting equation. Expanding both sides of the equation in inverse powers of p, we arrive at the system of linear equations, from which the expansion coefficients s^, and then the coefficients An(p) and Cn(p) are determined.

Since the values of the function Bn(p; p) and its derivative B'n(<p; p) 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 <ps(p) < <p < k/2, and the explicit dependence of its values upon the parameter p on the grid For the considered potential (5) we get the asymptotes of the g and u potential curves at n = j — n0 for n/p < 1:

kmax (k)

- £n

~7k

£n(p)=n2 + ^^ ; for ^-parity at odd n, or «-parity at even n. (44)

k=1

The first terms of these expansions are expressed as

41)=4.4862257n2, en2)=15.094666n2, en3)=-88.164324n2+0.2988408n4, en4)=-770.50344n2+3.3516683n4, e n5)=—13803.854n2+162.14850n4-0.2922734n6, e n6)=—68926.763n2+1584.0782n4-4.4329308n6.

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 p= 50 and the numerical results (NUM)

gerade ((/-parity)

n 1 3 5 7

n2 1.000000 9.00000 25.00000 49.00000

+£ (n)/p 1.089724 9.80752 27.24311 53.39650

+£ ^/p2 1.095762 9.86186 27.39405 53.69235

+£ ^/p3 1.095059 9.85570 27.37792 53.66353

+e i4)/p4 1.094936 9.85464 27.37517 53.65878

+£ ^/p5 1.094893 9.85428 27.37437 53.65775

+£ ^/p6 1.094888 9.85425 27.37432 53.65774

(NUM) 1.094887 9.85424 27.37431 53.65776

ungerade ( «-parity)

n 2 4 6 8

n2 4.000000 16.00000 36.000000 64.0000

+£ i1)/p 4.358898 17.43559 39.230082 69.7423

+£ ^/p2 4.383049 17.53219 39.447445 70.1287

+£ p3 4.380266 17.52152 39.425152 70.0934

+£ l4)/p4 4.379781 17.51968 39.421409 70.0877

+£ ^/p5 4.379613 17.51911 39.420448 70.0868

+£ {r6)/p6 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 p = 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 m = i — n0 and n2 = j — n0 of the same (g or u) parity at n0 = 5, i,j = n0 + 1,... for n/p < 1 have the form:

Qn1U2 (p)-

"-max

Qnin2

0fc+2

fc=0

Hnin2 (p) —

^max

^^ HniU2

k=0

0fc+4

(45)

First terms of these expansions read as

Q&—2.2431128^^2—2|y, Q&, —5.0315554^-^,

Q(4) —188.67822 -1-^oo,O-I-2(-2"-2)

+0.22413-

(-i--2) ' 22

H&,—10.06311 -1(-22(--+;2), Hi%2 —45.14538^iffl-+69.93509-1-2, ( n1 -2) ( n1 -2)

H,(°L—-1642.273-«+8,44Ä, +476,762„,„,

i - - 2 )

33

Hnim =0.6289444+2.0691442n?, ff^U =2.8215867+79.217746n2, = - 102.642068 + 138.328874n? + 0.826991n4.

Using Eq. (35) and Eq. (41), we get the asymptotic expansions for Qin(p) and Hin(p) between the cluster states i=1,... ,n0 and pseudostates n—n0=1, 2,..., for n/p<1

km ax ^(k+5/2) km ax ^ (k + 7/2)

Qin(p)= E , Hin(p)= E -k+T/-, (46)

k=0 F k=0 F

The first terms of these expansions read as

Ql5n/2)=-0.428911n, Ql7/2)=—1.443145n, Qi9/2)=-25.01600n+0.183791n3,

Q25/2)= 1.273900n, Q2n2)=+4.286254n, q29/2)= 75.66328n-0.546657n3,

Q35/2)=-2.497511n, q37/2)=-8.403301n, Q39/2)=-152.3919n+1.075421n3,

Q45/2)= 3.668834n, Q47/2)=+12.34441n, q49/2)= 231.9007n-1.599005n3,

Q55/2)=-4.098659n, Q5?n/2)=-13.79066n, Q59/2)=-247.9253n+2.018873n3,

^ln/2)= 22.53006n, ^ln/2)=+77.24936n, fijni/2)=1551.8678n-9.835892n3,

H2 n/2)=-26.57543n, H2 n/2)=-93.70381n, H{2 ni/2)=-2203.032n+11.89657n3,

^3 n/2)=-12.19890n, #3n/2)=-32.64199n, #3ni/2)=262.11512n+4.400493n3,

^4 n/2)= 85.06821n, H4 n/2)=+273.8820n, H4n/2)=4387.8346n-35.91925n3,

^5 n/2)=-120.4823n, ^5 n/2)=-391.5948n, H5( ni/2)=-8001.428n+54.83111n3.

5. Conclusion

The model for beryllium trimer in collinear configuration is formulated as a 2D boundary-value problem for the Schrodinger 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 Wij(p) in inverse powers of p. It can be used for calculation of the asymptotes of the fundamental solutions of the system of second-order ODEs at large values of p [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, A. GoZdZ, 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. GoZdZ, 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. Spirko, K. Szalewicz, On the Elusive Twelfth Vibrational State of Beryllium Dimer, Science 326 (2009) 1382-1384.

УДК 519.632.4

Алгоритмы решения краевых задач для атомных тримеров в коллинеарной конфигурации методом Канторовича

А. А. Гусев*, О. Чулуунбаатар*+, С. И. Виницкий**,

В. Л. Дербов§

* Объединённый институт ядерных исследований, г. Дубна, Россия ^ Институт математики, Монгольский государственный университет, Улан-Батор,

Монголия

* Российский университет дружбы народов, Москва, Россия S Саратовский государственный университет, г. Саратов, Россия

Модель атомных тримеров с парными молекулярными взаимодействиями в коллинеарной конфигурации формулируется в виде двумерной краевой задачи в якобиевских полярных координатах. Последняя сводится методом Канторовича к одномерной краевой задаче для системы ОДУ второго порядка, используя разложение искомого решения по угловым базисным функциям, зависящим от гиперрадиуса, как от параметра. Представлены алгоритмы решения параметрической краевой задачи методом конечных элементов и вычисления асимптотических разложений параметрических угловых базисных функций и эффективных потенциалов системы ОДУ при больших значениях параметра. Эффективность алгоритмов подтверждается сравнением асимптотических решений параметрической задачи на собственные значения и эффективных потенциалов с их численными значениями, полученных методом конечных элементов при больших значениях параметра. Применимость алгоритмов демонстрируются на примере расчетов асимптотических разложений решений параметрической краевой задачи и эффективных потенциалов, и собственных значений энергии связи тримера бериллия в коллинеарной конфигурации.

Ключевые слова: краевые задачи, метод Канторовича, системы ОДУ второго порядка, метод конечных элементов

Литература

1. Efimov V. Few-Body Physics: Giant Trimers True to Scale // Nature Phys. — 2009. — Vol. 5. — Pp. 533-534.

2. Observation of an Efimov Spectrum in an Atomic System / M. Zaccanti, B. Deissler, C. D'Errico, M. Fattori, Jona-Lasinio, S. M. Muller, G. Roati, M. Inguscio, M. G. // Nature Phys. — 2009. — Vol. 5. — Pp. 586-591.

3. Voigtsberger J. et al. Imaging the Structure of the Trimer Systems 4He3 and 3He4He2 // Nature Commun. — 2014. — Vol. 5. — P. 5765.

4. 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 / O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashke-vich // Computer Physics Communications. — 2009. — Vol. 180. — Pp. 1358-1375.

5. 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 / A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, A. G. Abrashkevich // Comput. Phys. Commun. — 2014. — Vol. 185. — Pp. 26362654.

6. Описание программы вычисления собственных значений и собственных функций и их первых производных по параметру для параметрической самосопряжённой системы эллиптических дифференциальных уравнений / А. А. Гусев, О. Чулуун-баатар, С. И. Виницкий, А. Г. Абрашкевич // Вестник РУДН, серия «Математика. Информатика. Физика». — 2014. — № 2. — С. 336-341.

7. KANTBP: A Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adiabatic Approach / O. Chu-luunbaatar, A. Gusev, A. G. Abrashkevich, A. Amaya-Tapia, M. Kaschiev, S. Larsen, S. Vinitsky // Comput. Phys. Commun. — 2007. — Vol. 177. — Pp. 649-675.

8. 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 / O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich // Comput. Phys. Commun. — 2008. — Vol. 179. — Pp. 685-693.

9. Gusev A. A., Hai L. L., Chuluunbaatar O, Vinitsky S. I. Program KANTBP 4M for Solving Boundary-Value Problems for Systems of Ordinary Differential Equations of the Second Order. — http://wwwinfo.jinr.ru/programs/jinrlib/kantbp4m.

10. Symbolic Numerical Algorithm for Solving Quantum Tunneling Problem of a Diatomic Molecule Through Repulsive Barriers / S. I. Vinitsky, A. A. Gusev, O. Chuluunbaatar, L. L. Hai, V. L. Derbov, P. M. Krassovitskiy, A. Gozdz // Lect. Notes Comp. Sci. — 2014. — Vol. 8660. — Pp. 472-490.

11. Krassovitskiy P., Pen'kov F. Contribution of Resonance Tunneling of Molecule to Physical Observables // J. Phys. B. — 2014. — Vol. 47. — P. 225210.

12. Wang J., Wang G, Zhao J. Density Functional Study of Beryllium Clusters with Gradient Correction // J. Phys. Cond. Matt. — 2001. — Vol. 13. — Pp. L753-L758.

13. Стренг Г., Фикс Г. Теория метода конечных элементов . — Москва: Мир, 1977.

14. Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов. — Москва: Стройиздат, 1982.

15. Березин И. С., Жидков Н. П. Методы вычислений. — М.: Физматлит, 1962. — Т. 1.

16. Гусев А. А., Хаи Л. Л. Алгоритм решения двумерной краевой задачи для модели квантового туннелирования двухатомной молекулы через отталкивающие барьеры // Вестник РУДН, серия «Математика. Информатика. Физика». — 2015. — № 1. — С. 15-36.

17. Pijper E., Fasolino A. Quantum Surface Diffusion of Vibrationally Excited Molecular Dimers // J. Chem. Phys. — 2007. — Vol. 126. — P. 014708.

18. Модели резонансного туннелирования составных систем через отталкивающие барьеры / А. А. Гусев, А. Гоздз, В. Л. Дербов, С. И. Виницкий, О. Чулуунбаатар, П. М. Красовицкий // Новости ОИЯИ. — 2014. — № 1. — С. 22-26.

19. Mitin A. V. Ab Initio Calculations of Weakly Bonded He2 and Be2 Molecules by MRCI Method with Pseudo-Natural Molecular Orbitals // Int. J. Quantum Chem. — 2011. — Vol. 111. — Pp. 2560-2567.

20. Merritt J. M., Bondybey V. E., Heaven M. C. Beryllium Dimer-Caught in the Act of Bonding // Science. — 2009. — Vol. 324. — Pp. 1548-1551.

21. Patkowski K., Spirko V., Szalewicz K. On the Elusive Twelfth Vibrational State of Beryllium Dimer // Science. — 2009. — Vol. 326. — Pp. 1382-1384.

© Gusev A. A., Chuluunbaatar O., Vinitsky S.I., Derbov V. L., 2016