Алгоритм решения двумерной краевой задачи для модели квантового туннелирования двухатомной молекулы через отталкивающие барьеры

Математическое моделирование

УДК 517.958:530.145.6

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

отталкивающие барьеры

А. А. Гусев*, Л. Л. Хай

* Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, г. Дубна, Московской обл., Россия, 141980

^ Белгородский государственный национальный исследовательский университет ул. Победы, д. 85, г. Белгород, Россия, 308015

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

Ключевые слова: квантовое туннелирование, двухатомные молекулы, отталкивающие барьеры, краевые задачи, метод Галёркина, метод Канторовича, асимптотические решения, метод конечных элементов.

1. Введение

Исследования туннелирования связанных частиц через отталкивающие барьеры [1] выявили эффект резонансной квантовой прозрачности барьеров: когда размер кластера сравним с пространственной шириной барьеров, имеют место механизмы, приводящие к большей прозрачности барьеров, подобные механизмам просветлённой оптики. Эти механизмы связаны с формированием барьерных ре-зонансов, обусловленных тем фактом, что потенциальная энергия составной системы имеет локальные минимумы, приводящие к возникновению метастабиль-ных состояний движущегося кластера [2]. В настоящее время этот эффект и его возможные приложения являются предметом интенсивных исследований различных квантовомеханических задач, например, квантовой диффузии молекул [3], резонансного прохождения экситонов через гетероструктурные барьеры [4], резонансное образование молекул из отдельных атомов [5], управление направлением диффузии в твёрдом теле [6], туннелирование ионов и кластеров через отталкивающие барьеры [7-9].

Для анализа подобных задач необходимо разработать модели и численно-аналитические методы, основанные на приближениях, обеспечивающих реалистическое описание взаимодействий как атомов в молекуле, так и с барьерами, и разработать эффективные алгоритмы и комплексы программ.

Статья поступила в редакцию 18 сентября 2014 г.

Авторы благодарят Л. А. Севастьянова, В. Л. Дербова, С. И. Виницкого, П. М. Красовиц-кого, Ф. М. Пенькова, О. Чулуунбаатара за сотрудничество и поддержку. Работа поддержана грантами РФФИ 14-01-00420 и 13-01-00668.

В настоящей работе представлена формулировка и исследование модели квантового туннелирования двухатомной молекулы с ядрами, связанными потенциалом Морзе, через гауссовские барьеры в s-волновом приближении, используя разложения Галёркина и Канторовича искомого решения. Даны формулировки двумерных краевых задач в декартовой и полярной системах координат. Используя различные базисные функции — решения вспомогательных краевых задач по поперечной переменной или по угловой переменной с параметрической зависимостью от радиальной переменной, краевая задача сводится к системе связанных обыкновенных дифференциальных уравнений второго порядка. В декартовых координатах матричные элементы убывают экспоненциально (ниже порога диссоциации), а в полярных координатах — убывают как обратные степени по независимой переменной. Поэтому в последнем случае требуется вычисление асимптотических разложений матричных элементов и фундаментальных решений системы связанных обыкновенных дифференциальных уравнений второго порядка. Для их вычисления разработаны и представлены символьно-численные алгоритмы, реализованные в системе компьютерной алгебры Maple.

Дан сравнительный анализ потенциальных матричных элементов (матриц переменных коэффициентов — эффективных потенциалов систем дифференциальных уравнений), вычисленных в декартовых и полярных координатах, которые использовались для решения задачи квантового туннелирования ниже порога диссоциации. Формулировка краевых задач в двух системах координат необходима для дальнейшего самосогласованного изучения задачи выше порога диссоциации, для которой корректные и удобные для эффективного численного решения задачи краевые условия формулируются в полярной системе координат. Представлен анализ эффекта квантовой прозрачности — резонансной зависимости коэффициента прохождения от энергии налетающей на барьер молекулы.

Структура работы следующая. В разделах 2 и 3 дана формулировка двумерных краевых задач в декартовых и полярных координатах и их редукция методами Галёркина и Канторовича. В разделе 4 представлены ведущие члены асимптотических разложений матриц переменных коэффициентов и фундаментальных решений, дано описание алгоритмов их вычисления до требуемого порядка точности. В разделе 5 анализируются решения краевых задач для модели квантового туннелирования молекулы и эффект квантовой прозрачности барьеров. В заключении указаны возможные применения развитых алгоритмов и комплексов программ.

2. Модель А. Редукция краевой задачи методом

Галёркина

Рассмотрим двумерную модель двух тождественных частиц с массой т, связанных парным потенциалом У(х2 — х{) и взаимодействующих с внешними барьерными потенциалами Уъ(х\) и Уъ(х2). Выполняя замену переменных х = х2 — хл, у = х2 + хл, у е та), х е (—та, та), получаем уравнение Шрёдингера

для волновой функции Ф(ж,у) в ¿¡-волновом приближении

( Н2 1 д . , .д Н2 1 д . . , д , Л т , , ^

шШ ЪуШ 8Гу — шШ ^

где У(х,у) = Ум(ж) + Уь(х{) + Vь(х2) — потенциальная функция, Е энергия системы, Н — постоянная Планка.

Уравнение, описывающее молекулярную подсистему, имеет вид (—+ 1> м м — У)«*>=а

— полная

Предполагается, что молекулярная подсистема имеет дискретный спектр, состоящий из конечного числа п связанных состояний с собственными функциями ф, (х), ] = 1,п, и собственными значениями ё, = — |ё, |, и непрерывный спектр собственных значений £ ^ 0 с соответствующими собственными функциями фц(х).

Асимптотические краевые условия, налагаемые на решения двумерной крае-

вой задачи з-волновым приближением Ф( у,х) = {Ф,(ув асимптотической

области ^^ = {(х, у)||х|/|у| ^ 1}, описывающей движение молекулы (ниже порога диссоциации Е < 0) в направлении V имеют вид [9]:

Ф,(У ^ —<х>,х) ^ ф0{х)--ХРЩ + ^ФМ^^Р-Ш)Я

, ( У ^ ^^^, х) г *ф, (х) -- | / *ф^ (х) -- —^, ,

л/Рз /2

Ф,

где /1(у) = /2(у) = 1, (Е) и Тц(Е) — амплитуды отражения и прохождения,

Nо ^ п — число открытых каналов, р^ = ^(т/Н2)(Е — ё^ > 0 — волновое число, Фз (х) и е, < 0 при ] = 1, п — собственные функции и собственные значения краевой задачи (2).

Решение уравнения (1) ищем в виде разложения Галёркина

^шах

Фга

Здесь (у) — неизвестные функции и ф, (х) — ортонормированные базисные функции, определяемые как собственные функции краевой задачи для уравнения

+

ТЗщёх и(х)1+ ум (х) — £'1ф'

на интервале 0 ^ х ^ хтах с краевыми условиями и условиями нормировки и ортогональности

^шах

ф,

0

где /з(х) = f4(х) = 1, V(х) = (т/Н2)У(х), е, = (т/Н^ё,.

С помощью программы ОБРЕУР численного решения краевых задач методом конечных элементов [10] вычисляем набор п связанных состояний, состоящий из собственных функций ф, (х) и собственных значений е,, '} = 1,п, и требуемый набор псевдосостояний, состоящий из собственных функций ф, (х) и собственных значений е , ^ 0, ] = п + 1,]тах. Набор псевдосостояний аппроксимирует набор собственных решений непрерывного спектра е ^ 0 краевой задачи для уравнения (2). Система связанных уравнений в форме Галёркина имеет вид

1 к( у)^ + £г — Е

Ь(у)ё у ё у

Хгг 0

=1

где эффективные потенциалы Vij(y) даются интегралами

Vb

о

В результате задача рассеяния (1)—(2) с асимптотическими условиями (3) сведена к краевой задаче для системы связанных уравнений в форме Галёркина (7) для /i( у) = f 2(у) = 1 с краевыми условиями при у = j/min и у = ?/max [8,11]:

dF

dy

yt)F

Здесь утгп) И у тах) симметричные матрицы размерностью ,^шах х .^шах;

зависящие от энергии е, ^(у) = (у)}^о°=1 = о(у)}7=а1х— искомая

матрица решений размерностью ]шах X где н0 = шах; ^ ^шах — число открытых каналов. Эти матрицы и матрицы амплитуд отражения и прохождения К и Т размерностью м0 X м0 вычисляются с заданной точностью с помощью программы КЛМТБР 3.0 [12], реализующей метод конечных элементов [13].

В качестве потенциала взаимодействия Vм (ж) был выбран потенциал Морзе

Ум

Решения дискретного спектра краевой задачи (5)—(6) с потенциалом (10), с требуемой точностью (< 10-10 при выбранных значениях параметров) аппроксимируется известным спектром краевой задачи для уравнения (2):

= -D

1 - S

2 Г 1

1 -1 i 1 , з = 1,...,п = ? +2

Собственные функции дискретного спектра фj (ж) (5)—(6) аппроксимируются решениями фj (Q) уравнения (2) по новой переменной Q:

2A.frt лл1.!?\ ( Л , „. 1 /О s2\

^(0 + 1#Ж) + / 1 , 3 + *з - 1/2 (Л) = 0

+ + Г 4+ с (с) = 0,

где = V-7/Р = -3 + 1/2 и

соответствуют расширенному интервалу до всей оси ж £ (-та, +то) и имеют вид

«C) = NexP

Набор псевдосостояний с собственными функциями ф7 (ж) и собственными значениями £7 ^ 0, ^ = п + 1,^шах, аппроксимируется набором решений непрерывного спектра фк (С) при фиксированных к = > 0 уравнения (2) по новой переменной

¿2фк(0 + ^(С) + / 1 + VЪ/р.

"Ж" + + Г4 + ~ + ё) фк(С) = 0.

При фиксированном в к = к эти решения имеют вид

Р

- exp-wK^F, + 1 + f, 1 + ™ А) ,

и-f))-

w = arg ( Г ( 1 + ^ ) ) + arg ( Г ( + 1 -

Асимптотически ф^(х ^ ж) = вш(кх + д(к)), ё(к) = —кхеч — 8к Ы(2л/И)/р) + т соответствует фазе рассеяния.

Для тестовых вычислений использовались следующие параметры для молекулы Ве2: приведённая масса р = т/2 = 4,506Ба, среднее расстояние между ядрами 2,47А, частота вибрационных колебаний молекулы в температурных единицах Нш = 398, 72К, основное состояние молекулы 1Х+, волновое число 277,124сш-1 для наблюдаемых переходов из возбуждённого в основное состояние (использовано соотношение 1К = 0, 69503476 сш-1 из [14]). Зная средние размеры молекулы и учитывая расстояния между уровнями энергии, параметризуем молекулярный потенциал, фитируя наблюдаемые величины: Б = 1280К, хея = 2,47А, р = 2, 968А.-1 определяется из условия (ё2 — ё1)/(2пНс) = 277,124 сш-1, ? =

рН/^тБ = 0,193 — безразмерный параметр задачи, и Б = ^тБ/Н^ =

(р/0,193)2 = (2, 968А-1/0,193)2 = 236, 5А-2. Согласно (11) энергия основного состояния молекулы Ве2 равна — £1 = —1044, 88К.

Поскольку химическая связь в молекуле Ве2 — ван-дер-вальсовского типа, то можно рассматривать каждый составляющий атом как независимо взаимодействующий с внешним барьерным потенциалом. Барьерный потенциал был выбран так, что его высота и ширина соответствует типичным барьерам в реальной кристаллической решётки. Более того, этот потенциал должен быть гладкой функцией, имеющей вторую производную, чтобы было можно применять высокоточные численные методы, такие как метод Нумерова или метод конечных элементов, при решении краевых задач для систем обыкновенных дифференциальных уравнений второго порядка. Барьерный отталкивающий потенциал был выбран в гауссовской форме:

Vb

Здесь Vo = 1280K, D = 236, 510003758401A"2 = (m/h2)Vo, a = 5, 23  10"2A2 определялись из требования, чтобы ширина отталкивающего потенциала при кинетической энергии, равной энергии основного состояния молекулы, равнялась 1A, так что среднее расстояние 2,47A между атомами молекулы Be было бы меньше, чем расстояние 2, 56A между атомами Cu atoms в плоскости (111) ячейки

кристаллической решётки. Высота потенциального барьера Vo порядка 200 meV оценивалась, следуя экспериментально наблюдаемой квантовой диффузии атомов водорода на поверхности меди [15].

2

На рис. 1 показаны графики гауссовского барьера и потенциала Морзе в А

На рис. 2 представлены сечения полной потенциальной энергии, вычисленные собственные функции краевой задачи (5) и эффективные потенциалы Vij(y) из уравнения (8).

Рис. 1. Гауссов барьер Vь(жг) = .Оехр

Ь = 236, 510003758401А-2 = (т/П2)Уо = (т/Ь2)В, Уо = В = 1280К, а = 5, 23 • 10-2А2 и двухчастичный потенциал взаимодействия

Vм(ж) = Ь{ехр[-2(ж - хеч)р] - 2ехр[-(ж - хеч)/3]}, хеч = 2, 47А, р = 2, 96812423381643А.-1

Рис. 2. Сечения полной потенциальной энергии V(у;х) = Vм(ж) + Уь(у;х) при

у = 2, 2; 2, 3; 2, 4; 2, 6; 2, 8; 3; 3, 5; 4

Заметим, что волновые функции ф7(ж) и собственные значения £7 связанных состояний '} = 1, 5 (непрерывные линии) аппроксимируют известные аналитические результаты краевой задачи для уравнения (2) с потенциалом Морзе (10) с четырьмя и семью значащими цифрами соответственно. Эти состояния локализованы в потенциальной яме, тогда как волновые функции псевдосостояний '} = 6,..., 12 аппроксимируются примерно с той же точностью, но локализованы вне потенциальной ямы. Из рисунка видно, что матричные элементы между связанными состояниями локализованы в окрестности барьеров, а матричные элементы между псевдосостояниями локализованы вне барьеров. Матричные элементы между связанными и псевдосостояниями малы. Решения краевой задачи

(5), (6) вычислялись с помощью фортрановской программы ODPEVP [10] на конечно-элементной сетке Q,x = [0(Neiem = 800)12} с Neiem = 800 лагранжевыми элементами четвёртого порядка р = 4.

3. Модель Б. Редукция краевой задачи методом

Канторовича

Используя замену переменных х = р 8Шф, у = р совф, перепишем уравнение

1 д д -Р

р др др р2 дф

1 д 2

-2 Я^ + V (РФ -Е) Ф(р,ф)=0,

где потенциальная функция V

Сечения потенциальной функции V(р, ф) при наборе значений медленной переменной р показаны на рис. 3. Заметим, что при больших р ширина потенциальных ям убывает с ростом р. Следовательно, при больших р потенциал задачи двух центров, симметричный по отношению ф = л/2, трансформируется в пару одноцентровых потенциалов Морзе.

Рис. 3. Сечения полной потенциальной энергии V(р; = Vм(р; + УЪ(р\ у) в полярных координатах при р = 2, 2; 2, 3; 2, 4; 2, 6; 2, 8; 3; 5; 10 (соответствующие кривые обозначены номерами 1,...,8). Горизонтальные линии — энергии уровней е(р)/р2 (в А-2) при р =10

Асимптотические краевые условия наложенные на искомое решение двумерной модели в ¿¡-волновом приближении Ф(р, ф) = (Ф¡(р, в асимптотической области = (( ф,р)\ф/р ^ 1} записываются в виде

Na

ф(р,ф,ф0) = V Фуio (Р,Ф)Фга (-'-Рв;Р ^ +Ж),

/ у г о

¿о = 1

2 *

Ф о

=1

Ф*о

exP

где угловая переменная ро определяет направление падающей волны, в частности, ро = 0 соответствует v и ро = ж соответствует v Sj¿о (Е) — элементы S-матрицы размерностью N0 х Na, Na — число открытых каналов,

Pi = \J (т/Н2)(Ё — ) > 0 — волновое число. Ниже порога диссоциации Е < 0,

= £ < 0 собственные функции,

локализованы в асимптотических областях s.

Решение уравнения (15) ищем в виде разложения Канторовича

jmax

Ф о

j=l

Здесь Xjiо (р) — неизвестные функции и фj (р; р) — ортонормированные базисные функции на интервале р £ [0, ж] определяются как собственные функции краевой задачи для уравнения

подчинённые краевым условиям первого рода и условиям ортогональности и нормировки

к

ф

Собственные функции фj(р;р), ] = 1,20 показаны на рис. 4 при р = 3 и р = 10. Учитывая упомянутую симметрию потенциала V(р, р) = V(т — р, р), набор собственных функций расщепляется на два поднабора, именно, чётные ф^=1(р; р) и нечётные ф^=-1(р;р) функции.

Чётные и нечётные решения вычислялись с помощью фортрановской программы ОБРЕУР [10] на конечноэлементных сетках, выбираемых в зависимости от значения параметра р3 = (8 + рхед)/(рр) > т/4. Если р3 > т/4, то использовалась равномерная сетка = {р1(Ме1 ет = 3200)т/2}, иначе использовалась квазиравномерная сетка 0 ^ = {р1( Ме1 ет = 1200)р2( Nеlет = 240)р4^е1ет = 160)рб(^еето = 400)т/2}. Здесь значения угловой переменной р1 = (—3+рхед)/(рр) и р2 = (4 + рхея) /(рр) определяются левой и правой границами потенциальной ямы (16), а значения р4 = т/4 — 4^[а/р и р5 = т/4 + 4^[а/р — левой и правой границами потенциального барьера (16).

Система связанных самосопряжённых уравнений в форме Канторовича имеет вид

i d d + ^(р) Е jmax

— 1~рт~ +--2--Е

d р d р р2

Хгго

=l

где матричный оператор Wij^) даётся выражением

^•

Рис. 4. Чётные и нечётные собственные функции параметрической задачи на собственные значения для быстрой подсистемы при р =3 и р =10

Рис. 5. Чётные эффективные потенциалы Нг^

Потенциальные кривые (термы) е^(р) (см. рис. 6) и эффективные потенциалы Яц(р) = (р), Нц(р) = (р) и VЪ(р) (см. рис. 7-8) даются интегралами

к

Я,ЛР) = -) МШ**^, Н,М = !йф¿^¿ША,

0 0

к

VIр ^¡Лфф^^и^р™(ф+л/4У\ + ^ (рМф-1/4>)) Ф,фр).

0

Рис. 6. Потенциальные кривые е^

Рис. 7. Чётные эффективные потенциалы Q¡j

Рис. 8. Чётные эффективные потенциалы Уг^

Линейные комбинации (р; р) = (ф^=1(р;р) ± ф^=-1(р;р))/\/2 при больших р имеют максимумы в окрестности р = 0 и р = т соответственно, так что

они соответствуют функциям, представленным на рис. 2. Учитывая это свойство, имеем выражения [16]

T =

которые связывают чётные S+i и нечётные S-i элементы матрицы S = em/4Sегж/4 из уравнения (18) с элементами матриц амплитудами прохождения T и отражения R из асимптотик искомого решения (3).

Таким образом, задача рассеяния для уравнения (15) с асимптотическими краевыми условиями (18) сведена к краевой задаче для системы связанных уравнений в форме Канторовича (17) с краевыми условиями при p = pmin и р = pmax [11]:

^^ =

dp p=P,

где '(p) — неизвестная jmax x jmax симметричная матричная функция, F(p) = [Xia(p)}j0=i = {{Xji0(p)}j=aix}ij=i — искомое матричное решение размерностью 3 max x N0 вычислялись, используя фортрановскую программу KANTBP 3.0 [12].

4. Асимптотики эффективных потенциалов и решений

Алгоритм 1. При больших р ширина потенциальной ямы убывает с ростом р (см. рис. 3). Этот факт позволяет линеаризовать аргумент р эшф — хея ^ р(ср — агс8ш(Хед/р)) при |х — хед|/р ^ 1 в выражении потенциальной функции Vм (р эшф) и переписать уравнение (20) на интервале ф = (0, л) в виде

Это уравнение по форме совпадает уравнением (5), (10), учитывая переобозначения

Б ^Бр2, р ^ )р, хед ^ агс8ш

В результате получаем приближенные собственные значения е, (р), которые зависят от р как от параметра,

(p) = p2ef, ef = -Ь

\_р(з - 2V

л/Ь

2

j = l,...,n

VD 1

p +2

Эти собственные значения имеют правильное асимптотическое поведение £ , (р) /р2 = ё ,, описывающее нижнюю часть дискретного спектра задачи (2). В рассмотренном случае они соответствуют первым пяти (^ = 1,..., 5) собственным значениям е 1,..., £5. Соответствующие собственные функции ф,(ф; р) при ] = 1, ...,п, параметрически зависящие от медленной переменной р через новую независимую

переменную С = ((ф; р) = 2р\[Ъехр[—рр(ф — агс8ш(жед/р))]/), С € [0, +го) принимают вид

ф,«;p) = N

где в, = \[~Б/) — ] + 1/2 — положительный параметр. В рассмотренном случае волновая функция вне потенциальной ямы при |х — хея|/р ^ 1 экспоненциально убывает. Этот факт позволяет интегрировать произведения функций ф,

и/или дфз (С(р;р);р)/др|^ =со^ по переменной £ на полубесконечном интервале С € (0, Вычисленные собственные функции при р =10 для ] = 1,..., 5, по-

казанные на рис. 4, также локализованы в потенциальной яме, что и функции связанных состояний, приведённые на рис. 2. Матричные элементы между состояниями нижней части дискретного спектра %,'} = 1,...,п = 5 с собственными

значениями е^ (р)/р2 ~ ^ представимы в виде асимптотических разложений по обратным степеням р:

vmax

J2fc)

в, (р) = ef У +

max

j Г 1 А^! р2к '

fc=l

п(2к-1) kmax Н(2к) V '

Qij

к=1 р к=1 р

и вычисляются до требуемого порядка fcmax согласно приведённому алгоритму, реализованному в системе компьютерной алгебры MAPLE. Например, вычисленные коэффициенты Q(1) и Н(2) разложений (30) представлены в табл. 1. Там же для сравнения приведены численные значения матричных элементов Qij(p) и Hij(p) при р = 100.

Таблица 1

Вычисленные коэффициенты Q(1) Н2 асимптотических разложений (30) (верхние строки) (в A-2) и соответствующие численные значения Qij и Hij

при р = 100

Q(j) Qij 1 2 3 4 5

1 0 0 10,320087 0,10 322049 -3,817896 -0,03 818264 1,831702 0,01 831845 -0,903311 -0,00 903906

2 -10,32008759 -0,10 322049 0 0 12,24187273 0,12 245108 -5,544019175 -0,05 545224 2,689862642 0,02 691957

3 3,817896 0,03 818264 -12,241872 -0,12 245108 0 0 11,509603 0,11 514688 -5,307448 -0,05 312602

4 -1,831702 -0,01 831845 5,544019 0,05 545224 11,509603 0,11 514688 0 0 7,994372 0,08 004568

5 0,903311 0,00 903906 -2,689862 -0,02 691957 5,307448 0,05 312602 -7,994372 -0,08 004571 0 0

Н(2) ij Hij 1 2 3 4 5

1 127,980 0,0128 021 -67,281 -0,0067 296 -85,381 -0,0085 433 73,391 0,0073 436 -44,488 -0,0044 545

2 -67,281 -0,0067 296 317,555 0,0317 675 -161,486 -0,0161 528 -40,519 -0,0040 598 46,871 0,0046 971

3 -85,381 -0,0085 433 -161,486 -0,0161 528 408,505 0,0408 705 -231,090 -0,0231 186 45,251 0,0045 224

4 73,391 0,0073 436 -40,519 -0,0040 598 -231,090 -0,0231 186 385,918 0,0386 221 -215,798 -0,0216 076

5 -44,488 -0,0044 545 46,871 0,0046 971 45,251 0,0045 224 -215,798 -0,0216 076 223,510 0,0224 050

Видно, что первые ненулевые коэффициенты этих разложений дают аппроксимацию численных значений матричных элементов с тремя значащими цифрами.

Коэффициенты Я^, Яр^, Н3 и Н^ разложений (30), вычислен-

ные с помощью интерполяции численных значений ез(р), Я^(р), Н^(р) при р = 50,100, 200, представлены в табл. 2.

Таблица 2

Коэффициенты Q<'^]\ , (верхние строки) и , Н(4\ (нижние строки) (в А-2), вычисленные с помощью интерполяции численных значений (р), (р), (р) при р = 50,100, 200

Я? я<3 ) 1 2 3 4 5

1 0 0 -10,32008 -19,61136 3,81789 3,63249 -1,83170 -1,59581 0,90330 6,37175

2 10,32008 19,61136 0 0 -12,24187 -32,23030 5,54401 12,57078 -2,68984 -22,23888

3 -3,81789 -3,63249 12,24187 32,23031 0 0 -11,50960 -52,07881 5,30742 54,16579

4 1,83170 1,59581 -5,54401 -12,57078 11,50960 52,07881 0 0 -7,99432 -105,8683

5 -0,90330 -6,37175 2,68984 22,23888 -5,30741 -54,16577 7,99432 105,8683 0 0

Н(2) г3 Н(4) 1 2 3 4 5

1 127,981 414,478 -67,282 -142,335 -85,381 -526,512 73,391 464,741 -44,488 -628,710

2 -67,282 -142,335 317,556 1204,009 -161,486 -420,383 -40,519 -812,200 46,871 1126,996

3 -85,381 -526,512 -161,486 -420,383 408,505 2009,349 -231,090 -948,323 45,253 -471,183

4 73,391 464,741 -40,519 -812,200 -231,090 -948,323 385,918 3046,143 -215,798 -2608,49

5 -44,488 -628,710 46,871 1126,996 45,253 -471,183 -215,798 -2608,49 223,510 5105,98

1 2 3 4 5

е(0) с 3 -193,066013 -119,392672 -63,338854 -24,904558 -4,089760

с3 -127,730560 -317,305630 -408,271216 -385,695053 -223,621354

Выполнение асимптотического соотношения е3 ) + Н^- = 1/4 соответствует правильному асимптотическому поведению решения в координатах (х, у) в 8-волновом приближении при отсутствии центробежных членов и служит критерием точности вычисления собственных функций задачи на собственные значения (20)—(21) и их производных по параметру р при его больших значениях.

Из табл. 2 видно, что для первых пяти состояний это соотношение при р =100 выполняется с точностью не хуже ~ 1 %. Однако с увеличением номера состояния

] абсолютная точность выполнения соотношения е3 + Н^ = 0, 25 ухудшается, и для ] = 1,..., 5 соответственно равна: 0,0000; 0,0000; 0,0161; 0,0267; 0,3715, что приводит к необходимости увеличения числа узлов конечноэлементной сетки на подынтервале р € (ср\(р),р2(р)).

Матричные элементы между состояниями с собственными значениями £з(Р ^ = — п)2 > 0, г,з = п +1 = 6, 7,... представимы в виде асимптотических разложений по обратным степеням р:

(Р) = U -п)2 + ^ ^ k=1 Р

kmax ,_

kmax Q 2k) km ax -

QM= E Q^, fii,-(p) = E , Vij(p) = E ik-^. k=l Р k=l Р k=l Р

Коэффициенты Q(2), Q(4), H(4) и H(6), и разложений (31), вычисленные с помощью интерполяции численных значений Qij(p), Hij(p) и Vij(p) при p = 50,100,200, представлены в табл. 3 и пригодны для экстраполяции с точностью < 5-10%.

Матричные элементы между состояниями с собственными значениями Ej (p) ^ 0, j = п + 1,... и £j(p)/p2 ~ ej0) < 0, j = 1,..., п = 5 представимы в виде асимптотических разложений по обратным степеням p:

kmax Q2k + 1/2) kmax —

Qij(p)= E ijk+1/2 , —ij(p)= E p2k+3/2 , Vij(p) = 0(exp(-p)). (Щ k=i p k=i p

Коэффициенты Q(.5/2), Q(j/2, Hij7/2) и —j11/2) разложений (32), вычисленные с помощью интерполяции численных значений Qi ( p), Hi ( p) при p = 50, 100, 200, представлены в табл. 4 и пригодны для экстраполяции с точностью ^ 1%. Корневое поведение матричных элементов объясняется тем, что амплитуда собственных функций дискретного спектра j = 1,...,п пропорциональна /p, а собственных функций дискретного спектра j = п + 1,..., выходит на константу.

На рис. 9 представлены волновые чётные и нечётные функции фj(p;p) краевой задачи (20), (21) c собственными значениями Sj(p) ^ 0 при j = п +1 = 6,..., 10. Положения нулей собственных функций, представленных на эпюре рис. 9, с увеличением их номера j > 6 незначительно смещаются в сторону р = 0.

Для состояний i,j = п +1, ..., Jmax с собственными значениями £j(p ^ 'X))/p2 = (j-п)2/p2 +0(1/p4) =s(0)/p2 + 0(1/p4) = k2 +0(1/p2), соответствующими псевдосостояниям краевой задачи (5), (6), сведённой на конечный интервал р £ (0, т/2), использовано приближение собственными функциями непрерывного спектра (см. уравнение (13) в обозначениях (28)) с помощью процедуры, реализованной в системе компьютерной алгебры MAPLE. Спектр энергии чётных и нечётных состояний вычислялся с учётом вышеуказанного поведения нулей собственных функций в окрестности минимума потенциала Морзе, соответственно подчинённых краевым условиям

d <j>k(p;p) | n i f ю \ n

dp |^=./2 =0 и <k (t/2;p) =

Поведение вычисленных собственных функций при p = 10 для = 6, ... , 10 согласуется с численными собственными функциями, показанными на рис. 4, и качественно согласуется с поведением собственных функций псевдосостояний, показанных на рис. 2. Таким образом, базисные функции разложения Галёркина (4) соответствуют асимптотическим базисным функциям разложения Канторовича (19) при больших значениях параметра p.

Таблица 3

Коэффициенты , Н(4\ (верхние строки) и , Н^6, (нижние

строки) (в А-2), вычисленные с помощью интерполяции численных значений Qij (р), (р), Уц (р) при р = 50,100, 200

9<?' ей' 6 7 8 9 10

6 0 0 -1,68690 -301,546 0,94773 -69,297 -0,65908 -61,462 0,50651 -46,423

7 1,69497 139,529 0 0 -4,26270 257,625 2,37296 218,332 -1,68318 -78,820

8 -0,94152 -79,974 4,22617 578,150 0 0 -6,59009 -640,149 3,59129 888,837

9 0,65906 58,134 -2,36875 -286,875 6,61712 -46,624 0 0 -8,77460 -3393,818

10 -0,50837 -49,683 1,69196 215,297 -3,64165 -100,580 8,89720 914,186 0 0

Я(4' г3 Я(6' 6 7 8 9 10

6 5,9167 13052,14 -10,8673 -38578,53 6,8113 61320,91 -6,1801 -85458,38 6,2733 111555,5

7 -10,8673 -38578,53 42,8126 120399,0 -47,3459 -186997,1 24,7308 258580,3 -21,2839 -334667,7

8 6,8113 61320,91 -47,3459 -186997,1 118,1730 299064,5 -107,5731 -432192,6 50,7212 553049,6

9 10 -6,1801 -85458,38 6,2733 111555,5 24,7308 258580,3 -21,2839 -334667,7 -107,5731 -432192,6 50,7212 553049,6 228,5693 634226,5 -189,1366 -833872,0 -189,1366 -833872,0 368,9474 1220438,

уА1' г з уА3' 6 7 8 9 10

6 121,5746 -13414,11 -124,27405 -53281,30 -118,94685 89202,5 126,9529 114968,4 116,5646 -175109,8

7 -124,2740 -53281,30 126,9742 126302,8 121,6545 -29390,91 -129,6553 -194078,6 -119,2950 122761,6

8 -118,9468 89202,5 121,6545 -29390,91 116,3012 -157521,5 -124,3399 -25644,41 -113,8856 235349,4

9 126,9529 114968,4 -129,6553 -194078,6 -124,3399 -25644,41 132,3411 267673,0 122,0015 -74954,4

10 116,5646 -175109,8 -119,2950 122761,6 -113,8856 235349,4 122,0015 -74954,4 111,4203 -304457,9

Диагональные и недиагональные барьерные матричные элементы У^(р) показаны на рис. 6 и 8. Их следует сравнивать с соответствующими матричными элементами, показанными на рис. 2. Из сравнения следует, что матричные элементы Уг](р) из (24) между состояниями дискретного спектра краевой задачи (20), (21) и матричные элементы У^(у) из (8) между состояниями дискретного спектра и псевдосостояниями (5), (6) демонстрируют качественное сходное поведение

Таблица 4

Коэффициенты н\у2 (верхние строки) и Q(iУ2\ (нижние

строки) асимптотического разложения (в А-2), вычисленные с помощью интерполяции численных значений Qij (р), Нг^ (р) при р = 50,100, 200

<$/2) 6 7 8 9 10

1 -0,61311 1,83562 -3,07274 4,29103 -5,48095

-191,631 640,795 -747,558 1270,066 -2354,089

2 1,82101 -5,45193 9,12628 -12,74469 16,27883

571,136 -1909,166 2229,887 -3785,644 7009,516

3 -3,57012 10,68864 -17,89229 24,98625 -31,91492

-1125,595 3760,664 -4399,770 7461,106 -13794,19

4 5,24449 -15,70152 26,28377 -36,70467 46,88247

1665,144 -5559,590 6515,254 -11032,27 20360,29

5 -5,85899 17,54104 -29,36412 41,00545 -52,37294

-1842,771 6157,353 -7127,912 12085,368 -22431,565

Н (7/2) 1 Ч Я(11/2) 1 «7 6 7 8 9 10

1 32,2121 -96,4400 161,4380 -225,4437 287,9541

10571,20 -35184,13 41714,53 -70159,27 128212,5

2 -38,0069 113,7880 -190,4818 265,9993 -339,7444

-13296,49 43998,75 -53208,53 88399,05 -158710,1

3 -17,4025 52,1043 -87,2089 121,7953 -155,6006

-3104,68 11140,26 -9918,36 20120,33 -45697,38

4 121,5501 -363,9136 609,1614 -850,6939 1086,628

35804,24 -120419,4 137603,1 -236753,9 446549,4

5 -172,1671 515,4534 -862,8458 1204,948 -1539,080

-53709,59 179539,0 -208786,7 354333,5 -656920,8

Рис. 9. Волновые функции ф^ (<р; р): чётные — левая панель, нечётные —

правая панель, краевой задачи (20), (21) с собственными значениями £з (р) ^ 0 при ] = п +1 = 6,..., 10 и значении параметра р = 100. На нижних панелях даны эпюры на интервале 6 [0, 02; 0, 05] в окрестности минимума

потенциала Морзе

по переменным у и р. Поскольку р = \/х2 + у2 > у, то, очевидно, потенциалы Уг^(р) сильно делокализованы по отношению к потенциалам Уг^(у).

Благодаря медленно убывающему кинематическому поведению потенциалов 9(р) и Нц(р) как р-1 и р-2 соответственно, по сравнению с экспоненциальным убыванием Уц(у), необходимо учитывать лидирующие члены разложения их асимптотических разложений при решении краевых задач (22)—(24), генерируемых разложением Канторовича (17) при вычислении искомых величин задачи рассеяния с пятью открытыми каналами даже при энергиях ниже порога диссоциации.

5. Алгоритм вычисления асимптотик фундаментальных

решений

Алгоритм 2. Вычисляем асимптотическое решение системы N обыкновенных дифференциальных уравнений второго порядка при больших значениях независимой переменной р ^ 1

— —У—. ++ Кгг(р) - 2Е

Хгг '

1±р А + £Ар)

р Пр Пр р2

N

п I п

Хц'

Е

3 = 1, 3=г

-9гз

Входные параметры. Коэффициенты уравнений (33), где = У-^ + Нц представлены в виде разложений по обратным степеням (30). ТТТяг 1. Ищем решение уравнений (33) в виде:

Хзг'

где ф^г'(р) и ^^г /(р) — неизвестные функции, Яг' (р) — известная функция. Выбираем Яг' (р) как решение вспомогательной задачи для эталонного уравнения:

1 п п , ^А2' 2

р п р п р р2

яг'

где 7(2' = р(2' + Ъ(2'

где 7г = £г' + нг' г'.

Шаг 2. Вычисляем коэффициенты фг/ (р) и фг> (р) разложения (34) в форме усечённого разложения по обратным степеням р (ф^к'<° = ГФ%<=0):

^шах ф(к'' ^шах „/,(к''

Фзг

к'=1 р к'=1 р

После подстановки выражений (34)—(36) в уравнение (33), используя уравнение (35), получаем систему рекуррентных соотношений:

(4°' - 2Е + р2 ) ф%' - 2р2 (к' - 1)$-1' = - £',

г(к' (к'

где правые части г и г даются соотношениями

№) = (-( - 2)2 - 2) + Е -?фГ k)+

\2 у

k=2

k' W

+ z(2)(2k' - -3) + £ £ (ZQ^Z^-k-2)-

k=1 j=1,j=i

- 2^Q(kVjk'-k) + Q(k)(-2k' + k + 3)0jk-k-1) + -k));

= (-( kk -1)2 - zp^c;-2) +^ -k)+

7

k=2

W k'

+ E E

j=1, j=i k=1

с начальными условиями = 2 Е - е, ф((/ = fe/, ф((

0.

ТТТаг 3. Из рекуррентных соотношений (37) следует рекуррентная формула

4k ), при 2 Е = б(,0), % = i' и k' = 2,. . . , kmax:

i-') + 2р2( k' - 1)ф£ '-1); ,

- (k) - 2( k' - 1)Й '-1);, Ф&-1) = - [2(k' - 1)]-1 д§, -1) = l~2(k' - 1) Г2Е - ef

оэффициентов ф

) = Гр(0) „(0)1 е i - е ('

) = |р(0) р(0) £( - £с

-1

1

Вышеописанный алгоритм, реализованный в системе компьютерной алгебры MAPLE, генерирует подпрограммы на языке FORTRAN для вычисления искомых

коэффициентов разложения ф^,) и фЦ?,) до требуемого порядка kmax.

Выбор подходящего значения pmin и pmax для генерации разложений линейно независимых решений при р(о > 0 контролировался выполнением условия на вронскиан с требуемой точностью ewr:

2

Wr(Q(p); Х*(Р), Х(Р)) = ~Iоо,

ж

w (Q, X*, х) = х*т ( - Qx) - хт( - Qx*

dp

6. Анализ модели квантового туннелирования

Решение краевых задач (7)-(14) и (22)—(26) выполнялось на конечно-элементных сетках Пу = {-12(^егет = 120)12} и Пр = {0(^егет = 1200)120} соответственно, с Ме1ет лагранжевыми элементами четвёртого порядка р = 4 между узлами, используя фортрановскую программу КАМТБР 3.0 [12]. Разложения (4) и (19) искомого решения по построенным ортогональным базисам при (]тах = 15) только с десятью закрытыми каналами позволило вычислять приближенные решения двумерных краевых задач (1) и (15) при Е < 0 с точностью (10-4).

Рис. 10 демонстрирует резонансное поведение полных вероятностей прохождения барьеров из первых каналов, имеющих энергии Е1 = -1044,879649, Е2 =

-646,1570935, Е3 = -342, 7919791, Е4 = -134, 7843058, Е5 = -22,13407384 (в К), с переходами во все открытые каналы, вычисленные с помощью разложений (4) и (17) в декартовых и полярных координатах.

Рис. 10. Полные вероятности прохождения барьеров из открытых каналов

Полные вероятности прохождения барьеров из открытых каналов демонстрирует резонансное поведение, т.е. эффект квантовой прозрачности. С ростом энергии начальных возбуждённых состояний пики смещены в сторону больших энергий, набор положений пиков примерно тот же, что для переходов из основного состояния. Например, левый эпюр показывает, что положения 13-го и 14-го пиков из первого состояния совпадают с положениями 1-го и 2-го пиков из второго состояния, тогда как правый эпюр показывает, что положения 25-го и 26-го пиков из первого состояния совпадает с положениями 1-го и 2-го пиков из третьего состояния.

Из рис. 2 видно, что диагональные матричные элементы потенциала У^. (у) имеют форму двойных барьеров, а недиагональные матричные элементы У^. (у) примерно в четыре раза меньше, чем У^.(р) и У^-(р) на рис. 6 и 8. Это означает, что положение пиков соответствует вещественной части энергии метастабильных состояний, погруженных в непрерывный спектр, которые локализованы между двойными барьерами.

7. Заключение

Продемонстрирована эффективность разработанных символьно-численных алгоритмов и комплексов программ, реализующих метод конечных элементов для решения краевых задач, описывающих модели квантового туннелирования двухатомных молекул в 8-волновом приближении, связанных реалистическими молекулярными потенциалами ниже порога диссоциации. Представлен сравнительный анализ базисных функций, матриц переменных коэффициентов и решений систем уравнений в форме Галёркина и Канторовича в декартовых и полярных системах координат. Выявлен эффект квантовой прозрачности при резонансном туннели-ровании двухатомных молекул через отталкивающие барьеры за счёт метаста-бильных состояний, погруженных в непрерывный спектр. Предложенная модель и разработанные алгоритмы и комплексы программ применимы для описания фрагментации лёгких ядер, подбарьерных реакций ионов и молекулярной квантовой диффузии.

Литература

1. Пеньков Ф. М. Квантовая прозрачность барьеров для структурных частиц // ЖЭТФ. — 2000. — Т. 118. — С. 806-815.

2. Pen'kov F. M. Metastable States of a Coupled Pair on a Repulsive Barrier // Physical Review A. — 2000. — Vol. 62. — P. 044701.

3. Pijper E., Fasolino A. Quantum Surface Diffusion of Vibrationally Excited Molecular Dimers // Journal of Chemical Physics. — 2007. — Vol. 126. — P. 014708.

4. Kavka J. J., Shegelski M. R. A., Hong W. P. Tunneling and Reflection of an Exciton Incident Upon a Quantum Heterostructure Barrier // Journal of Physics: Condensed Matter. — 2012. — Vol. 24. — P. 365802.

5. Shegelski M. R. A., Hnybida J., Vogt R. Formation of a Molecule by Atoms Incident Upon an External Potential // Physical Review A. — 2007. — Vol. 78. — P. 062703.

6. Bondar D. I., Liu W.-K., Ivanov M. Y. Enhancement and Suppression of Tunneling by Controlling Symmetries of a Potential Barrier // Physical Review A. — 2010. — Vol. 82. — P. 052112.

7. Symbolic-Numerical Algorithm for Generating Cluster Eigenfunctions: Tunneling of Clusters Through Repulsive Barriers / A. Gusev, S. Vinitsky, O. Chuluunbaatar et al. // Lecture Notes in Computer Science. — 2013. — Vol. 8136. — Pp. 427442.

8. 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.

9. Гусев А. А. Модель туннелирования кластеров через отталкивающие барьеры в представлении симметризованных координат // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2014. — № 1. — С. 54-73.

10. 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 // Computer Physics Communications. — 2009. — Vol. 180. — Pp. 1358-1375.

11. Gusev A. A. Algorithm for Computing Wave Functions, Reflection and Transmission Matrices of the Multichannel Scattering Problem in the Adiabatic Representation using the Finite Element Method // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2014. — № 2. — С. 93-114.

12. KANTBP 3.0: New Version of a Program for Computing Energy Levels, Reflection and Transmission Matrices, and Corresponding Wave Functions in the CoupledChannel Adiabatic Approach / A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, A. G. Abrashkevich // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2014. — № 2. — С. 342-349.

13. 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 // Computer Physics Communications. — 2008. — Vol. 179. — Pp. 685-693.

14. Fundamental Physical Constants. — http://physics.nist.gov/constants.

15. Lauhon L. J., Ho W. Direct Observation of the Quantum Tunneling of Single Hydrogen Atoms with a Scanning Tunneling Microscope // Physical Review Letters. — 2000. — Vol. 85. — Pp. 4566-4569.

16. 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. // Journal of Physics A. — 2007. — Vol. 40. — Pp. 1148511524.

UDC 517.958:530.145.6 Algorithm for Solving the Two-Dimensional Boundary Value Problem for Model of Quantum Tunneling of a Diatomic Molecule Through Repulsive Barriers

A. A. Gusev*, L. L. Hait

* Laboratory of Information Technologies Joint Institute for Nuclear Research 66, Joliot-Curie str., Dubna, Moscow region, Russia, 141980

^ Belgorod State National Research University 85, Pobedy str., Belgorod, Russia, 308015

Algorithm for solving the boundary value problems that describe the model of quantum tunneling of a diatomic molecule through repulsive barriers in s-wave approximation is presented. The boundary value problems are formulated and reduced to the one-dimensional ones for systems of coupled second-order differential equations by means of the Galerkin and Kantorovich methods. The description of elaborated algorithms and the calculated asymptotes of parametric basis functions, matrices of variable coefficients, and fundamental solutions of the systems of the coupled second-order differential equations needed for solving the boundary problems on a finite interval are given. The BVPs were solved by the elaborated set of programs implementing the finite element method. Analysis of benchmark calculations of quantum tunneling of a diatomic molecule model with the nuclei coupled by the Morse potential through Gaussian barriers and quantum transparency effect induced by metastable states embedded in continuous spectrum below dissociation threshold are presented.

Key words and phrases: quantum tunneling problem, diatomic molecule, repulsive barriers, boundary-value problems, Galerkin method, Kantorovich method, asymptotic solutions, finite element method.

References

1. F. Pen'kov, Quantum transmittance of barriers for composite particles, ZHETF 91 (2000) 698-705, in Russian.

2. F. M. Pen'kov, Metastable States of a Coupled Pair on a Repulsive Barrier, Physical Review A 62 (2000) 044701.

3. E. Pijper, A. Fasolino, Quantum Surface Diffusion of Vibrationally Excited Molecular Dimers, Journal of Chemical Physics 126 (2007) 014708.

4. J. J. Kavka, M. R. A. Shegelski, W. P. Hong, Tunneling and Reflection of an Exciton Incident Upon a Quantum Heterostructure Barrier, Journal of Physics: Condensed Matter 24 (2012) 365802.

5. M. R. A. Shegelski, J. Hnybida, R. Vogt, Formation of a Molecule by Atoms Incident Upon an External Potential, Physical Review A 78 (2007) 062703.

6. D. I. Bondar, W.-K. Liu, M. Y. Ivanov, Enhancement and Suppression of Tunneling by Controlling Symmetries of a Potential Barrier, Physical Review A 82 (2010) 052112.

7. A. Gusev, S. Vinitsky, O. Chuluunbaatar, V. A. Rostovtsev, L. L. Hai, V. Derbov, P. Krassovitskiy, Symbolic-Numerical Algorithm for Generating Cluster Eigenfunctions: Tunneling of Clusters Through Repulsive Barriers, Lecture Notes in Computer Science 8136 (2013) 427-442.

8. A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar, V. P. Gerdt, V. A. Rostovtsev, Symbolic-Numerical Algorithms to Solve the Quantum Tunneling Problem for a Coupled Pair of Ions, Lecture Notes in Computer Science 6885 (2011) 175-191.

9. A. Gusev, A. Model of tunneling of clusters through repulsive barriers in symmetrized coordinates representation, Bulletin of PFUR. Series "Mathematics. Information Sciences. Physics" (1) (2014) 54-72, in Russian.

10. 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.

11. A. A. Gusev, Algorithm for Computing Wave Functions, Reflection and Transmission Matrices of the Multichannel Scattering Problem in the Adiabatic Representation using the Finite Element Method, Bulletin of PFUR. Series "Mathematics. Information Sciences. Physics" (2) (2014) 93-114.

12. A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, A. G. Abrashkevich, KANTBP 3.0: New Version of a Program for Computing Energy Levels, Reflection and Transmission Matrices, and Corresponding Wave Functions in the CoupledChannel Adiabatic Approach, Bulletin of PFUR. Series "Mathematics. Information Sciences. Physics" (2) (2014) 342-349.

13. 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, Computer Physics Communications 179 (2008) 685-693.

14. Fundamental Physical Constants.

URL http://physics.nist.gov/constants

15. L. J. Lauhon, W. Ho, Direct Observation of the Quantum Tunneling of Single Hydrogen Atoms with a Scanning Tunneling Microscope, Physical Review Letters 85 (2000) 4566-4569.

16. O. Chuluunbaatar, A. A. Gusev, V. L. Derbov, M. S. Kaschiev, L. A. Melnikov, V. V. Serov, S. I. Vinitsky, Calculation of a Hydrogen Atom Photoionization in a Strong Magnetic Field by Using the Angular Oblate Spheroidal Functions, Journal of Physics A 40 (2007) 11485-11524.