Применение гауссовых волновых пакетов для расчетов квантовомеханических систем тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат физико-математических наук Айдагулов, Галлям Раильевич

• Айдагулов, Галлям Раильевич
• кандидат физико-математических наук
• 2004

Оглавление диссертации кандидат физико-математических наук Айдагулов, Галлям Раильевич

Введение

1 Метод гауссовых волновых пакетов

1.1 Когерентные состояния.

1.2 Метод гауссовых волновых пакетов.

1.2.1 Квадратичный потенциал.

1.2.2 Потенциал: общий случай.

1.2.3 Модификация метода.

1.3 Вычисление наблюдаемых.

2 Построение численной процедуры

2.1 Предварительные замечания.

2.2 Доказательство теоремы

2.3 Численный алгоритм.

2.3.1 Фиксированная сетка.

2.3.2 Подвижная сетка.

2.4 Связь с методом частиц.

3 Результаты расчетов

3.1 Аналитическое решение

3.2 Гармонический осциллятор.

3.3 Равноускоренное движение частицы.

3.4 Двухъямный потенциал.

3.5 Потенциал Морзе.

3.6 Рассеяние на потенциальном барьере.

Введение диссертации (часть автореферата) на тему «Применение гауссовых волновых пакетов для расчетов квантовомеханических систем»

Предметом настоящей работы является разработка приближенного метода решения задачи Коши для нестационарного уравнения Шрёдингера дтЬ h2 ih-^(x,t) = -—AxP(x,t) + V(x)iP(x,t), teR]., (1) ф(х,0 + 0) = ф0{х), x e (2) где ip(-,t) E L2(Rd,dx), t G R+. В работе предлагается новый численный метод решения задачи. Также проводится ряд численных экспериментов с использованием нового метода, на основе результатов которых демонстрируется его практическая применимость.

Уравнение Шрёдингера (1) — одно из основных уравнений математической физики, описывающее эволюцию квантовой системы во времени. Из аксиом квантовой механики [9, 12] следует, что волновая функция системы, находящейся во внешнем поле с потенциалом V = V(x), есть решение задачи Коши (1)-(2), где — волновая функция в начальный момент времени. Задача (1)-(2) возникает, например, при теоретическом описании явлений в наноэлектронике, полупроводниках, квантовой химии и физике твердого тела. В общем случае уравнение (1) может быть решено только численно. В виду интенсивного развития математического моделирования в этих областях возникает проблема разработки эффективных численных методов решения задачи (1)-(2).

Обзор существующих методов

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

15 4 более употребительные в настоящее время подходы к решению задачи (1)-(2). В контексте настоящего исследования это удобно сделать следующим образом.

Постоянная Планка h в уравнении (1) играет фундаментальную роль во всех квантовых явлениях. Ее относительная величина (по сравнению с другими величинами той же размерности) определяет "степень квантовости" той или иной физической системы. В настоящей работе под h мы будем понимать параметр, входящий в уравнение (1), которое было обезразмерено при выборе системы единиц, наиболее всего отвечающей масштабу рассматриваемого физического явления. При этом формально квантовая механика содержит в себе классическую в качестве предельного случая при h —* 0. Например, при описании движения электрона вокруг ядра, где характерное расстояние имеет порядок Ю-10 м, в соответствующей (атомной) системе единиц величина h равна 1. В задачах физики твердого тела при длинах порядка 10~7 м, h имеет порядок 10-2. С дальнейшим уменьшением h квантовые эффекты проявляются все слабее, и система начинает подчиняться классическим законам. Случай малого значения параметра h (h -С 1) называется квазиклассическим [12] и рассматривается особо при изучении свойств задачи (1)-(2).

Оказывается, что и численные методы решения задачи (1)-(2) могут быть разделены в зависимости от величины h на квазиклассические (h «С 1) и квантовые (h ~ 1). При этом метод одного класса часто оказывается неэффективным в применении к задачам другого. Поясним это, кратко рассмотрев существующие методы.

Квазиклассические методы (случай h 1)

В основе всех квазиклассических методов лежит тот факт, что при /г —* 0 динамика рассматриваемой системы становится близка к классической. Эта информация позволяет свести исходную задачу (1)-(2) к более простым. В результате получаются приближенные методы, требующие существенно меньших ресурсов для получения решения, чем квантовые методы (см. ниже), основанные на "прямом" решении исходной задачи.

К квазиклассическим методам относится, прежде всего, известный метод ВКБ1, в котором решение задачи ищется с помощью разложения по степеням малого па

1 Метод назван по имени ученых - физиков Венцеля, Крамерса и Бриллюэна (Wentzel, Kramers, Brillouin, WKB). раметра h [9, 33]. Остальные методы, как правило, заключаются в рассмотрении классической системы, в которую в качестве поправок вносятся основные квантовые эффекты. Наиболее разработанным в этой области является метод гауссовых волновых пакетов (Gaussian wave packet dynamics, GWPD) [32, 45, 34]. В нем решение в любой момент времени предполагается в форме гауссиана, центр которого подчиняется классическим уравнениям движения. Такой выбор формы решения обусловлен тем, что гауссовы волновые пакеты, как состояния, в которых достигается минимум соотношения неопределенностей между координатой и импульсом, являются наиболее близкими квантовыми аналогами классического понятия точки в фазовом пространстве. Строгим математическим изложением этой техники могут служить работы [30, 31], где построены асимптотические разложения решения задачи (1)-(2) по базису гауссовых волновых пакетов при h—* 0.

Основным недостатком квазиклассических методов является именно их асимптотический характер. Оценки точности, как правило, имеют вид [31] O(hm)n-*o, т > 0. Это не позволяет получить решение с произвольной точностью в практических расчетах, в которых значение h хоть и мало, но все же фиксировано. Поэтому ф границы применимости получаемых приближений, вообще говоря, неясны, и требуют дополнительного исследования в каждом конкретном случае. Кроме того, по мере приближения к квантовому режиму (/г ~ 1) погрешности проявляются сильнее, что делает метод неприменимым в этом случае.

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

Для получения более подробной информации по квазиклассическому приближению и методу гауссовых волновых пакетов можно порекомендовать работы [42, 31, 34], обзор [33], а также цитируемую там литературу.

Квантовый случай (h ~ 1)

В квантовом случае (h ~ 1) наиболее употребительными являются сеточные методы решения (см. обзор Kosloff [43] и ссылки там), когда пространственная и временная части уравнения (1) аппроксимируются на сетке. Это представляется вполне естественным, так как величина постоянной Планка в данном случае определяет такой масштаб рассматриваемых явлений, что, вообще говоря, необходимо "использовать" всю информацию, которую дает задача (1)-(2), для адекватного описания процесса. Например, непригодность в квантовом случае квазиклассических методов, может быть обусловлена именно тем, что с самого начала не принимаются во внимание некоторые свойства решения, которыми можно пренебречь на макроуровне (h 1), но которые начинают играть существенную роль при переходе к квантовому режиму (h ~ 1).

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

Для удобства будем считать, что уже сделана аппроксимация пространственной части на первом этапе, и рассмотрим сначала временную составляющую. Пусть оператор Tid - некоторый сеточный аналог2 гамильтониана3 Н в правой части уравнения (1). Тогда переход на следующий временной слой формально определяется с помощью аппроксимации оператора эволюции

Здесь можно отметить следующие подходы

2Индекс d здесь означает именно дискретизацию, и, естественно, никак не связан с размерностью пространства в записи задачи (1)-(2).

3 Оператор правой части в уравнении (1)

К2 называется оператором Гамильтона и выделяется особо. В квантовой механике этот оператор соответствует величине полной энергии системы.

3)

• Двуслойная схема Кранка-Николсона (Crank-Nicholson, CN). Эта схема может быть интерпретирована как Паде-аппроксимация экспоненты (3) где Та — единичный оператор, действующий в пространстве сеточных функций. Если оператор 7id - самосопряжен, то схема Кранка-Николсона - эрмитова, абсолютно устойчива, сохраняет L2—норму решения. Основным недостатком схемы является ее неявность.

• Центральная разность (second order differencing, SOD). SOD-схема получается, если производную по времени в левой части (1) аппроксимировать центральной разностью. Получается явная трехслойная схема того же порядка по времени, что и CN. Но она является условно устойчивой и лишь приближенно унитарной.

• Схема расщепления (split operator technique, SP). При построении схем расщепления в применении к уравнению Шрёдингера4, как правило, рассматривают приближенную факторизацию оператора эволюции (3) с помощью различных вариаций формулы Троттера. Это хоть и не дает законченного алгоритма перехода на следующий временной слой, но позволяет свести исходную задачу к цепочке более простых. Впервые этот подход был предложен в работе Feit [26], где выделялись кинетическая и потенциальная составляющие гамильтониана К. = и У = V(x): П = К, + V. Тогда exp {-itn/h} = ехр{—0.5 itIC/h} exp{—itV/h} ехр{—0.5 itIC/h} + 0(t3), (4) что позволяет свести решение исходной задачи к последовательному решению двух стандартных подзадач: свободное двио/сение и локальное взаимодействие волны частицы с потенциалом. Забегая вперед, скажем, что каждая из них может быть решена точно в пространственном или импульсном представлении соответственно. Необходимые переключения между представлениями осуществляются с помощью дискретного преобразования Фурье, использование же алгоритма быстрого преобразования Фурье (БПФ) позволяет сделать эти переходы чрезвычайно эффективными и точными. Это делает метод расщепления

4Общее изложение метода расщепления при построении численных алгоритмов можно найти t + o(t% в [13].

Feit одним из наиболее часто используемых и в настоящее время, см., например, [28]. Также с его помощью получено большинство результатов численного моделирования в книге [8].

В качестве модификации метода можно отметить работу [20], где вместо (4) используется более точная факторизация, состоящая из семи множителей и имеющая погрешность порядка 0(t4). Кроме того, выбор того или иного способа расщепления гамильтониана определяется спецификой задачи и предполагаемого дальнейшего способа решения получаемых подзадач. В [24] требовалось найти пропагатор задачи в представлении когерентных состояний гармонического осциллятора. Поэтому в исходном гамильтониане выделялась гармоническая и негармоническая части. В этом случае эволюция когерентного состояния под действием гармонического потенциала столь же тривиальна, как и свободное движение плоской волны в стандартном методе расщепления Feit (см. главу 1 и метод гауссовых волновых пакетов).

• Глобальный пропагатор (global propagator). Все рассмотренные до сих пор приближения оператора эволюции (3) строились за счет выбора достаточно малого шага по времени и могут быть получены из ряда Тейлора

Можно же рассмотреть задачу о построении так называемого глобального про-пагатора, который бы позволил находить приближенное решение задачи в любой, не обязательно малый, момент времени t "за один шаг". Отмечается, что использование с этой целью разложения Тейлора (5) оказывается неустойчивым. Поэтому в [53] был предложен метод, в котором для аппроксимации оператора эволюции используется отрезок ряда по системе полиномов Чебышева (СН) вида

Коэффициенты ап выражаются через значения функций Бесселя, а необходимые при реализации алгоритма выражения Рп (—itTid/h) фо, n = 0,l,.,N могут быть найдены последовательно в ходе вычислений по известным рекуррентным формулам для многочленов Чебышева. В виду быстрой сходимости ряда, пред

5) п=0

6) ложенный метод имеет высокую точность, которая определяется выбором не шага по времени, а количеством членов удержанных в разложении (6).

Одним из принципиальных недостатков глобального пропагатора является то, что он не дает возможности находить решение задачи в промежуточные моменты времени. Но такой подход оказывается чрезвычайно удобным когда требуется определить решение "на больших временах". Например, при нахождении основного состояния системы.

Рассмотрим теперь методы аппроксимации пространственной части. Комбинируя их с представленными выше методами перехода на следующий слой по времени, можно получить конкретные численные алгоритмы решения задачи (1)-(2).

• Конечные разности (finite differencing, FD). Считается, что численное решение нестационарного уравнения Шрёдингера началось именно с метода конечных разностей. Одной из первых в этой области считается работа [48]. В ней для аппроксимации пространственной производной в правой части уравнения (1) используется вторая разностная производная, а шаг по времени делается с помощью рассмотренной выше схемы Кранка-Николсона. Таким образом получается двуслойная разностная схема, которая является абсолютно устойчивой, сохраняет норму решения и имеет второй порядок аппроксимации по пространству и времени. В отечественной литературе обсуждение разностных схем для уравнения Шрёдингера можно найти в книге [15]. Здесь рассматривается общий случай двуслойной схемы с весами и показывается, что устойчивые схемы есть только среди неявных, причем норму решения сохраняет только схема Кранка-Николсона. По этой причине эта схема в течение некоторого времени была наиболее часто используемым методом. Необходимость использовать неявные схемы усложняет решение задачи в многомерном случае. Тем не менее такие расчеты проводились. В работе [46] двумерный лапласиан аппроксимируется на тринадцатиточечном шаблоне, а возникающая при этом система уравнений решается итерационным методом. Один из способов сократить объем вычислений - использовать экономичные методы решения [15]. Например, в [38, 50] для решения задачи в двумерном случае используется схема переменных направлений (Peaceman-Rachford), которая может быть рассмотрена как приближенная факторизация (расщепление) схемы Кранка-Николсона.

Трехслойная явная устойчивая схема предложена в работе [19]. Для аппроксимации по пространству здесь используется, как и выше, вторая разность, а по времени - центральная, т.е. SOD-пропагатор. Эта схема оказывается условно устойчивой, что накладывает жесткие ограничения на возможные шаги сетки. Кроме того, она является лишь приближенно унитарной (общее свойство SOD). Но можно показать, что отклонение нормы решения от 1, на самом деле, на несколько порядков меньше погрешности самой схемы и этот недостаток незначителен, в сравнении с возможностью использования в многомерном случае [19]. Заметим, кстати, что в применении к уравнению теплопроводности такая схема оказывается неустойчивой.

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

В применении к уравнению Шрёдингера наиболее употребительным является разложение по тригонометрической системе и для аппроксимации пространственной производной используется связь между образами Фурье самой функции и ее производной. Тогда оператор Лапласа Аф(х), х 6 Rd, может быть вычислен с помощью последовательного нахождения (/-мерного преобразования Фурье, последующего умножения полученного результата на функцию —(Щ + к% + . + к%), к = (ki,.,kd) 6 Rd, и вычисления обратного Фурье преобразования окончательного выражения. Дискретный аналог преобразования Фурье позволяет, применяя эту процедуру к вектору узловых значений функции ф, получить сеточную аппроксимацию6 лапласиана функции в правой

5Который в свою очередь можно рассматривать как частный случай метода коллокаций, см. [42,44].

6В [7, 6, 11] дискретное преобразование Фурье рассматривается в контексте тригонометрической интерполяции. Последняя тем точнее, чем глаже функция. Это замечание позволяет рассматривать предлагаемый способ аппроксимации производной как еще один вариант использования интерполирования при построении формул численного дифференцирования. Для подробностей Kosloff отсылает к [23]. Также см. ссылки по спектральному методу ниже. части уравнения (1) и, как результат, приближенный аналог всего гамильтониана. Использование алгоритма БПФ позволяет существенно уменьшить объем вычислений.

К самому раннему использованию описанного метода для решения уравнения Шрёдингера можно отнести работу Feit [26], где с его помощью, по сути, решалась возникающая в ходе расщепления подзадача о свободном движении пакета. Обобщение на случай потенциала общего вида было предложено Kosloff в [39] и успешно опробовано при моделировании конкретной реакции в [40]. Здесь для аппроксимации пространственной части использовалась описанная процедура, а шаг по времени делался по явной SOD-схеме. Следствием последнего является приближенная унитарность и условная устойчивость.

В [39] указывается основное ограничение на использование метода: решение задачи должно быть "представимо" на сетке. Это означает, что волновой спектр решения должен лежать внутри полосы, ширина которой определяется шагом пространственной сетки

1*1 < *с = jr. (7)

Природа этого ограничения та же самая что и у частоты Найквиста в теории обработки сигналов. Хотя условие (7) представляет собой серьезное ограничение на гладкость решения, но на практике оно является следствием ограниченности диапазона возможных энергий пакета и, поэтому, его можно считать выполненным, выбирая достаточно малый пространственный шаг7 h. С другой стороны, если параметры выбраны правильно и решение "представимо", то метод Фурье аппроксимации производной оказывается гораздо точнее разностных отношений. Следствием этого является возможность использования более разреженных сеток по сравнению с разностными схемами (до пяти раз в каждом направлении [40]), что очень важно при рассмотрении многомерного случая. Другим важным преимуществом спектрального метода является возможность устранения численного нефизичного уширения пакета, которое всегда возникает при использовании конечных разностей.

Представленный метод может быть рассмотрен с другими пропагаторами. В [41,

Пример подобного алгоритма, понятие частоты Найквиста, а также краткий обзор теории дискретного преобразования Фурье, можно найти в [51, р.490].

49] вместо SOD в расчетах используется глобальный СН-пропагатор [53]. Причем в [41] непосредственно используется "глобальный" характер аппроксимации: ведется поиск основного состояния - т.е. как раз случай, когда решение в промежуточные моменты времени не важно.

В заключение отметим, что вместо тригонометрической системы функций может быть использована любая другая система, наиболее подходящая в каждом конкретном случае. Так, в работе [22] рассматривается уравнение Шрёдингера в полярной (сферической) системе координат. Для вычисления радиальной составляющей лапласиана в этом случае удобнее использовать специальное преобразование Ханкеля, для эффективного вычисления которого здесь же и был разработан аналог алгоритма БПФ. Также, для аппроксимации производных может быть использовано разложение по системе многочленов Чебышева [25,36]. Этот способ имеет преимущества при рассмотрении непериодических функций.

• Finite Basis Representation (FBR). На примере уже изложенных методов видно, что для реализации того или иного способа перехода на следующий временной слой всегда возникает необходимость вычисления результата действия оператора гамильтона на решение задачи. В связи с этим вполне естественно строить пространственные дискретизации волновой функции и оператора таким образом, чтобы выполнение этой операции было максимально эффективным. Эта задача может быть решена за счет выбора специального представления (базиса в пространстве состояний), в котором оператор имеет наиболее простой вид. В иностранной литературе этот подход часто называется как finite basis representation (FBR) или discrete variable representation (DVR) technique. Изложенный выше псевдоспектральный метод Kosloff полностью укладывается в эту схему, так как по сути, используемая в нем тригонометрическая система является системой собственных функций гамильтониана (свободной частицы), который в ее представлении имеет диагональный вид.

Примерами работ в этом направлении могут являться [37, 35] (и ссылки там). В [37] волновая функция представляется в виде разложения по базису состоящему из произведений собственных функций гармонического осциллятора. Выбор такого представления был обусловлен физическими особенностями моделируемого явления. Учет специфики задачи позволил использовать меньше членов разложения, чем в общем методе Фурье. Это в свою очередь позволило рассмотреть системы с большим количеством степеней свободы. Для аппроксимации по времени использовалась модификация SOD-схемы четвертого порядка за счет удержания большего числа членов в разложении (5) и использования квадрата гамильтониана.

Итак, мы представили краткий перечень наиболее часто употребляемых методов решения задачи (1)-(2) в квантовом случае. Для получения более подробной информации можно порекомендовать читателю обратиться к процитированным выше работам и ссылкам в них, а также к обзорам Kosloff [42, 47,43, 44].

Предпосылки создания нового метода

Рассмотрим в целом описанные в предыдущем разделе приближенные методы решения задачи (1)-(2) и представим их недостатки как предпосылки разработки нового метода и сформулируем те требования, которым он должен удовлетворять.

• Квазиклассические методы. Основой для построения этих методов является малость постоянной Планка h 1. При этом типичная качественная картина процесса будет следующей: система представляет собой узкий в проекции на фазовое пространство (х,р) волновой пакет, центр которого перемещается с большой точностью согласно классическим законам движения. Эта ситуация схематически изображена на Рис. 1 слева. Изложенные выше квазиклассические методы имеют асимптотический характер и их погрешность имеет вид 0(Лт)до, т > 0. В реальных задачах значение постоянной Планка хоть и мало, но фиксировано. Поэтому практически невозможно построить универсальный метод, с помощью которого можно было бы получить решение широкого класса задач с любой наперед заданной степенью точности. Кроме того, например, в методе гауссовых волновых пакетов [34, 31] предполагается, что решение всеща имеет гауссов профиль, что, строго говоря, верно лишь для квадратических потенциалов. По мере приближения к квантовому случаю h ~ 1 погрешность связанная с этим ограничением начинает проявляться отчетливее. Далее, по построению, центр пакета движется вдоль классических траекторий, что делает невозможным описание большинства неклассических Р к<: t X X

Рис. 1: Качественная картина поведения рассматриваемой системы в проекции на фазовую плоскость. Слева — ситуация, характерная для квазиклассического случая h 1; справа — квантовый режим h ~ 1. явлений, например, туннелирования. Несмотря на постоянно проводимую работу по модификации метода [34], квазиклассические методы остаются, как правило, неэффективными при рассмотрении квантового случая.

• Квантовые методы. С другой стороны, использование изложенных выше квантовых подходов наталкивается на трудности при реализации в квазиклассическом случае. Например, погрешность аппроксимации пространственной части конечно-разностными методами определяется градиентом решения (скоростью роста), который становится большим при малых h. Это приводит к необходимости использования мелких сеток. Применение основных сеточных методов к квазиклассическому случаю было исследовано в работе [21]. Здесь показывается, что даже устойчивые сеточные методы требуют жестких ограничений на шаги сетки /гиг для того чтобы адекватно воспроизвести возникающие особенности решения. При этом для спектрального метода эти ограничения менее жесткие, чем для разностной схемы, но тем не менее имеют вид: h = 0(h) ит = o(h), h —* 0. Кроме того, обычно нет возможности использовать тот факт, что в каждый отдельный момент времени решение отлично от нуля лишь на малом подмножестве этой подробной сетки. Это приводит на практике к большому объему лишних вычислений и затратам памяти. Таким образом, сеточные методы эффективны в случае h ~ 1, когда для того, чтобы учесть все квантовые эффекты, необходимо непосредственное решение задачи (1)-(2), а само решение представляет собой медленно меняющуюся в области, где рассматривается процесс, функцию (см. на Рис. 1 справа).

Но часто возникает промежуточная ситуация, изображенная на Рис. 2. Здесь мы

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

Целью настоящей работы является разработка нового приближенного метода решения задачи (1)-(2), который справился бы с описанными трудностями и удовлетворял бы следующим требованиям Р X

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

• являлся бы точным методом, т.е. параметры метода можно всегда выбрать таким образом, чтобы получить решение задачи с произвольной заданной точностью;

• позволял бы эффективно рассматривать динамику волновых пакетов (в случае, например, инфинитного движения).

Краткое содержание работы. Описание метода

Итак, поставлена проблема построения приближенного метода решения задачи (1)-(2). С этой целью в настоящей работе используется модификация уже упомянутого выше метода гауссовых волновых пакетов, получившего широкое распространение при рассмотрении квазиклассического случая [32, 31]. В частности нами используется результат работы [4], где эта методика была применена для аппроксимации функции Грина оператора Шрёдингера и была получена соответствующая оценка погрешности. Сразу же отметим, что всюду в настоящей работе задача (1)-(2) рассматривается в одномерном по пространству случае d, = 1. Рассмотрение одномерного случая не является существенным ограничением, так как все конструкции допускают естественное обобщение на многомерный случай. В то же время появляется возможность наглядно и с достаточной строгостью описать новый метод и иллюстрировать его применение.

Изложение материала в диссертации разбито на три главы. В первой главе приведены необходимые сведения и результаты из теории гауссовых волновых пакетов. Во второй главе строится численный метод решения задачи (1)-(2) и дается его теоретическое обоснование. Работу завершает третья глава, где обсуждаются результаты расчетов новым методом некоторых модельных задач.

В методе гауссовых волновых пакетов уравнению (1) ставится в соответствие система обыкновенных дифференциальных уравнений = Р(0, Q( о) = 9о i f—P(0) = pc,

HA <B(i), Л(0) = a, (8) i,

§ = 5 - 5(0) = o.

Здесь начальные данные qQ и po произвольны, a a > 0 - параметр, значение которого будет дано ниже. На решении этой системы определяется функция

•»> = «ш: ех» t^1 + '£<« -« + '!} ■ (9)

Метод гауссовых волновых пакетов основывается на том, что функция W(x,t\qo,po) удовлетворяет уравнению (1) в случае квадратичного потенциала точно на любом временном интервале (см., например, [30]). В случае же произвольного потенциала — W может быть рассмотрена как приближенное решением уравнения (1) на достаточно малом отрезке времени, удовлетворяющее начальному условию специального вида

М*\Яо,Ро) = 2(7^)3/^/2 ехр + 1Т{Х - 9о)} (10) гауссову волновому пакет,у8. Таким образом, функция W есть гауссов волновой пакет с начальной дисперсией о и центром, перемещающимся по классической траектории, исходящей из точки (<?о,Ро)• Функции (10) образуют базис в пространстве состояний, что выражается формулами для прямого

F : - Fh { <р | a, q,р } = J+°° ехр ~ " 9)} Ф) dx (11) и обратного преобразования Фурье-Гаусса (FBI-преобразование в [27]; преобразование Габора в [17]) 2(^)3/^/2 JR2 еХР {—inf + »%(* " ?)} Р К g, Р } dqdp, (12)

8Строго говоря, функция в правой части формулы (10) отличается от гауссова волнового пакета нормировочным множителем. Но это в данном случае не изменяет сути.

Преобразование (11) определяется на функциях ip G Co°(R), р, q € R,cr > 0. На функциях из L2(R, dx) оно распространяется на основе равенства Парсеваля. Далее, на Cq°(R) определим оператор fC{t, а) : 1р0(х) K(t, а)ф0(х) = W{x, t\q0,p0)Fh { ф0 | a, qQ,p0 } dq0 dp0. (13)

Пусть теперь U(t,h) — оператор эволюции задачи (1)-(2). Под потенциалом V(х)

Основным результатом работы [4] является оценка в операторной норме погрешности аппроксимации оператора эволюции задачи (1)-(2) U(t, h) степенью оператора К

Теорема 1 ([4]). Если потенциал удовлетворяет условию (14), то существует такое е > 0, зависящее только от потенциала V(x), и для любого фиксированного Т < оо такая не зависящая от h,0 < h < const, константа а(Т) < оо, что при всех целых п и 0 ^ t ^ тт(ещТ) выполнена оценка где || • || - операторная норма в L2(R, dx).

Теорема 1 утверждает аппроксимацию оператора эволюции задачи степенью оператора 1С и гарантирует, что при этом можно добиться произвольной точности, выбирая достаточно большое п. Последнее является основанием для поиска приближенного решения задачи (1)-(2) по временным слоям, последовательно применяя оператор /С. Но на практике расчеты непосредственно по этой схеме потребовали бы постоянного переключения между пространственным представлением и представлением Фурье-Гаусса с помощью преобразования (11). Чтобы этого избежать, предлагается рассмотреть FBI—образ функции W, который можно вычислить явно, используя интеграл Пуассона (см. параграф 1.3) будем понимать функцию из Z/(R, dx), образ Фурье V(k) которой удовлетворяет неравенству

IC{t/n,t/n)n-U(t,h)|| ^ a{T)til2n-z'2,

15) G(t,q,p,q0,p0,a). (16)

Тогда, если задаться FBI-преобразованием Ф0(<7>р) начального условия чро(х), то, согласно теореме 1, мы можем всегда выбрать: достаточно маленький9 шаг по времени т, FBI-параметр о и соответственно достаточно большую степень п, чтобы получить с любой заданной точностью FBI-образ Ф решения задачи в произвольный фиксированный момент времени Т = пт

П pai ---ч

T,Sn)= J ••• JG(t,Zn,Sn-1).G(t,Zl,Zo)MSo)dSo.dtn-l, (17) r2" где

Zj = {Qj,Pj), \Sj\2 = Qj+Pp = dqj dpj, j = 0,1,. ,n - 1.

Таким образом, функция G является аппроксимацией функции Грина задачи в представлении Фурье-Гаусса. При этом, для того чтобы вычислить функцию G в заданных точках, необходимо решить на отрезке [0, t] систему обыкновенных дифференциальных уравнений (8) для параметров пакета W и подставить значения полученного решения (взятые в точке t) в формулу (16) для G. Процедура построения функции G и необходимые сведения из теории гауссовых волновых пакетов подробно изложены в первой главе диссертации.

Вторая глава диссертации посвящена построению численного метода решения задачи (1)-(2) и его теоретическому обоснованию. Содержание главы полностью отражено автором в публикациях [1, 2, 3].

Формула (17) выступает в работе в качестве отправной точки при построении численной процедуры. Она дает основание искать приближенное решение задачи на фиксированном отрезке времени [0,Т] последовательно по временным слоям {tk = кт, к = 0, l,.,n; tn = пт = Т} / Ы1 ей2, к = 0,гг — 1. (18)

J r2

Бесконечная область интегрирования в (17),(18) является главным препятствием при построении численного алгоритма. Но эта проблема может быть решена при условии быстрого убывания на бесконечности начального условия Фо и в силу свойства функции G, которое устанавливается в работе (см. [2]). А именно: на основании проведенного в работе исследования функции G, в качестве основного результата,

9т не должен превосходить некоторого е, которое определяется только потенциалом. доказана следующая теорема о характере приближенного решения, получаемого по формуле (17)

Теорема 2. Пусть потенциал V(x) удовлетворяет условию (14), а функция Фо(£о)> £о £ R2 такова, что

Фо(£о)| ^ С0е-°оЫ2, для всех & е R2, где Со > 0, Qo > 0 - некоторые постоянные. Тогда

1. С любой степенью точности можно заменить интегрирование в (17) по бесконечной области интегрированием по некоторой конечной области. Причем эта область может быть выбрана одной и той же для всех £п 6 R2. (равномерная по параметру £n € R2 сходимость интеграла (17))

2. Приближенное решение Ф(£п, £п). полученное по формуле (17), также является быстро убывающей на бесконечности функцией: существуют положительные постоянные С„, ап такие, что

Ф(*п,£пЖ Cne"^2, V£neR2. (19)

Данная теорема позволяет непосредственно перейти к построению численного алгоритма решения задачи. В работе сначала был предложен численный алгоритм расчетов на фиксированной сетке [1], а затем его модификация, использующая подвижную сетку [3].

Не ограничивая общности, можно считать, что в качестве начального условия ф0 выступает гауссов волновой пакет вида (10). Его образ Фурье-Гаусса Ф0(<7о,Ро) является быстро убывающей на бесконечности функцией. Поэтому, согласно теореме 2, решение, получаемое по формуле (18), будет также быстро убывающей на бесконечности функцией, которую с любой заданной точностью можно считать отличной от нуля лишь в некоторой конечной области - "носителе". При этом интегрирование в (18) можно сколь угодно точно заменить интегрированием по достаточно большой, но конечной и фиксированной области плоскости (q,p), которая бы содержала носители приближенного решения на каждом временном слое tk из [0,Т]. В качестве такой области в работе был выбран прямоугольник D на плоскости (q, р) со сторонами, параллельными координатным осям. После введения в D равномерной сетки Г2 с шагами hq и hp, и приближенного интегрирования в (18) по квадратурной формуле, была получена следующая численная процедура решения задачи (1)-(2) hghp^Gfatk+i,Sk,Sk+i е П, к = О^Т, (20) аеп где {VI,к(£к), € ft} — вектор (матрица) узловых значений приближенного решения задачи на сетке П на к-ом слое по времени. Параметрами метода, определяющими точность вычислений, являются: шаг по времени г, FBI-параметр а, размер и расположение прямоугольника D, шаги сетки Q. Эти параметры должны фиксироваться до начала вычислений и выбираться эмпирически. Таким образом, получается следующий алгоритм расчетов на фиксированной сетке: шаг 1) Задаемся образом Фурье-Гаусса начального условия Ф0; шаг 2) Фиксируем параметры метода. Для выбранных значений параметров определяем по формуле (16) функцию G(r,в узлах сетки; шаг 3) Далее приближенное решение находится последовательно по временным слоям по формуле (20).

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

Gh(r) = € П}, которая соответствует элементарному шагу по времени т. В работе описывается процедура возведения в квадрат матрицы Gu- После тг—кратного применение этой операции получается матрица пропагатора с эффективным шагом 2пт, которая может быть сохранена и использована в дальнейшем. Это позволяет существенно повысить общую эффективность вычислительного процесса в случаях, когда величина начального шага г много меньше отрезков времени, на которых рассматривается задача.

Описанная только что численная процедура решения задачи (1)-(2) на фиксированной сетке оказывается неэффективной, когда решением является волновой пакет, перемещающийся в фазовой плоскости (q,p). В этом случае требуется выбор изначально большого прямоугольника D, который смог бы "покрыть" все возможные положения пакета на рассматриваемом отрезке времени [О,Т]. При этом в каждый отдельный момент времени использование такой большой области неоправданно, т.к. решение отлично от нуля лишь на малом подмножестве D. Поэтому применение процедуры (20) потребует большого количества лишних вычислений и памяти. В связи с этим возникает необходимость в модификации алгоритма за счет введения подвюкной сетки. Перейдем к изложению алгоритма.

Вместо прямоугольника D рассмотрим на каждом слое по времени к наименьший прямоугольник D(k) со сторонами, параллельными координатным осям, и содержащий множество всех точек плоскости, где значения решения по абсолютной величине превышают некоторый малый порог S. В прямоугольнике D(k) вводится сетка Q(k), состоящая лишь из "существенных" на текущем слое узлов сетки Q. Вместо (20) используется следующая схема вычисления по слоям

4+1(&+i) = mp е <7(т,й+1,&)ф*(&), к = 0^1. (21) en(fc)

С течением времени прямоугольник D(k) может перемещаться в плоскости (q,p) и изменять свой размер. Соответствующие преобразования сетки Q(k), связанные с добавлением и исключением узлов при переходе на следующий слой по времени, предлагается выполнять согласно следующему алгоритму шаг 1) Пусть известно решение vPfc(£fc) на сетке Q(k). Вычисляем по формуле (21) решение на следующем временном слое во всех узлах Q(k); шаг 2) Анализируем найденные значения и удаляем лишние узлы из сетки. Для этого значения в узлах сетки Q(k) просматриваются вглубь по строкам или столбцам, в зависимости от выбранной стороны прямоугольника D{k). Из сетки удаляются строки (столбцы) узлов, в которых все значения решения меньше порогового по модулю: |Ф*+1(&)| < 6, до тех пор пока не встретится строка (столбец), где есть хотя бы одно значение, превышающее 5. После этого указанная процедура повторяется, начиная со следующей стороны. Таким образом рассматриваются все стороны прямоугольника в определенной последовательности (например, по часовой стрелке, начиная с левого края); шаг 3) Рассматривается возможность добавления в сетку новых узлов. Сетка может быть расширена на сторонах, где прямоугольник не был обрезан. Если такие стороны есть, то они рассматриваются последовательно. При этом значения решения vI'fc+1 в потенциально новых узлах, не принадлежащих сетке П(к), рассчитываются, как и прежде, по формуле (21). Сетка расширяется за счет строк (столбцов) узлов, которые добавляются с соответствующих сторон, до тех пор пока в них есть хоть одно значение решения, превышающее по модулю порог S.

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

Третья глава диссертации содержит результаты расчетов модельных задач, на примере которых доказывается практическая применимость предложенного метода. Были рассмотрены следующие постановки, которые определяются выбором конкретного потенциала V(x) в уравнении (1)

- гармонический осциллятор;

- равноускоренное движение частицы;

- движение частицы в двухъямном потенциале;

- движение в потенциале Морзе;

- взаимодействие с потенциальным барьером, туннелирование.

Перечисленные задачи являются характерными примерами финитного и инфинит-ного движения системы.

В качестве результата, помимо самого решения, находились некоторые наблюдаемые, такие как среднее значение координаты частицы, норма волновой функции10. Также отслеживалось перемещение "носителя" приближенного решения (сетки) по фазовой плоскости. Получаемые значения сравнивались (где это возможно) с известными аналитическими решениями. Расчеты проводились для различных значений параметров метода с целью продемонстрировать их влияние на точность. Содержание третьей главы полностью отражено в работах [1] и [3].

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

Место метода и научная новизна

Обозначим место предложенного метода среди уже существующих подходов. Итак, в работе предложен новый сеточный метод решения задачи (1)-(2), который может быть рассмотрен как

• модификация метода гауссовых волновых пакетов и нестационарного метода Хартри-Фока, позволяющая за счет выбора параметров метода получать решение задачи с любой наперед заданной точностью;

• метод подвижной сетки для решения нестационарного уравнения Шрёдингера.

Основные результаты

Сформулируем еще раз основные результаты работы, которые выносятся на защиту

1. Предложен новый сеточный метод решения задачи Коши для нестационарного уравнения Шрёдингера, использующий подвижную сетку.

2. Сформулирована и доказана теорема, являющаяся теоретическим обоснованием предложенного в п. 1 метода. Проведены численные расчеты, результаты которых доказывают практическую применимость предложенного метода.

В заключение автор выражает глубокую признательность своему научному руководителю Алексею Алексеевичу Арсеньеву за постановку задачи и постоянное внимание к работе, а также коллективу кафедры вычислительных методов за поддержку и ценные рекомендации в процессе работы.

Список литературы диссертационного исследования кандидат физико-математических наук Айдагулов, Галлям Раильевич, 2004 год

1. Айдагулов Г.Р. Применение преобразования Фурье-Гаусса к решению задачи Ко-ши для уравнения Шрёдингера // ЖВМиМФ. 2002. Т.42. №12. С.1810-1815.

2. Айдагулов Г.Р. Применение преобразования Фурье-Гаусса к решению задачи Ко-ши для уравнения Шрёдингера: теоретический анализ численного алгоритма // Вычислительные методы и программирование. 2003. Т.4. №2. С. 16-20.

3. Айдагулов Г.Р. Метод подвижной сетки для решения нестационарного уравнения Шрёдингера // Вычислительные методы и программирование. 2004 (январь). Т.5. Раздел 1. С.18-30. (http://num-meth.srcc.msu.su/)

4. Арсенъев А.А. Оценка функции Грина оператора Шрёдингера // Теор. и матем. физ. 1998. Т.115. №1. С.85-92.

5. Арсенъев А.А. Аппроксимация функции Грина уравнения Шредингера // Дифф.уравнения 1999. Т.35. №3. С.319-324.

6. Бахвалов Н.С. Численные методы. М.: Наука, 1975.

7. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987.

8. Волкова Е.А., Попов A.M., Рахимов А.Т. Квантовая механика на персональном компьютере. М.:УРСС, 1995.

9. Давыдов А.С. Квантовая механика. М.: Наука, 1973.

10. Жидков Е.П., Лобанов Ю.Ю. Метод приближенного континуального интегрирования и некоторые его приложения. // Матем. моделирование. 1999. Т.Н. №5. С.37-83.

11. Калиткин Н.Н. Численные методы. М.: Наука, 1978.

12. Ландау Л Д., Лифшиц Е.М. Квантовая механика. Нерелятивистская теория. М.: Наука, 1989.

13. Марчук Г.И. Методы вычислительной математики. М.:Наука, 1980.

14. Рихтмайер Р. Принципы современной математической физики. T.l. М.: Мир, 1982.

15. Самарский А. А. Теория разностных схем. М.:Наука, 1989.

16. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: МГУ, 1999.

17. Чуй К. Введение в вэйвлеты. М.: Мир, 2001.

18. Askar A.,Cakmak A.S. Explicit integration method for the time-dependent Schrodinger equation for collision problems // J.Chem.Phys. 1978. Vol.68. N6. P.2794-2798.

19. Bandrauk A.D., Shen H. Improved exponential split operator method for solving the time-dependent Schrodinger equation // Chem.Phys.Lett. 1991. Vol.176. N5. P.428-432.

20. Bao W., Jin S., Markowich PA. On time-splitting spectral approximations for the Schrodinger equation in the semiclassical regime // J.Comp.Phys. 2002. Vol. 175. P.487-524.

21. Bisseling R., KosloJfR. The fast Hankel transform as a tool in the solution of the time dependent Schrodinger equation // J.Comp.Phys. 1985. Vol.59. P.136-151.

22. Bracewell R.N. The Fourier transform and its application. McGraw-Hill, New-York, 1978.

23. Burghardt В., Stolze J. Numerical evaluation of coherent-state path integrals in quantum dynamics // J.Phys.A: Math.Gen. 1999. Vol.32. P.2075-2089.

24. Dong S., Meng H. Chebyshev spectral method and Chebyshev noise processing procedure for vorticity calculation in PIV post-processing // Experimental Thermal and Fluid Science. 2001. Vol.24. P.47-59.

25. Feit M.D., Fleck J.A Jr., SteigerA. Solution of the Schrodinger equation by a spectral method // J.Comp.Phys. 1982. Vol.47. P.412-433.

26. Folland G.B. Harmonic analysis in phase space. Princeton Univ., Princeton, 1989.

27. Grossmann F. Inelastic resonant tunneling with wavepackets // Chem.Phys. 2001. Vol.268. P.347-353.

28. Hagedorn4 G.A. Semiclassical quantum mechanics I: The h —* 0 limit for coherent states. // Commun. Math. Phys. 1980. Vol.71. P.77-93.

29. Hagedorn G.A. Raising and lowering operators for semiclassical wave packets // Ann. Phys. 1998. Vol.269. P.77-104.

30. Hagedorn G.A., Joye A. Exponentially accurate semiclassical dynamics: Propagation, localization, Ehrenfest times, scattering, and more general states // Ann. Henri Poincare 2000. Vol.1. P.837-883.

31. Heller E.J. Generalized theory of semiclassical amplitudes // J.Chem.Phys. 1977. Vol.66. N12. P.5777-5785.

32. Herman M.F. Dynamics by semiclassical methods // Annu. Rev. Phys. Chem. 1994. Vol.45. P.83-111.

33. Huber D., Heller E.J. Generalized Gaussian wave packet dynamics // J.Chem.Phys. 1987. Vol.87. N9. P.5302-5311.

34. Iracane D., Bernis L. Finite-basis-set technique for solving time-dependent Schrodinger equation // Phys.Rev.A. 1990. Vol.42. N11. P.6831-6844.

35. Manthe U., Koppel H. New method for calculating wave packet dynamics: strongly coupled surfaces and the adiabatic basis //J.Chem.Phys. 1990. Vol.93. N1. P.345-356.

36. Koonin S.E., Davies K.T.R, Maruhn-Rezwatii V., Feldmeier II, Krieger S.J. Time-dependent Hartree-Fock calculations for 160+160 and 40Ca+40Ca reactions // Phys.Rev.C. 1977. Vol.15. N4. P. 1359-1374.

37. KosloJfD., Kosloff5 R. A Fourier method solution for the time dependent Schrodinger equation as a tool in molecular dynamics // J.Comp.Phys. 1983. Vol.52. P.35-53.

38. Kosloff R., Kosloff D. A Fourier method solution for the time dependent Schrodinger equation: A study of the reaction tf+ + H2, D+ + HD, D+ + H2 II J.Chem.Phys. 1983. Vol.79. N4. P. 1823-1833.

39. Kosloff R., Tal-Ezer H. A direct relaxation method for calculating eigenfunctions and eigenvalues of the Schrodinger equation on a grid // Chem.Phys.Lett. 1986. Vol. 127. N3. P.223-230.

40. Kosloff R. Time-dependent quantum-mechanical methods for molecular dynamics // J.Phys.Chem. 1988. Vol.92. N8. P.2087-2100.

41. KosloffR. Propagation methods for quantum molecular dynamics // Annu. Rev. Phys. Chem. 1994. Vol.45. P.145-178.

42. Kosloff R. Dynamics of Molecules and Chemical Reactions, In R. E. Wyatt and J. Z. Zhang, editors, Dynamics of Molecules and Chemical Reactions, P. 185-230, Marcel Dekker, New York, 1996.

43. Lee S.-Y., Heller E.J. Exact time-dependent wave packet propagation: Application to the photodissociation of methyl iodide//J.Chem.Phys. 1982. Vol.76. N6. P.3035-3044.

44. Leforestier C. Competition between dissociation and exchange processes in collinear A + ВС collision. I. Exact quantum results. // Chem.Phys. 1984. Vol.87. P.241-261.

45. McCullough E.A., WyattR.E. Dynamics of the collinear #+#2 reaction I. Probability density and flux // J.Chem.Phys. 1971. Vol.54. N8. P.3578-3591.

46. Mowrey R.C., Sun X, Kouri D.J. A numerically exact full wave packet approach to molecule-surface scattering//J.Chem.Phys. 1989. Vol.91. N10. P.6519-6524.

47. Pindzola M.S., Robicheaux E, Gavras P. Double multiphoton ionization of a model atom // Phys.Rev.A. 1997. Vol.55. N2. P. 1307-1313.

48. Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical recipes in Fortran 77: The art of scientific computing. Vol.1. Cambridge University Press, 1992.

49. Robinett R.W. Quantum wave packet revivals // Phys.Pep. (to appear) (arXiv:quant-ph/0401031)