CC BY
21
UDC 517.958:530.145.6
Algorithm for Computing Wave Functions, Reflection and Transmission Matrices of the Multichannel Scattering Problem in the Adiabatic Representation using the Finite Element
Method
A. A. Gusev
Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow region, 141980, Russia
In adiabatic representation the multichannel scattering problem for a multidimensional Schrodinger equation is reduced to the boundary value problem (BVP) for a system of coupled self-adjoined second-order ordinary differential equations on a finite interval with homogeneous boundary conditions of the third type at the left and right boundary points in the framework of the Kantorovich method using adiabatic basis of surface functions depending on longitudinal variable as a parameter. The homogeneous third-type boundary conditions for the desirable wave functions of the BVP are formulated using the known set of linear independent regular and irregular asymptotic solutions in the open channels of the reduced multichannel scattering problem on an axis which involve the desirable reflection and transmission amplitude matrices, and the set of linear independent regular asymptotic solutions in the closed channels. The economical and stable algorithm for numerical calculation with given accuracy of reflection and transmission matrices, and the corresponding wave functions of the multichannel scattering problem for the system of equations containing potential matrix elements and first-derivative coupling terms is proposed using high-order accuracy approximations of the finite element method (FEM). The efficiency of the proposed algorithm is demonstrated by solving of the two-dimensional quantum transmittance problem for a pair of coupled particles with oscillator interaction potentials penetrating through repulsive Coulomb-type potentials and scattering problem of electron in a Coulomb field of proton and in the homogeneous magnetic field in the framework of the Kantorovich and Galerkin-type methods and studying their convergence.
Key words and phrases: multichannel scattering problem, reflection and transmission amplitude matrices, multidimensional Schrodinger equation, adiabatic representation, the Kantorovich method, boundary value problem, a system of coupled self-adjoined second-order ordinary differential equations, the finite element method.
1. Introduction
In the adiabatic representation [1] the multichannel scattering problem for a multidimensional Schrodinger equation describing three-dimensional tunneling of a diatomic molecule incident upon a potential barrier [2], fission model of collision of heavy ions [3] or Coulomb scattering with transversal confinement produced by uniform magnetic field or channeling ions in crystal [4-6] is reduced by separating the longitudinal coordinate, labeled as z, from the transversal variables to the boundary value problem (BVP) for a system of self-adjoined second-order ordinary differential equations containing the potential matrix elements and first-derivative coupling terms. Such reduction of the problem is performed in the framework of the Kantorovich method [7] using basis of surface eigenfunctions by transversal variables of auxiliary BVP depended on the longitudinal variable as a parameter [8]. In order to guarantee high-order accuracy of numerical solution the BVP the relevant potential matrix elements should be evaluated with the same level of accuracy as the approximate solutions. The corresponding algorithms of numerical solution with given accuracy for the parametric two-dimensional boundary-value problem and calculation of the solution derivative with respect to the parameter and the matrix elements using
Received 18th February, 2014.
The author thanks Drs. O. Chuluunbaatar, V.L. Derbov, L.A. Sevastyanov and S.I. Vinitsky for collaboration. The work was supported partially by grants 14-01-00420 and 13-01-00668 RFBR..
the finite-element method (FEM) [9-11] have been elaborated [8]. The variable coefficients of the ordinary differential equations and the corresponding solutions can have a long-range asymptotic behavior by inverse powers of the independent variable. In this case to reduce the problem to the finite interval of integration the new algorithm for evaluation of asymptotic expansions of desirable solutions, in series by inverse powers of the independent variable based on using an appropriate etalon equation has been elaborated [12]. For problems similar to the penetration of composite system through the long-range repulsive barrier potentials, solving the auxiliary BVP depended on the center mass variable as a parameter is rather complicate problem [13]. In this case we will include the barrier potentials after averaging over basis of surface parametric eigenfunctions of composite system as additional potential matrix into the systems of coupled self-adjoined differential equations with respect to center mass variable derived in the framework of the Kantorovich method. Preliminary we studied this problem using averaging barrier potentials over the basis of surface eigenfunctions of the Galerkin-type method [14,15]. Thus, to solve the above problems we need to formulate an appropriate BVP and develop method, algorithms and software together with benchmark calculations which will reveal specific features of realization of the Kantorovich and Galerkin-type methods and their combinations.
The purpose of this paper is to present a suitable formulation of the multichannel scattering problem for a multidimensional Schrodinger equation based on the Kan-torovich and Galerkin-type methods implemented as an economical and stable algorithm based on high-order accuracy approximations of the boundary value problem using the FEM. The third-type boundary conditions are formulated for the considered scattering problem with respect to the desirable wave functions, reflection and transmission matrices by using the known set of linear independent asymptotic regular and irregular solutions in the open channels, and the set of linear independent regular asymptotic solutions in the closed channels, respectively. An essential part of the resulting algorithm consists in economical formulations of nonhomogenous algebraic problems using the matrices of logarithmic derivatives of asymptotic regular solution at one boundary point and the matrices of logarithmic derivatives of the solution calculated in the finite interval at another boundary point used to determine the reflection matrix. As a benchmark calculation, the algorithm implemented in the form program KANTBP 3.0 [16] is applied to computing the transmission coefficient for the 2D-model of a pair of particles, coupled by the oscillator interaction potential, describing the penetration through symmetric or nonsymmetric barriers, as well as the long-range repulsive truncated Coulomb [13] and Coulomb like barriers [12], and the Coulomb scattering with the transversal confinement oscillator potential produced by uniform magnetic field [4].
The paper is organized as follows. In Section 2 we give a brief overview of the problem. In section 3 the nonhomogeneous algebraic problem using FEM is formulated. In section 4 the description of the auxiliary algorithm for the calculation of linear independent asymptotic regular and irregular solutions in open channels, and linear independent regular asymptotic solutions in closed channels is given. In section 5 the benchmark calculations of penetration coefficient and analysis of convergence within the framework of Kantorovich and Galerkin-type methods are presented. In conclusion a brief summary is given and further applications are described.
A wide class of quantum-mechanical problems are reduced to the solution of multidimensional Schrodinger equation for the wave function ^(r, x):
2. Statement of the Problem
Here L(x; z) = — AX + U(z, x) is a self-adjoin elliptic differential operator with partial derivatives in a finite region X c Rd , x = {xj G X is a set of independent (fast) variables, z G (zm-m, zmax) G B c R1 is independent (slow) variable, X = B®X c Rd
is a finite region of configuration space Rd ; E is a spectral parameter, corresponding to the energy of a quantum system. It is assumed that the functions f1(z) > 0, f2(z) > 0, f3(z) > 0, dzf2(z), U(z,x) and dzU(z,x) are continuous and bounded for all ( , x) G X. It is also assumed that the self-adjoin operator L( x; ) for each value of z G (z1 = zmin, z2 = zmax) G B c R1 has only a discrete real spectrum e(z).
The solutions z,x) G L2(X) of Eq. (1) are subject to the boundary conditions of the third kind
Vi 9^, x) — \i V(z,x) = 0, z = Zi, x G dX UX,
dz (2) dz,x) , , , T , , _ - r ,
a—--b(z)w(z ,x) = 0, x GdX, z G [zmin, zmax\,
where a are real constants; Xi = Xi(zi) are real functions depending on zi;
V2 + ^2 = 0; b(z) and dzb(z) are continuous and bounded functions; a2 + b2(z) = 0; n is unit vector normal to the bounds dX of domain X.
In Kantorovich method [7\ the partial wave function ^i(r, x) sought in the form of an expansion over the set of the one-parameter basic functions{^j (x;z)}^=1 G Tz ~
L2(X):
N
*i( z,x) = (x;z)xf (z). (3)
j=i
In the expansion (3) the functions {%ji)( z)}1N=1 are unknown. The basis functions {^j(x;z)}1N=1 are defined as parametric solutions of the eigenvalue problem
d n
L(x;z)^j (x; z) = £j (z)ipj (x;z), - b(z)ipj(x;z) = 0, x G dX, z G [zmin, ¿max].
Here e1(z) < ■ ■ ■ < £N(z) < ■ ■ ■ G e(z) is the set of the desired N real-valued eigenvalues, arranged in the ascending order. The basis functions {^j(x; z)}1N=1 form
an orthonormal basis for a set of variables x = {xj }d=i G X for each value of G (zmin, zmax) G B, which is regarded as a parameter. The basis functions satisfy the orthogonality and normalization conditions
(^i(x;z)tpj(x;z)) = / t^i(x;z)i^j(x;z)dxd'-1 = (5)
\ / x Jx
where i is the Kronecker symbol.
Substituting the expansion (3) in Eq. (1) and averaging over the orthogonal basis (4), (5), the multidimensional Schrodinger equation is reduced to a finite set of N ordinary second-order differential equations on the finite interval [zmin, zmax\ for the
partial solution z) = (xl^O2),..., Xn\z))
/ i d d (D - 2£ I) ^^ S(-I^ dAW^ + VW +
+i|Q<4+m1^-2*')№o. (6)
Here I, V(z) and Q(z) are the unit, symmetric and antisymmetric N x N matrices:
Vij (z)
ei{z) + £j( z) + f 2(z) n (
2f3(z) ^ + fi(z)Hi (Z), = ^ '
Hij (z) = Hji(z) =
' d^i(x;z)
Qij(z) = -Qji(z) = - ( l(x; z)
dz
di^j (x; z)
dz dij (x; z) d
(7)
The basis functions {l (x;z)}j^=1 form an orthonormal basis, and on the bound-
aries of the interval z G [zmin, zmax] the following conditions are satisfied:
(
N
(i) d^x) -A« ,x)) =0, * = *,
from which the uniform matrix boundary conditions of the third type follow
№ (id - Q(*)) x«(z) - A(f)X(l)(r) = 0, * =
(8)
(9)
We assume that V( z) and Q( z) matrices have the following asymptotic behaviour at large z — z± ^
Vij( Z±)=\ €j + ,utJ
( 2 7 ±\ i A 1 ,±)
U + 2^ % + E It-, ^±) = E V, (10)
V *± / l=2 Z± l = 1 z±
where e1 < ... < e n are the threshold energy values.
In the present work a scattering problem is solved using the boundary conditions at Z — Zm[n and Z — 2maX:
d§(z)
dz
d§(z)
d
= ¿max)$( Zmax), (11)
where z) is an unknown N xN matrix function, $(z) — (z)}Ij=L1 is the required N x N0 matrix solution and N0 is the number of open channels, N0 — max2^^ej j ^ N.
From this we obtain the quadratic functional (similar to Eq. (23) in [17] and Eq. (5) in [18])
E($,E, Zmin, Zmax ) = / (z) (D - 2E I) $(z)dz = n($, E, Zmin, *max)-
- f 2 (^max)^T(¿max)G(Zmax)$(¿max) + ¡2 (¿min)$T(¿min)G(Zmin)$(¿min), (12)
where n($,E, zmin, zmax) is the symmetric functional
d$T(z) d§(z)
П($, E, ^min ^max) J
M z)-
d d
+ fi(z)GT (z)Y(z)*(z) +
X
z=z
z=z
max
Z
Z
+h(z)*T( z)Q(z)- h(z)^^Q(*)*(*) - fi(z)2E$T(Z)G(Z)
dz
dz
dz, (13)
and G(z) = V,(z) — Q( z) is a N x N matrix function which should be symmetric according to the conventional R-matrix theory [19].
2.1. The Physical Scattering Asymptotic Forms of Solutions in Longitudinal Coordinates and the Scattering Matrix
The matrix solution ( z) = $( z) describing the incidence of the particle and its scattering, which has the asymptotic form "incident wave + outgoing waves" (see Fig. 1a), is
$v (z ^ ±ro) = <
X(+)(z)Tv, z> 0,
X(+)(z) + X(-)( z)Rv, z< 0,
X(-)(4 + X(+)(z)Rv, z> 0,
X(-)(z)Tv, z < 0,
= ,
= ,
(14)
where R^ and T are the reflection and transmission Na x Na matrices, v and v denote the initial direction of the particle motion along the z axis. Here the leading term of the asymptotic rectangular matrix functions X(±)(z) has the form [20,21]
(*) ^ Ï-X'2 exp ((p,z - Z ln(2pj|z|))) %
Pj = yj2 E - €j i = 1,...,N, j = 1,...,N0,
(15)
where Zj = Z+ at z > 0 and Zj = Z- at z < 0.
(z "+œ)
(z "+œ)
(z "+œ)
(z "+œ)
X(+)( Z) X(!)(z) X(+)(z)Tt_ X(-)(z)Rt_ X(+)( z)R# X(!)(z)T#
z > 0
X(!)(Z)R" z < 0 X(+)(Z)T" X(!)(Z)T# z > 0 Z < 0 X(+)( z)R# z > 0 z < 0 X(+)( Z) X(!)(Z) z > 0 z < 0
(a)
(b)
Figure 1. Schematic diagrams of the continuum spectrum waves having the asymptotic form: (a) "incident wave + outgoing waves", (b) "incident waves +
outgoing wave"
The matrix solution $v ( z,E) is normalized by the condition
CO
y '( s,E')*v(z, E)h(z)dz = 2ttö(E' - E)^Io
(16)
where I oo is the unit Na x Na matrix. Let us rewrite Eq. (14) in the matrix form at ^ and z- ^ —to as
0 X
0 X
(Z-)J \xM(z-) 0 J ' \X(-)(z-)
(-0( +(
S, (17)
)
0
where the scattering matrix S
(R+ T^\
VT^ R^J
(18)
is composed of the reflection and transmission matrices.
It should be noted that the functions X(±)(z) satisfy the relations
Wr(Q(z); X(^)( z), X(±)( z)) = ±2iI00, Wr(Q( z); X(±)( z), X(±)(z)) = 0, (19)
where Wr(»; a( z), b( z)) is the generalized Wronskian with the long derivative defined as
Wr(.; a(z), b(z)) = aT( z)
^ -b(^) -{^ -«) b(«). (20)
This Wronskian is used to estimate the desirable accuracy of the above expansion.
Let us prove that the .scattering matrix (18) is symmetric and unitary. Using Eqs. (14) and (19), we arrive at the following relations
+2«Tt T^, z > 0,
Wr(Q(z); *%(z), ( z)) =
Wr(Q(z); *t(z), ( z)) =
Wr(Q(z); *%(z), ( z)) =
Wr(Q(z); *t(z), ( z)) =
Wr(Q(z); &^(z), ( z)) =
Wr(Q(z); &^(z), ( z)) =
Wr(Q(z); *^(z), ( z)) =
+2i(I00 - R^R^), z< 0,
-2iTl_T^, z < 0,
-2i(Ioo - Rj_R^), *> 0,
+2îTlR^, z > 0, -2îRlT^, z < 0,
+2aRt_Ti, z > 0, -2îT|_RI, z < 0,
-2iTl, z > 0, -2ÎTI, z < 0,
(21)
+0,
> 0,
+2î(Rl - Ri), z< 0, -2î(RT - R^), z> 0,
+0, z < 0,
where the asterisk denotes complex conjugation. From here, we obtain the following properties of the reflection and transmission matrices:
T^Ti + R^R^ = I 00 = TtT^ + RtR^, T^R^ + R^T^ = 0 = Rj_Ti + Tj_Ri, TT = T^, RT = Ri, RT = R^.
This means that the scattering matrix (18) is symmetric and unitary.
(22)
Another type of the matrix solution $v ( z) = $( z) describing the incidence of the particle and its scattering, which has the inverse asymptotic form "incident waves +
S
outgoing wave" (see Fig. 1b), is
<!V (z — ±ro) = <
X(+)(z) + X(-)(z)Rt, 0,
X(+)(z)TI, z< 0,
X(-)(z)fV, z> 0,
X(-)(z) + X(+)(z)Rt, z< 0,
,
v —<-
(23)
Note, that the equality *± (z) = (z) should be valid from which we obtain
R—> — R^-, R-s—
R_
TV. Therefore we consider below only matrix-solution
2.2. Calculation of Matrices G(zmax) at v and G(zmjn) at v Suppose that the set of linear independent regular square-solutions (z) =
1 for the problem under consideration with the components Xreg(z) = (xT^z),..., Xng(Z))T is known at z > 0, v and at z < 0, v i.e.,
^(z) = X(+)(z), z> 0, $r!s(z) = X(-)(z), z < 0,
(j),
Xj±\z) — Xt±\z)
1,
,N, j — 1,
,Nn
If some channels are closed, we should use additional linear independent regular asymptotic functions at z > 0 and z < 0, respectively:
Xif (z) - q-1'2 exp ( ^ ( q3z + ln(2q3 |z|) I I ,
z±
(24)
qj — y/ej - 2 E, i — 1,...,N, j — Na + 1,...,N.
In this case the matrix G( z) at z — zmax > 0, v —■— or at z — zmin < 0, v —^ can be expressed via the known set of linear independent regular solutions $Veg( z)
d$>reg (z)
G( z) — n( z) - Q( z)— vdz ( )
(■Veg(z))-1 - Q(z).
(25)
The matrix G( z) at z = zmin, v and at z = zmax, v will be used in the next sections.
3. Formulation of the Algebraic Problem using the FEM
Computational schemes having high order of accuracy for the solution of the multichannel scattering problem (6)-(11) at a fixed value of energy E in open channels are derived from the variational functional (12), (13) on the basis of the FEM. The general idea of the FEM in one-dimensional space is to divide the interval [zmin, zmax] into many small domains referred as elements. The size of elements can be defined very freely so that the physical properties or qualitative behavior of the desired solutions can be taken into account.
The interval A = [zmin, zmax] is covered by a system of n subintervals Aj = [Zj-T, Zj] in such a way that A = (J"=
J"=1 Aj. In each subinterval Aj the nodes
zp — 3,r
h
Zj—1 \ T* , h1 - Z q Z q —
i-1,
0,p,
(26)
V
and the Lagrange elements {(fp r(¿0}p of the order p
ttr (*) = n
( * - zji)
(ZP _ ZP .)
¿=0,i=r (Ùl,r )
(27)
for the approximated function are determined. Using the Lagrange interpolation elements (z), we define the set of local functions Ni(z) as follows:
N[(z) =
<Pp,o(z), 0,
z G Ai,
ZG Ai, z e
z G Aj, z e A,
(Pp+l,o(z), ze Aj+1, 0, ^ AjUAj+i,
Vpn,p(z), z e An,
o, z<e An,
Vj,r (z), 0,
(pp,p(z), j
1 = 0,
l = r+p(j — 1), r = 1,p—1,
(28)
l = ЗP, j = 1,n — 1,
I = np.
The functions { Nf(z)}f= 0, L — np form a basis in the space of the piecewise polynomials of the (p+ 1)-th order. Now, each component of the vector functions x^to)(z) G H1(Qhz) is approximated by a finite sum of local functions Np(z)
1=0
(29)
i.e. the vector function x^t°)(z) has a generalized first-order partial derivative and belongs to the Sobolev space H1(Qhz) [9].
After substituting the expansion (29) into the variational functional (12), (13) the solution of the multichannel scattering problem (6)-(11) at a fixed value of energy E in open channels similar to [17] is reduced to the solution of the algebraic problem with respect to the matrix solution $h = ((x^ 1))h,..., (x^No))h):
Gp$h = (Ap — 2EBp)$h = (Mmax — Mmin) $
which subjects to the boundary conditions d$h(z)
dz
(G( z) + Q( z))$h( z), z = zn
z =
(30)
(31)
Here Ap and Bp are symmetric ( L N) x (L N) matrices, L is the number of the nodes of the finite element grid on the interval [zmin, zmax], Mmax and Mmin are ( L N) x (L N) matrices with zero elements except the right-lower and left-upper N x N matrices equal to ^(^max)G(^max) and /2(zmin)G(^min), respectively, G(z) — K(z) - Q(z) is the N x N matrix function, Ap is the stiffness matrix; Bp is the positive definite mass matrix; (x^o))h is the vector approximating the solution on the finite-element grid. The matrices Ap and Bp have the following form:
Ap
j=i
Bp = £ b
(32)
j=i
h
the local matrices ap- and b^ are calculated as
+i
1 4 d tpP ( z)dtpp (z)
K)^ =J{V^ *) at ) dz ) + ^^( ( z)^P'r(
4 , , vd ¥p,„ ( z)d¥P,r ( Z)
-1 J
2 f ( \ p ( \n f 2 ( / W / x Hi , /oq^
+ Y-fo^P'i--—(z),'pp'r(z) | 2 dr1, (33)
+1
(bp)^ — s^j fi(z)<p?>q(z)^Pr(z)hjdi], (34)
-i
z — Zj-1 + 0.5hj(1 + rj), q,r — 0,p.
In particular, the integrals (34) are evaluated using the Gaussian quadrature. Note, that in this approach maximum value of the half-band of the matrices Ap and Bp is small compared to their dimension and is not greater than N( + 1).
Let D(z) be a continuous and bounded positively defined operator in the space H1 with the energy norm, x(t)(z) G H2 are the exact solutions of Eqs. (6)-(11), and (x(t))h(z) G H1 are the corresponding numerical solutions. Then the following estimates are valid [9]:
||X(%) - (*W)h||o < chp' ,c> 0, (35)
where ||x(i )(
^)ll2 — Iz max dz(x(l^(z))TX(t)(z), h is the grid step, (p+1) is the order of finite elements, is the number of the corresponding eigensolution, and the constant does not depend on the step h.
3.1. Calculation of the Matrix Solution ( z)
First, we consider the numerical algorithm for calculating the matrix solution $h — . In this case Eq. (30) can be rewritten in the following form
(GP + MP ) r= )(*
(G + Mmin)^J - G^ ) UU
— Mw)[0 G( ^maxOQ^) . (36)
Here Gp and Mmin are the matrices determined in accordance with Eqs. (30), (31) and (25), while G(zmax) is the unknown matrix of the dimension N x N and and = £max) are the required matrix solutions of the dimension ( LN — N) x N0 and N x N0, respectively. From here, we obtain the explicit expressions
— —(Gta )-1Gi G(^max) — /2-1(w)(Gf — (Gta)-1Gt6). (37)
From Eqs. (31) and (37) we can obtain the relation between and its derivative
= *max)$t> Zmax) = G( Zmax) + Q( Zmax). (38)
Note, that the matrix G(zmax) is determined via the inverse of the submatrix G^ calculation that requires substantial computer resources. For evaluating Eq. (38) without such calculation of the inverse submatrix G^, let us consider the following auxiliary system of nonhomogeneous algebraic equations
( Gr Gt\ / FM ) ( 0 X
I G ba Qbb I I F*l_ J _ J Zm a I J ' ( )
The problem (39) is solved with the help of LDLT-factorization of the stiffness matrix and the back substitution method [22]. As the determinant of the matrix Gp + Mmin is nonzero, the above equation has a unique solution
F t = -(Gta)-1Gt6F t, F t = f 2(zmaGbl - Gt!(Gta)-1Gt6) . (40) Taking this into account, the required zmax) matrix is equal to
K( Zmax) = (Ft )-1 + Q( ¿max), (41)
and the required solution $t is calculated by the formulae (37) and (40)
= Ft (Ft)-1 , = X(-)(Zmax) + X(+)(Zmax)R^. (42)
Using the calculated zmax) from Eq. (41) and Eq. (14) for the asymptotic solution with the unknown reflection R^ and transmission T^ matrices, we obtain the following matrix equations for their calculation:
Y(+)( 7 = -Y(-)(r ) X(-)( ^ • = ( ^ • )
(^ max)R^ — (^ max), ■X (^min) T ^ — (^ mm),
Y(±)(, ) _ dX(±)(z)
A (^max)
dz
(43)
_T?( r )X(±)( z )
/v( ^ max)X ( ^ max)-
Note, that when some channels are closed, Y(±)(z) and X(-)( z) are rectangular N xNa matrices. Therefore, using the pseudoinverse matrices of Y(+)( z) and X(-)(z), we obtain the following expressions:
T \ — 1 T
Rt _ — i (Y(+) ( Zmax)) Y(+)( Zmax) (y(+) ( Zmax)) Y(—) ( Zmax),
V J i (44)
/ T \ — 1 T
Tt _( (X( —) ( Zmin)) X( —)( Zmin ) j (x(—)( Zmin)) $t( *min)-
Now we proceed to a brief description of the calculational scheme for the matrix solution $h _ $t. The required zmin) matrix is equal to
n(Zmin)_(Ft) —1 + Q(¿min) , (45)
and the required solution $t is calculated as
$t _ Ft (Ft) —1 $t, *t _ X(+)(Zmin)+ X( —)(Zmin)R^. (46)
Here $t = ^t(¿min) and $t are the matrix solutions of the dimension N x N0 and (LN — N) x Na. F^ and Ft are the matrices of dimension N xN and (LN — N) x N,
which are solutions of the auxiliary system of nonhomogeneous algebraic equations
i F^ \ ( G^ Gf \ ( F^ \ ( IN
(Gp - MmaxW = i ^ j i = - h( zmm)^0) . (47)
Finally, using the calculated zmin) from Eq. (45) and Eq. (14) we obtain the following matrix equations for calculating the reflection R^ and transmission T^ matrices
Y(-)(r . )R^ — _y(+)(r . ) X(+)( z )T^ — $h ( z )
(^min/R^ (^min)j (^max) T—V ^max/j
dX(±)(z)
Y(±)( r . ) —
A —> (^min)
dz
¿min)X( )( ¿min)-
The reflection R« and transmission T« matrices are evaluated using the pseudoinverse matrices of Y(-)(zmin) and X(+)(zmax)
i t \ — 1 T
R« = -( (y(— )(¿min)) Yi—)( ZminH (y(—)(¿min)) Y« (Zmin),
V J i (48)
/ T \ — 1 T
T« = fx (+)( ^ax)) X <+>(z maxM fa (+) ( ¿max)) ( ¿max).
4. Algorithm for Calculating the Asymptotic Forms of Regular and Irregular Solutions in the Longitudinal
Coordinate
We calculate the asymptotic solution of a set of N coupled ordinary differential equations (ODE) at large values of independent variable |z| ^ 1 for the particular case when f1(z) = f2(z) = zd— 1:
( 1 d zd—1 d + Vu(z) - 2E)Xw(z)
zd-1 dz~ dz
N
)
(Vij(z) + (z)d + -J-1 ^zd-1Qt,(z^ Xii'(z). (49)
We assume that the coefficients of Eqs. (49) can be represented in the general asymptotic form as
( f(1)\ V.W Q0
Vu (z) — 40) + Jr) * + £ Vir, Qu (*) — Z Q-. (50)
V J 1=2 1 = 1
Usually in the case d — 2,3,... the asymptotic solution is calculated in one interval zmax ^ z < In the case d — 1 the asymptotic solution is calculated in two
intervals —œ < z < zmin and zmax < z < where, in the general case, the
coefficients (1), V/P and Q{sl) are different for z > 0 and z < 0. Below we will consider
j lj lj
only the case > 0.
Step 1. We construct the solution of Eqs. (49) in the form:
X0%' (z) — ^ (z)Ri, (z) + ^i, (*)dR^), (51)
where ^/(z) and Jji/(z) are unknown functions, R^ (z) is a known function. We choose R^ (z) as the solutions of the auxiliary problem treated like an etalon equation
[z (k<i) = Z (k>k'max) = 0^.
K.
d i d , Z( ,
Ak)
,„-i dT ^ dr^-j- — Pi I* = <52>
k=i
Remark 1. If Z( = 0, then the solutions of the last equation are expressed in terms of the hypergeometric functions, exponential, trigonometric, Bessel, Coulomb functions, etc. For example, if the leading terms of the asymptotic solutions are given by the a formula
(z) = -j== exp (±1 (pi'Z — Z^ ln(2pe , (53)
the coefficients of potential in the etalon equation (52) have the form:
Z<i) = 2Zf, Zf = — <d — 3><d — 1> ± 7— — Z2. (54)
' ' 4 pv pi
Step 2. At this step we compute the coefficients ^/(z) and Jji/(z) of the expansion (51) in the form of a series over the inverse powers of z:
kmax Ak') kmax J (k')
**'(*)= £ -jjr, Jj^)= £ Jkr-. (55)
k =0 k =0
After the substitution of Eqs. (51), (55) into Eq. (49) with the use of Eq. (52) and equating the coefficients at z~k R^ (z) and z~k ^dJ^, we arrive at the set of recurrent relations at k' < kmax:
(40) — 2 E + pi) ) + — Z^) $-i) — 2pi(k' — 1)J2'-i) = — ^),
, \ , , / , , (56)
— 2 E + p?,)^ + 2(k' — 1)^'-i) + — 4^) Jii'-i) = — ( )
where the right-hand sides ¡¡k, ) and g[k, ) are defined by the relations
k
k') - _ib'^ib' rhJk-2)
fS') = —( k' — 2)(k' — -2) + £ (v<k) — $-k) +
k=2
k' N k'
+ £(z^ (2k' — 2 — k) Jk-k-i) + £ ( £ 2^Zf -k~k") — k=i j=i,j=i k" = i
— 2p|Qiî)jjk-k) + Q(k)(—2k' + k + d + !)$'-k-i) + vgUjï-k))); (57)
$) = —( k' —1)( k' — 3 + d)jf-2) + £ (vik) — Z[k)) -k)+
k=2
N k'
+ E E(2QV4-fc) -Q^W + d - 3 -fc-1) + vg)*%-fc)) j=1,j=i k=1
with the initial conditions p2 = 2E — e(,0), = Sw, = 0, at z' = z0 span over
the open channels i0 = 1,..., N0 and pv = iqv, qv > 0, q2 = e^,0) — 2E at i' = ic span over the closed channels ic = Na + 1,..., N that followed from (15) and (24). Also from Eq. (56) at k' = 1 and i = z',
(41) -Z^ — o, (41) 0)-
0,
(58)
we obtain the condition zj,1) — e^1.
(k') (k')
Step 3. Here we perform the calculation of the coefficients ) and ) using step-by-step procedure of solving Eqs. (56) for 2 E — e(0), i — z' and k' — 1,..., kmax:
—
,(0) ,(0)
1
-/f0 - 41) -Z< ^-1) + 2(k' - 1)tf-1)
(k') _
,(0) J0)
and for 2 E — e(0), z — z' and k' — 2,... ,km
1
-9¡fc) - 2(k' - 1)^'-1) -(41) -Zi1))
(1A „/.(fc'-1)'
(59)
4*-1) — - [2(k' - 1)]-1 ),
t
(fc'-1)
2(k' - 1) [2E -e
(0)
1
f(k' ) J i'i' .
(60)
The algorithm described above was implemented in Maple and Fortran. The resulting output provided the evaluation of the required solutions Xji' (z) and their
derivatives
Xji' 0) dz
This algorithm has been examined with the results of Ref. [4].
Remark 2. The choice of the appropriate values zmin and zmax for the constructed expansions of the linearly independent solutions for pio > 0 is controlled by the fulfillment of the Wronskian condition (19), (20)
Wr(Q( z); x*(z), X(z)) — ±2zI0 up to the prescribed precision eWr •
(61)
5. Benchmark Calculations of Transmission Coefficient
The wave function ^(x, z) of two particles (or ions) labeled by i = 1,2 coupled by the oscillator potential, penetrating through repulsive (Coulomb) barriers U(xi) in the center-of-mass coordinate frame satisfies the two-dimensional Schrodinger equation [12]:
h2 d2
2Md z2 2.dx2 2
h2 ß2 . _ _ \ _ h U + .Ù2X2 + ^1(X1) + ^2(X2) -Ej *(X, Z)
0,
(62)
where Q is the oscillator frequency, E is the energy, x1 = z + s1x, x2 = z — s3x are variables in the laboratory system of coordinates. The parameters s1 = m2/M, s3 = m1/M are defined via the masses of particles m1 and m2, and their total mass M = m1 + m2 and the reduced mass ^ = m1m2/M.
After the transformation of variables
x = x-}cx, z = M/px-}cz, (63)
with the oscillator unit of length xosc = h/(pz), the corresponding Eq. (62) leads
to the following dimensionless equation in the form of Eq. (1) with fi( z) = f2(z) =
h(z) = f 4 (x) = h(x) = 1:
( 1 d d 1 \
/i( ¿Odz az h(z) J (64)
1 d d (6) L(x;z) = -U(x) dxH(x)~d~x + /3(")y(x,^
where £ = 2E = E/Eosc and V(x, z) are the dimensionless energy and barrier potential in units of energy Eosc = hu/2
V(x, z)=x2 + Ui(xi) + U2(x2) = Eosc-1 (Vz2x2/2 + Ui (xi) + U . (65)
where xi = s2z + six and x2 = s2z — s3x with s2 = \Jp/M.
Model A. We choose the barrier potentials Ui(xi) with the effective charges > 0 in the form of the repulsive truncated Coulomb potential cut off at small 0 < xmin < 1 and large xmax > 1 distances from xi = 0 as [13]
Ui(xi) = —t-(2Z\n --V — . (66)
min(max(xmin, |xi|),xmax) xmax
Model B. We define the Coulomb-like potentials Ui(xi) that depend on the integer parameter s > 2 and the truncation parameter xmin > 0 [12]:
2 Z-
Ui(xi) = . * . (67)
V |xi| + xmin
Model C. For the ionization problem of the hydrogen atom in a magnetic field we choose the total potential U(x, z) of the three-dimensional boundary value problem for Eq. (64) in cylindrical coordinates ( z,x = (p)) with fi(z) = f2(z) = f3(z) = 1 and f4(x) = f5(x) = p at fixed magnetic quantum number m, given by the sum of circular harmonic oscillator potential with frequency z = 7/2 and the three-dimensional Coulomb potential with charge Z = —1, V(x, z) = m2 /p2 + z2p2 — 1 /\/p2 + z2.
The asymptotic boundary conditions for the solution z,x) = {^io(z,x)}^=i with the direction = can be written in the obvious form
(0) exP U (pioz - sign( z) ffl2 ln(2pio |z|) (z ^ -œ,x) ^ BfJ(x)-^--LL+
yPio
(0) exp ( — (pjz - sign(^fi2 ln(2Pjk|)),
+ £Bj0)(x)—--J_ARjio, (68)
3=1 3
No
^ (0) exp(^\p3z - signWln(2Pjk|))) o(z ^ +œ,x) ^Y.B(0)(x)-^--11T3
3=1 VPi
Gusev A. A. Algorithm for Computing Wave Functions, Reflection and . . . 107
0 ( z,x ^ ±ro) ^ 0.
Here Na is the number of open channels at fixed energy 2 E = p2+e(0) > 0; Z12 = 0 for model A and Zi2 = (Z + Z2)/s2 for model B; Rj¿o and Tjio are unknown reflection and transmission amplitudes; ßj0) (x) are the basis functions of the oscillator corresponding to the energy ej0) = 2j — 1 at j > 1
1 9 " — 9 ,2 f0) \ R(0), v _
:h(x)- + x2 - 40) B)0)(X) = 0
)
f4 (x) 9x 9x
Xraax(Z) (g9)
y Bt(0)(x)B¡0)(x) /4(x)dx = Sij.
Note, that the desired solution of the BVP (64), (68) can be constructed in the form of Galerkin-type expansion:
N
^(x, Z) = ^B(0\x)Xji'(z). (70)
j=i
over the basis functions (69). In this case we also obtain the BVP for a set of N coupled ODEs (6), but instead of the effective potentials (7) we have the following ones:
^max (z)
Vij(z) = ef èl3 + J B(0)(x)V(x, z)B(0)(x)f4(x)dx, Ql3(z) = 0. (71)
^min (z)
The effective potentials Vjj and Vjj-2 from (7) and (71) are shown in Fig. 2. One can see that the effective potentials have similar behavior at increasing j and |z| and the diagonal ones Vjj converge to the harmonic oscillator eigenvalues ej0). The nondiagonal effective potentials in the Kantorovich method is smaller in absolute value within the interval of the longitudinal variable, except for the sharp peaks, which are produced by the avoiding crossings of the potential curves £j (z) of the parametric problem (4) with respect to the longitudinal variable.
The eigenvalue problem (4), (5) at a = 0, b(z) = 1 was solved using the ODPEVP program [23] for Z G [¿min, ¿max ] that yielded well-separated eigenvalues | s¿(z) — e¿-i(z)l > e > 0 where e ~ 0.05 for the double-precision arithmetic. This condition is valid for the accepted values of parameters of the considered models. In the case of poorly-separated eigenvalues, i.e. if 0 < |e¿(z*) — £i-1(z*)l < e, one should generate a more dense grid in the vicinity vz = |z — z*| < e* of the avoided crossing points z* and/or use multi-precision arithmetic. For long-range potentials one should construct the appropriate asymptotic expansion for the eigenvalues and the corresponding eigenfunctions z G (—to, +to)\[zmin, ¿max] to build up the asymptotic effective potentials with the leading terms
Vij(z) = (ef + sign(z)2^) ôZj + 0(z-3), Qij(z) = 0(z-3), (72)
and the asymptotic expansions of solutions of the above set of equations (49) at d = 1 (see more details in [12]).
Figure 2. The even effective potentials Vjj and Vjj-2 for Kantorovich (left panel) and Galerkin-type (right panel) expansions for Coulomb-like potential
barriers (67)
For the given number N of (4), the values xmin and xmax of the grid fi^ {xmin, xmax} were chosen in the region |x| > x0 = \/2N + 1, where the Hermite polynomial [24] (or the basis function Bj(x; z) in a general case) has no zeros. These values are computed from the condition
exp | — J dxy x2 — x0 | < e
x0
which in the given case leads to the inequality
J dx^Jx2 - x0 I < enum, (73)
2-,0/2ÏÏx + \/x2 -x2) 0/ xoœ0/2 < enum. (74)
exp ^-x^x2 - x2/2^ ^x + \Jx2 - x0^
To find an approximate solution, at the first step we choose the initial approximation xmax = x0, after which it is increased with step 1 until the condition (74) is satisfied. Values zmin < xmin and zmax > xmax were chosen from the condition that the potential (66) or (67) is negligible within the interval
x^Qin << x << xmax.
The matching points ^match and 2match of the numerical (7) and asymptotic (50), (72) effective potential were calculated as follows:
zmatch = min (zf, z^ , ^match = max (zf, z2^ ,
\r)(fcmaxh / |T/(fcmax) I (75)
zf = sign(z) fcmai/ i| , z2 = sign(z) fcmai/,
X
since < |Q
I T/"(^max)
«33' i - i^NN-1b IVjj>
satisfied by the inequalities zmin < z1 should be calculated from the conditions
I < |VÄ-x)|. So, the values z,
< ®min and ¿max > ¿2
and zmax were
match " - ■ and Zmax > ¿match > Imax that
¿min = min
max
match
n ,
min
^match, max j
(
j^ o
x)|
( ^max) jio
, mm | —
max
(
j^ o
(
j^o
(76)
In the considered examples we used the grids Qhx {xmin,xmax} = {-10(768)10} and {¿min, ¿max} = {-125(200)-25(100)-6(200)6(100)25(200)125} with the Lagrange elements of the order p = 4 between the nodes. In the above notations the number of grid elements for grids Qhx and Qhz is shown in the parentheses.
For the calculation of asymptotic solutions of model B, we used the etalon equation (52) at d = 1, k'max = 1 and Z(1) = 2sign(y)Z12, which corresponds to the known solutions in the open channels
R± (Pio, z) = P
-1/2 { (Go(Pio, +¿) ± ïFo(Pio, +¿))exp(^г<5ío)/2, I (Go(Pio, — z) T iF0(pio, — z))exp(±^o)/2,
z > 0,
z < 0,
and in the closed channels
- „-1/2
Rlc (qlc, z) = q-1/2t exp(—t/2) U(1 + Zi2/qic, 2, t), t = 2qlc |z|
(77)
(78)
Here F0(pio, z) and G0(pio, z) are the regular and irregular continuum zero-order Coulomb functions calculated by the subroutine RCWFNN [4] which is a modified version of the subroutine RCWFN [25] for the DOUBLE PRECISION accuracy; aio = argr (1 + iZ12/pio ) is the Coulomb phase shift [24]; and U(a, b, c) is the confluent hypergeometric function of the second kind, calculated by the subroutine CHGU [26]. Note, that for the numerical calculation we have neglected the exponentially small factor exp(-1/2) in Ric (qic, z) and its first derivative, since this factor is canceled during evaluation of the z) matrix in Eq. (25). At the boundary points zmin and zmax the absolute accuracy sWr of the calculated Wronskian was less then 10-11.
Below we use the parameters: m1 = m2 = 1, Xmin = 0.1, Z1 = Z2 = 0.5 and Z1 = Z2 = 1. Also Xmax = 5 for model A, and s = 8 for model B.
The total probabilities T = rT11 = N^ 1^1'|2 of penetration through truncated Coulomb (66) and Coulomb-like (67) potential barriers of models A and B are shown in Fig. 3. These figures illustrate the important peculiarity that a more realistic nontruncated Coulomb-like barrier, being wider than the truncated one, leads to a set of the probability maximal having a bigger half-width. This fact can be used for distinguishing between the models of type A and B by observing the quantum transparency effect.
The results given by Kantorovich and Galerkin-type expansions differ by less the 0.1%. The absolute maximum value max |%iio | = max |%iio( z)l of components xu0(¿) displayed in Fig. 4 show similar convergence of expansions (3) and (70). For the chosen approximation of Lagrange elements of the order = 4 for their eigenfunctions and transmission coefficient we obtain numerical estimations of Runge coefficients within 4.6 ^ 4.7, and in the range 7.4 ^ 7.7, respectively, which correspond to the theoretical error estimates (35) at fixed number N of Eqs. (6).
Note, that in this example, the Galerkin and Kantorovich expansions have approximately similar rate of convergence, which gives preference to the first method, since
k
k
max
1.0 0.80.60.40.20.0
Zl=Z2=05
m1=1 m2=1 x =0.1
min
x =5
max
1
3 5 7 9 11 13 15 17 19 21 23 2E
1.0 0.80.60.40.20.0
Z1=Z2=1
m1=1
m2=1
x =0.1 min
x =5 max 1 J. ll i A
3 5 7 9 11 13 15 17 19 21 23 2E
1.U-,
Z1=Z2= =0.5
0.8- m1=1
m2=1
0.6- x =0 min 1
0.4-
0.2-
00- J w
1.U-, z1=z2=1
0.8- m1=1
0.6- m2=1 x =0.1 min
0.4-
0.2-
0 0-
5 7 9 11 13 15 17 19 21 23 2E
5 7 9 11 13 15 17 19 21 23 2E
Figure 3. The total probabilities T = \T\\1 of penetration vs energy E = 2E through truncated Coulomb (66) (upper panel) and Coulomb-like (67) (lower
panel) potential barriers
Figure 4. The absolute maximum value max lxu01 vs of number i component of continuum spectrum solution by Kantorovich (K) and Galerkin-type (G)
expansions
the procedure for calculating the matrix elements by direct integration with known basis functions is less expensive than the procedure for calculating the matrix elements with the help of the basis functions and their first derivatives, obtained by solving the parametric boundary-value problems arising in the Kantorovich method. This is due to the fact that in the above examples, these basic functions of a finite set used in the Galerkin-type and Kantorovich expansions are transversely localized in the region, bounded by the scope of the oscillator potential.
For problems similar to the ionization of a hydrogen atom in the uniform magnetic field the total potential have spherical symmetry in a vicinity of the Coulomb center where the oscillator potential is negligible and cylindrical symmetry at a large distance from the center with increasing of the oscillator potential. In this case Galerkin-type expansion has a rather slow convergence rate and we apply the decomposition of the desired solution in cylindrical or spherical coordinates over the surface orthogonal basic functions depending on the longitudinal variable z or radial variable r as a parameter, respectively.
As an example, Fig. 5 shows the results of calculations of the transmission coefficient of the first open channel obtained by the algorithm discussed above and implements its program KANTBP 3.0 [16] in cylindrical and spherical coordinates of model C.
Figure 5. The transmission coefficient |T|2 = |T|n vs (E2 - 2E)-1/2 for 10, 20, 40
cylindrical adiabatic basis functions (dashed lines) and for 10 spherical adiabatic basis functions (solid line) at to =0; the arrow points at the energy of the first Landau threshold (E2 - £1 )-1/2 = « 2.236 (left panel) and the bounce ionization cross-section o~d(w) vs energy E from lower states of a hydrogen atom |j = 1,k = 1, to, a = +1} to continuum states IE, to) to = 0, -1,..., -12 (right panel).
Here 7 = 0.1 (B = 2.35 • 104T), Z = 1
One can see that for small values of the magnetic quantum number |m| the wave function localized in a region of the Coulomb center and the convergence rate of the Kantorovich expansion in spherical coordinates is higher than in cylindrical ones. However, with increasing of the magnetic quantum number |m| the wave function is displaced from the center of the Coulomb and the situation is reversed [27]. Nonmonotonic behavior of the transmission coefficient and the ionization cross-section vs energy shown Fig. 5 are produced by the series of metastable states embedded in the continuum spectrum due to contribution of close channels up to |rn| < 12. Decreasing of ionization cross-section at |rn| > 12 follows from adiabatic separation of transversal and longitudinal variable [27] at values of adiabatic parameter (wp/wz,i=i)4/3 = |to|71/3 > 6, where = 7/2, in agreement with [6]. The analysis of the models shows that solution of the multichannel scattering problem for different compound potentials is achieved with help of the proposed algorithms and software complexes [8,16,17].
6. Conclusion
A new formulation of the problem related to the multichannel scattering problem for a system of the second-order ordinary differential equations containing the potential matrix elements and first-derivative coupling terms with respect to the independent variable belonging to axis is proposed. The boundary problems for the coupled second-order differential equations are solved by the finite element method using high-order accuracy approximations with Lagrange elements. The order of approximation depends on the smoothness of required solution. Taking into account matrix of logarithmic derivatives zmjn) (or zmax)) determined by asymptotic regular solutions (25), algebraic problems (30)-(31) with respect to pair unknowns zmax), (or zmin), together with the reflection R^ (R^) and transmission T^ (T^)
matrices, arising after the corresponding replacement of the scattering boundary problem in open channels at fixed energy value, E, was reduced to the auxiliary system of nonhomogeneous algebraic equations (39) (or (47)). The later has been solved by the L D LT factorization of symmetric matrix and back-substitution methods using the modifications of the DECOMP and REDBAK programs [28], respectively without calculation of the inverse of submatrix G^ (or G^). Presented algorithm is implemented as Fortran program KANTBP 3.0 [16].
This approach can be used in the calculations of quantum transmittance of barriers for composite particles [14,15], channeling problem [5], quantum well and wire problems with Hydrogen-like impurities [29,30], and resonant molecule formation in waveguides [31].
References
1. Born M., Huang K. Dynamical Theory of Crystal Lattices. — New York, Oxford Clarendon Press, 1964.
2. Goodvin G. L., Shegelski M. R. A. Three-Dimensional Tunneling of a Diatomic Molecule Incident Upon a Potential Barrier // Phys. Rev. A. — 2005. — Vol. 72. — Pp. 042713-1-7.
3. Ринг Р., Расмуссен Д., Массман Г. Проблемы проницаемости неодномерных барьеров // ЭЧАЯ. — 1976. — Т. 7. — С. 916-951. [Barrier Penetration Theory in more than One Dimension / R. Ring, J. Rasmussen, H. Massman // Sov. J. Part. Nucl. — 1976. — Vol. 7. — P. 916-951. — (in russian). ]
4. POTHMF: a Program for Computing Potential Curves and Matrix Elements of the Coupled Adiabatic Radial Equations for a Hydrogen-Like Atom in a Homogeneous Magnetic Field / O. Chuluunbaatar, A. A. Gusev, V. P. Gerdt et al. // Comput. Phys. Commun. — 2008. — Vol. 178. — Pp. 301-330.
5. Channeling Problem for Charged Particles Produced by Confining Environment / O. Chuluunbaatar, A. A. Gusev, V. L. Derbov et al. // Phys. Atom. Nucl. — 2009. — Vol. 72. — Pp. 811-821.
6. Guest J. R., Raithel G. High-|m| Rydberg States in Strong Magnetic Fields // Phys. Rev. A. — 2003. — Vol. 68. — Pp. 052502-1-9.
7. Kantorovich L. V., Krylov V. I. Approximate Methods of Higher Analysis. — New York: Wiley, 1964.
8. Gusev A. A. The Algorithms of the Numerical Solution to the Parametric Two-Dimensional Boundary-Value Problem and Calculation Derivative of Solution with Respect to the Parameter and Matrix Elements by the Finite-Element Method // Bulletin of PFUR. Series "Mathematics. Information Sciences. Physics". — 2013. — No 4. — Pp. 101-121.
9. Strang G., Fix G. J. An Analysis of the Finite Element Method. — New York: Prentice-Hall, Englewood Cliffs, 1973.
10. Ramdas Ram-Mohan L. Finite Element and Boundary Element Applications in Quantum Mechanics. — New York: Oxford Univ. Press, 2002.
11. Finite-Element Solution of the Coupled-Channel Schrodinger Equation using High-Order Accuracy Approximations / A. G. Abrashkevich, D. G. Abrashke-vich, M. S. Kaschiev, I. V. Puzynin // Comput. Phys. Commun. — 1995. — Vol. 85. — Pp. 40-64.
12. Symbolic-Numerical Algorithms to Solve the Quantum Tunneling Problem for a Coupled Pair of Ions / A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar et al. // Lecture Notes in Computer Science. — 2011. — Vol. 6885. — Pp. 175-191.
13. Пеньков Ф. М. Квантовая прозрачность барьеров для структурных частиц // ЖЭТФ. — 2000. — Т. 118. — С. 806-815. [Pen'kov F. M. Quantum Transmittance of Barriers for Composite Particles // ZHETF. — 2000. — Vol. 91. — Pp. 698705. — (in russian). ]
14. Гусев А. А. Новый метод построения осцилляторных функций квантовой системы тождественных частиц в симметризованных координатах // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2013. — № 3. — С. 5267. [Gusev A. A. New Method for Constuctiond the Oscillator Functions of a Quantum System of Identical Particles in Symmetrized Coordinates // Bulletin of PFUR. Series "Mathematics. Information Sciences. Physics". — 2013. — No 3. — Pp. 52-67. — (in russian). ]
15. Гусев А. А. Модель туннелирования кластеров через отталкивающие барьеры в представлении симметризованных координат // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2014. — № 1. — С. 5473. [Gusev A. A. The Model of Tunneling of Clusters Through Repulsive Barriers in Symmetrized Coordinates Representation // Bulletin of PFUR. Series "Mathematics. Information Sciences. Physics". — 2014. — No 1. — Pp. 54-72. — (in russian). ]
16. Chuluunbaatar O, Gusev A. A., Vinitsky S. I., Abrashkevich A. G. A Program Package for Solution of Two-Dimensional Discrete and Continuum Spectra Boundary-Value Problems in Kantorovich (Adiabatic) Approach. — 2013. — http://wwwinfo.jinr.ru/programs/jinrlib/kantbp/indexe.html.
17. KANTBP: A Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adiabatic Approach / O. Chuluunbaatar, A. Gusev, A. Abrashkevich et al. // Comput. Phys. Commun. — 2007. — Vol. 177. — Pp. 649-675.
18. KANTBP 2.0: New Version of a Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adia-batic Approach / O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich // Comput. Phys. Commun. — 2008. — Vol. 179. — Pp. 685-693.
19. Macek J. Hermitian R Matrix in the Presence of First-Derivative Couplings // Phys. Rev. A. — 1984. — Vol. 30. — Pp. 1277-1278.
20. Calculation of a Hydrogen Atom Photoionization in a Strong Magnetic Field by using the Angular Oblate Spheroidal Functions / O. Chuluunbaatar, A. A. Gusev, V. L. Derbov et al. // J. Phys. A. — 2007. — Vol. 40. — Pp. 11485-11524.
21. A Symbolic-Numerical Algorithm for Solving the Eigenvalue Problem for a Hydrogen Atom in the Magnetic Field: Cylindrical Coordinates / O. Chuluunbaatar, A. Gusev, V. Gerdt et al. // Lecture Notes in Computer Science. — 2007. — Vol. 4770. — Pp. 118-133.
22. Golub G. H., Van Loan C. F. Matrix Computations. — New York, Johns Hopkins Univ. Press, 1996.
23. 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. Abrashkevich // Comput. Phys. Commun. — 2009. — Vol. 180. — Pp. 13581375.
24. Абрамовиц М., СтиганИ. Справочник по специальным функциям. — М.: Наука, 1979. [Abramovits M., Stigun I. A. Handbook of Mathematical Functions. — New York: Dover, 1972. ]
25. Coulomb Wave Functions for All Real r] and p / A. R. Barnett, D. H. Feng, J. W. Steed, L. J. B. Goldfarb // Comput. Phys. Commun. — 1974. — Vol. 8. — Pp. 377-395.
26. FORTRAN Routines for Computation of Special Functions. — http://jin.ece. illinois.edu/routines/routines.html.
27. Symbolic-Numerical Calculations of High-|m| Rydberg States and Decay Rates in Strong Magnetic Fields / A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar et al. // Lecture Notes in Computer Science. — 2012. — Vol. 7442. — Pp. 155-171.
28. Bathe K. J. Finite Element Procedures in Engineering Analysis. — New York: Englewood Cliffs, Prentice Hall, 1982.
29. Symbolic-Numerical Algorithms for Solving Parabolic Quantum Well Problem with Hydrogen-Like Impurity / S. I. Vinitsky, O. Chuluunbaatar, V. P. Gerdt et al. // Lecture Notes in Computer Science. — 2009. — Vol. 5473. — Pp. 334349.
30. The Application of Adiabatic Method for the Description of Impurity States in Quantum Nanostructures / A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky et al. // J. Phys. Conf. Ser. — 2010. — Vol. 248. — Pp. 012047-1-8.
31. Melezhik V. S., Schmelcher P. Quantum Dynamics of Resonant Molecule Formation in Waveguides // New J. Phys. — 2009. — Vol. 11. — Pp. 073031-1-10.
УДК 517.958:530.145.6
Алгоритм вычисления волновых функций, матриц отражения и прохождения многоканальной задачи рассеяния в адиабатическом представлении методом конечных элементов
А. А. Гусев
Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, г. Дубна, Московской обл., Россия, 141980
В адиабатическом представлении многоканальная задача рассеяния для многомерного уравнения Шрёдингера сведена к краевой задаче для системы самосопряжённых обыкновенных дифференциальных уравнений второго порядка на конечном интервале с однородными граничными условиями третьего типа в левой и правой граничных точках в рамках метода Канторовича, используя адиабатический базис поверхностных функций, зависящих от продольной переменной как от параметра. Для искомых решений краевой задачи сформулированы однородные условия третьего рода, используя известные наборы линейно-независимых регулярных и нерегулярных асимптотических решений в открытых каналах редуцированной многоканальной задачи рассеяния на оси, в которые входят искомые матрицы амплитуд прохождения и отражения, и набор линейно независимых регулярных асимптотических решений в закрытых каналах. Предложен экономичный и устойчивый алгоритм численного расчёта с заданной точностью матриц отражения и прохождения и соответствующих волновых функций многоканальной задачи рассеяния для системы самосопряжённых обыкновенных дифференциальных уравнений второго порядка с матрицами потенциалов и матрицами, содержащими первые производные, используя аппроксимацию высокого порядка точности методом конечных элементов (МКЭ). Эффективность предложенного алгоритма продемонстрирована решением двумерной квантовой задачи прохождения пары частиц с осцилля-торным потенциалом взаимодействия через отталкивающие потенциалы кулоновского типа и задачи рассеяния электрона в кулоновском поле протона и в однородном магнитном поле в рамках методов Канторовича и галёркинского типа, а также анализом их сходимости.
Ключевые слова: многоканальная задача рассеяния, матрицы амплитуд прохождения и отражения, многомерное уравнение Шрёдингера, адиабатическое представление, метод Канторовича, краевая задача, система самосопряжённых обыкновенных дифференциальных уравнений второго порядка, метод конечных элементов.
