Содержание к диссертации
Введение
1 Метод гауссовых волновых пакетов 26
1.1 Когерентные состояния 26
1.2 Метод гауссовых волновых пакетов 28
1.2.1 Квадратичный потенциал 28
1.2.2 Потенциал: общий случай . 31
1.2.3 Модификация метода 34
1.3 Вычисление наблюдаемых 39
2 Построение численной процедуры 47
2.1 Предварительные замечания 47
2.2 Доказательство теоремы 50
2.3 Численный алгоритм 55
2.3.1 Фиксированная сетка 55
2.3.2 Подвижная сетка 56
2.4 Связь с методом частиц 59
3 Результаты расчетов 61
3.1 Аналитическое решение 61
3.2 Гармонический осциллятор 63
3.3 Равноускоренное движение частицы 66
3.4 Двухъямный потенциал 69
3.5 Потенциал Морзе 71
3.6 Рассеяние на потенциальном барьере 74
Библиография
- Когерентные состояния
- Численный алгоритм
- Равноускоренное движение частицы
- Рассеяние на потенциальном барьере
Введение к работе
Предметом настоящей работы является разработка приближенного метода решения задачи Коши для нестационарного уравнения Шредингера
ih^(x1t) = ~Atl>(x,t) + V(x)iP(x,t)t ieRrf, teRl+1 (1)
^(х,0 + 0) = фа[х), x Є Rd, (2)
где i}(-,t) Є L2(Rd, dx), t Є R.+. В работе предлагается новый численный метод решения задачи. Также проводится ряд численных экспериментов с использованием нового метода, на основе результатов которых демонстрируется его практическая применимость.
Уравнение Шредингера (I) — одно из основных уравнений математической физики, описывающее эволюцию квантовой системы во времени. Из аксиом квантовой механики [9, 12] следует, что волновая функция системы, находящейся во внешнем поле с потенциалом V = V(x), есть решение задачи Коши (1)-(2), где -ф0 — волновая функция в начальный момент времени. Задача (1)-(2) возникает, например, при теоретическом описании явлений в наноэлектронике, полупроводниках, квантовой химии и физике твердого тела. В общем случае уравнение (1) может быть решено только численно. В виду интенсивного развития математического моделирования в этих областях возникает проблема разработки эффективных численных методов решения задачи (1)-(2)-
Обзор существующих методов
Для того, чтобы продемонстрировать необходимость разработки нового метода, сформулировать четкие требования, которым он должен удовлетворять, рассмотрим паи-
Введение
более употребительные в настоящее время подходы к решению задачи (1)-(2). В контексте настоящего исследования это удобно сделать следующим образом.
Постоянная Планка h в уравнении (1) играет фундаментальную роль во всех квантовых явлениях. Ее относительная величина (по сравнению с другими величинами той же размерности) определяет "степень квантовости" той или иной физической системы. В настоящей работе под h мы будем понимать параметр, входящий в уравнение (1), которое было обезразмерено при выборе системы единиц, наиболее всего отвечающей масштабу рассматриваемого физического явления. При этом формально квантовая механика содержит в себе классическую в качестве предельного случая при h —* 0. Например, при описании движения электрона вокруг ядра, где характерное расстояние имеет порядок Ю-10 м, в соответствующей (атомной) системе единиц величина h равна 1. В задачах физики твердого тела при длинах порядка 1СГ7 м, h имеет порядок Ю-2. С дальнейшим уменьшением h квантовые эффекты проявляются все слабее, и система начинает подчиняться классическим законам. Случай малого значения параметра h (/і
Оказывается, что и численные методы решения задачи (1)-(2) могут быть разделены в зависимости от величины h на квазиклассические (h <; 1) и квантовые (h -— 1). При этом метод одного класса часто оказывается неэффективным в применении к задачам другого. Поясним это, кратко рассмотрев существующие методы.
Квазиклассические методы (случай h <; 1)
В основе всех квазиклассических методов лежит тот факт, что при /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] 0(/іт)д_оі т > 0. Это не позволяет получить решение с произвольной точностью в практических расчетах, в которых значение h хоть и мало, но все же фиксировано. Поэтому границы применимости получаемых приближений, вообще говоря, неясны, и требуют дополнительного исследования в каждом конкретном случае. Кроме того, по мере приближения к квантовому режиму (h ~ 1) погрешности проявляются сильнее, что делает метод неприменимым в этом случае.
Важным преимуществом квазиклассических методов является меньший объем вычислительных затрат, что делает их менее чувствительными к размерности задачи по сравнению с квантовыми методами. А получаемые при этом результаты, выраженные в терминах классических траекторий, допускают простое интуитивное толкование.
Для получения более подробной информации по квазиклассическому приближению и методу гауссовых волновых пакетов можно порекомендовать работы [42, 31, 34], обзор [33], а также цитируемую там литературу.
Квантовый случай (ft ~ 1)
В квантовом случае (ft ~ 1) наиболее употребительными являются сеточные методы решения (см. обзор Kosloff [43] и ссылки там), когда пространственная и временная
Введение
части уравнения (1) аппроксимируются на сетке. Это представляется вполне естественным, так как величина постоянной Планка в данном случае определяет такой масштаб рассматриваемых явлений, что, вообще говоря, необходимо "использовать" всю информацию, которую дает задача (1)-(2), для адекватного описания процесса. Например, непригодность в квантовом случае квазиклассичсских методов, может быть обусловлена именно тем, что с самого начала не принимаются во внимание некоторые свойства решения, которыми можно пренебречь на макроуровне (h <2С 1), но которые начинают играть существенную роль при переходе к квантовому режиму (h ~ 1).
Задача (1)-(2) это прежде всего нестационарная задача, описывающая развитие некоторого процесса во времени. Поэтому построение сеточных методов решения задачи принято разбивать на два основных и независящих друг от друга этапа. Во-первых, это выбор представления, в котором ищется приближенное решение (например, пространственное или импульсное) и его последующая дискретизация за счет введения сетки. По сути, на этом этапе решается вопрос о том, как аппроксимировать пространственную часть уравнения. Далее, решение ищется как сеточная функция последовательно по временным слоям. Выбор способа перехода на следующий слой по времени, или аппроксимация временной части, н является второй составляющей метола.
Для удобства будем считать, что уже сделана аппроксимация пространственной части на первом этапе, и рассмотрим сначала временную составляющую. Пусть оператор Ид - некоторый сеточный аналог2 гамильтониана3 Ті в правой части уравнения (1), Тогда переход на следующий временной слой формально определяется с помощью аппроксимации оператора эволюции
і«,(0 = ехр{-ф<Л. (3)
Здесь можно отметить следующие подходы
^Индекс d здесь означает именно дискретизацию, и, естественно, никак не связан с размерностью
пространства в записи задачи (1)-(2).
эОператор правой части в уравнении (I)
называется оператором Гамильтона и выделяется особо, В квантовой механике этот оператор соответствует величине полной энергии системы.
Введение
Двуслойная схема Кранка-Николсона (Crank-Nicholson, CN). Эта схема может быть интерпретирована как Паде-аппроксимация экспоненты (3)
где Z(f — единичный оператор, действующий в пространстве сеточных функций. Если оператор W.& - самосопряжен, то схема Кранка-Николсона - эрмитова, абсолютно устойчива, сохраняет L2—норму решения. Основным недостатком схемы является ее неявность.
Центральная разность (second order differencing, SOD). SOD-схема получается, если производную по времени в левой части (1) аппроксимировать центральной разностью. Получается явная трехслойная схема того же порядка по времени, что и CN. Но она является условно устойчивой и лишь приближенно унитарной.
расщепления (split operator technique, SP). При построении схем расщепления в применении к уравнению ШрёдНигера4, как правило, рассматривают приближенную факторизацию оператора эволюции (3) с помощью различных вариаций формулы Троттера. Это хоть и не дает законченного алгоритма перехода на следующий временной слой, но позволяет свести исходную задачу к цепочке более простых. Впервые этот подход был предложен в работе Felt [26], где выделялись кинетическая и потенциальная составляющие гамильтониана 1С = -Цд и V = V(x): H = K + V. Тогда
exp{-it?{/h} = ехр{-0.5ШС//І} exp{-i(V/ft} exp{-0.5#/C/ft} + 0(t3), (4)
что позволяет свести решение исходной задачи к последовательному решению двух стандартных подзадач: свободное движение и локальное взаимодействие волны частицы с потенциалом. Забегая вперед, скажем, что каждая из них может быть решена точно в пространственном или импульсном представлении соответственно. Необходимые переключения между представлениями осуществляются с помощью дискретного преобразования Фурье, использование же алгоритма быстрого преобразования Фурье (БПФ) позволяет сделать эти переходы чрезвычайно эффективными и точными. Это делает метод расщепления
4Обшее изложение метода расщепления при построении численных алгоритмов можно найти в [13].
Введение
Feit одним из наиболее часто используемых и в настоящее время, см., например, [28]. Также с его помощью получено большинство результатов численного моделирования в книге [8].
В качестве модификации метода можно отметить работу [20], где вместо (4) используется более точная факторизация, состоящая из семи множителей и имеющая погрешность порядка C(f4). Кроме того, выбор того или иного способа расщепления гамильтониана определяется спецификой задачи и предполагаемого дальнейшего способа решения получаемых подзадач. В [24] требовалось найти пропагатор задачи в представлении когерентных состояний гармонического осциллятора. Поэтому в исходном гамильтониане выделялась гармоническая и негармоническая части. В этом случае эволюция когерентного состояния под действием гармонического потенциала столь же тривиальна, как и свободное движение плоской волны в стандартном методе расщепления Feit (см, главу 1 и метод гауссовых волновых пакетов).
Глобальный пропагатор (global propagator). Все рассмотренные до сих пор приближения оператора эволюции (3) строились за счет выбора достаточно малого шага по времени и могут быть получены из ряда Тейлора
Ud{t) ^Td~ijKd - ~-Н% +.... (5)
Можно же рассмотреть задачу о построении так называемого глобального про-пагатора, который бы позволил находить приближенное решение задачи в любой, не обязательно малый, момент времени t "за один шаг". Отмечается, что использование с этой целью разложения Тейлора (5) оказывается неустойчивым. Поэтому в [53] был предложен метод, в котором для аппроксимации оператора эволюции используется отрезок ряда по системе полиномов Чебышева (СН) вида
Щі) ^Y.^Pni-iyTi^j . (6)
Коэффициенты ап выражаются через значения функций Бесселя, а необходимые при реализации алгоритма выражения Рп (—itTid/h) ^от п = 0,1,..., N могут быть найдены последовательно в ходе вычислений по известным рекуррентным формулам для многочленов Чебышева. В виду быстрой сходимости ряда, пред-
Введение
ложенный метод имеет высокую точность, которая определяется выбором не шага по времени, а количеством членов удержанных в разложении (6).
Одним из принципиальных недостатков глобального пропагатора является то, что он не дает возможности находить решение задачи в промежуточные моменты времени. Но такой подход оказывается чрезвычайно удобным когда требуется определить решение "на больших временах". Например, при нахождении основного состояния системы.
Рассмотрим теперь методы аппроксимации пространственной части. Комбинируя их с представленными выше методами перехода на следующий слой по времени, можно получить конкретные численные алгоритмы решения задачи (1)-{2).
Конечные разности (finite differencing, FD). Считается, что численное решение нестационарного уравнения Шрёдингера началось именно с метода конечных разностей. Одной из первых в этой области считается работа [48]. В ней для аппроксимации пространственной производной в правой части уравнения (1) используется вторая разностная производная, а шаг по времени делается с помощью рассмотренной выше схемы Кранка-Николсона. Таким образом получается двуслойная разностная схема, которая является абсолютно устойчивой, сохраняет норму решения и имеет второй порядок аппроксимации по пространству и времени. В отечественной литературе обсуждение разностных схем для уравнения Шрёдингера можно найти в книге [15]. Здесь рассматривается обший случай двуслойной схемы с весами и показывается, что устойчивые схемы есть только среди неявных, причем норму решения сохраняет только схема Кранка-Николсона. По этой причине эта схема в течение некоторого времени была наиболее часто используемым методом. Необходимость использовать неявные схемы усложняет решение задачи в многомерном случае. Тем не менее такие расчеты проводились. В работе [46] двумерный лапласиан аппроксимируется на тринадцатиточечном шаблоне, а возникающая при этом система уравнений решается итерационным методом. Один из способов сократить объем вычислений - использовать экономичные методы решения [15]. Например, в [38, 50] для решения задачи в двумерном случае используется схема переменных направлений (Peaceman-Rachford), которая может быть рассмотрена как приближенная факторизация (расщепление) схемы Кранка-Николсона.
Введение
Трехслойная явная устойчивая схема предложена в работе [19]. Для аппроксимации по пространству здесь используется, как и выше, вторая разность, а по времени - центральная, т.е. SOD-пропагатор. Эта схема оказывается условно устойчивой, что накладывает жесткие ограничения на возможные шаги сетки. Кроме того, она является лишь приближенно унитарной (общее свойство SOD). Но можно показать, что отклонение нормы решения от 1, на самом деле, на несколько порядков меньше погрешности самой схемы и этот недостаток незначителен, в сравнении с возможностью использования в многомерном случае [19]. Заметим, кстати, что в применении к уравнению теплопроводности такая схема оказывается неустойчивой.
Псевдоспектральный подход (pseudospectral method). Здесь мы рассмотрим подход5 к аппроксимации пространственной части уравнения (1), который часто используется при решении дифференциальных уравнений и состоит в том, чтобы с помощью некоторого интегрального преобразования или разложения по специальной системе функций перейти к новому представлению, в котором вычисление требуемых характеристик (например, производных) оказывается эффективнее, чем в исходном представлении.
В применении к уравнению Шрёдингсра наиболее употребительным является разложение по тригонометрической системе и для аппроксимации пространственной производной используется связь между образами Фурье самой функции и ее производной. Тогда оператор Лапласа Аф{х), х Є Rd, может быть вычислен с помощью последовательного нахождения (/-мерного преобразования Фурье, последующего умножения полученного результата на функцию —(fc? + kj + ... + k7d), к - {ki,...tkd) Є Rd, и вычисления обратного Фурье преобразования окончательного выражения. Дискретный аналог преобразования Фурье позволяет, применяя эту процедуру к вектору узловых значений функции ф> получить сеточную аппроксимацию6 лапласиана функции в правой
5Который в свою очередь можно рассматривать как частный случай метода коллокаций, см. [42,44], 6В [7, 6, ] 1] дискретное преобразование Фурье рассматривается в контексте тригонометрической
иитерпачяции. Последняя тем точнее, чем глаже функция. Это замечание позволяет рассматривать предлагаемый способ аппроксимации производной как еще один вариант использования интерполирования при построении формул чисяенкого дифференцировакия. Для подробностей Kosloff отсылает к [23]. Также см. ссылки по спектральному методу ниже.
Введение
части уравнения (I) и, как результат, приближенный аналог всего гамильтониана. Использование алгоритма БПФ позволяет существенно уменьшить объем вычислений.
К самому раннему использованию описанного метода для решения уравнения Шрёдингера можно отнести работу Feit [26], где с его помощью, по сути, решалась возникающая в ходе расщепления подзадача о свободном движении пакета. Обобщение на случай потенциала общего вида было предложено KosIoffB [39] и успешно опробовано при моделировании конкретной реакции в [40], Здесь для аппроксимации пространственной части использовалась описанная процедура, а шаг по времени делался по явной SOD-схеме. Следствием последнего является приближенная унитарность и условная устойчивость.
В [39] указывается основное ограничение на использование метода: решение задачи должно быть "представимо" на сетке. Это означает, что волновой спектр решения должен лежать внутри полосы, ширина которой определяется шагом пространственной сетки
Природа этого ограничения та же самая что и у частоты Найквиста в теории обработки сигналов. Хотя условие (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: Качественная картина поведения рассматриваемой системы в проекции на фазовую плоскость. Слева — ситуация, характерная для квазиклассического случая ft < 1; справа — квантовый режим ft ~ 1.
явлений, например, туннелирования. Несмотря на постоянно проводимую работу по модификации метода [34], квазиклассические методы остаются, как правило, неэффективными при рассмотрении квантового случал.
Квантовые методы. С другой стороны, использование изложенных выше квантовых подходов наталкивается на трудности при реализации в квазиклассическом случае. Например, погрешность аппроксимации пространственной части конечно-разностными методами определяется градиентом решения (скоростью роста), который становится большим при малых h. Это приводит к необходимости использования мелких сеток. Применение основных сеточных методов к квазиклассическому случаю было исследовано в работе [21]. Здесь показывается, что даже устойчивые сеточные методы требуют жестких ограничений на шаги сетки Лиг для того чтобы адекватно воспроизвести возникающие особенности решения. При этом для спектрального метода эти ограничения менее жесткие, чем для разностной схемы, но тем не менее имеют вид: h = 0(h) и т = o(h), h —* 0. Кроме того, обычно нет возможности использовать тот факт, что в каждый отдельный момент времени решение отлично от нуля лишь на малом подмножестве этой подробной сетки. Это приводит на практике к большому объему лишних вычислений и затратам памяти, Таким образом, сеточные
Введение
методы эффективны в случае Й ~ 1, когда для того, чтобы учесть все квантовые эффекты, необходимо непосредственное решение задачи (1)-(2), а само решение представляет собой медленно меняющуюся в области, где рассматривается процесс, функцию (см. на Рис. 1 справа).
Но часто возникает промежуточная ситуация, изображенная на Рис. 2. Здесь мы
Рис. 2: Промежуточная ситуация между глубокой квазиклассикой и чисто квантовым режимом: рассматриваемая система в проекции на фазовое пространство уже не является точкой, но и не покрывает всю область своего перемещения.
видим динамику волновых пакетов, проекции которых на фазовое пространство уже не точки, но и не покрывают всю область их перемещения. В этом случае квазиклассические методы оказываются слишком грубыми. Использование же основных квантовых методов, таких как разностные схемы, оказывается не самым эффективным, так как их универсальный характер не использует физическую особенность данного случая: большой объем точек рассматриваемого фазового объема, где решение все-таки мало отличается от пуля. Кроме того, и в квантовом случае общепринятые сеточные методы испытывают трудности при рассмотрении инфинитного движения, например, задач рассеяния, так как, как правило, используют фиксированную сетку, а при этом решение неограниченно перемещается в фазовом пространстве.
Целью настоящей работы является разработка нового приближенного метода решения задачи (1)-(2), который справился бы с описанными трудностями и удовлетворял бы следующим требованиям
Введение
применимость в промежуточном случае между глубокой квазиклассикой и чисто квантовым режимом. В этом случае использование квазиклассических методов не позволяет получить решение с требуемой точностью, а квантовые методы оказываются неэффективными;
являлся бы точным методом, т.е. параметры метода можно всегда выбрать таким образом, чтобы получить решение задачи с произвольной заданной точностью;
позволял бы эффективно рассматривать динамику волновых пакетов (в случае, например, инфинитного движения).
Краткое содержание работы. Описание метода
Итак, поставлена проблема построения приближенного метода решения задачи (1)-(2). С этой целью в настоящей работе используется модификация уже упомянутого выше метода гауссовых волновых пакетов, получившего широкое распространение при рассмотрении квазиклассического случая [32, ЗІ]. В частности нами используется результат работы [4], где эта методика была применена для аппроксимации функции Грина оператора Шрёдингера и была получена соответствующая оценка погрешности. Сразу же отметим, что всюду в настоящей работе задача (1)-(2) рассматривается в одномерном по пространству случае d = 1. Рассмотрение одномерного случая не является существенным ограничением, так как все конструкции допускают естественное обобщение на многомерный случай. В то же время появляется возможность наглядно и с достаточной строгостью описать новый метод и иллюстрировать его применение.
Изложение материала в диссертации разбито на три главы. В первой главе приведены необходимые сведения и результаты из теории гауссовых волновых пакетов. Во второй главе строится численный метод решения задачи (1)-(2) и дается его теоретическое обоснование. Работу завершает третья глава, где обсуждаются результаты расчетов новым методом некоторых модельных задач.
В методе гауссовых волновых пакетов уравнению (1) ставится в соответствие система обыкновенных дифференциальных уравнений
Введение
"% = ПО. «С) = ч,,
г/Л
f-'^WWW). 5(0) = 1,
^ = 5 P2{t) - V(Q(t)), 5(0) = 0.
Здесь начальные данные qa и ро произвольны, а о* > 0 - параметр, значение которого
будет дано ниже. На решении этой системы определяется функция
1 Ґ В(х — О)"2 Р БЛ
W{x, t\qo>Po) = 2{7Thy/2A1/2 ^{-^2М + VX * Q) + 'ft ) (9)
Когерентные состояния
При малых k, т.е. при переходе к квазиклассическому режиму, эта аналогия становится все точнее. Другим важным для нас свойством гауссовых волновых пакетов является то, что они образуют базис в пространстве состояний [9, стр. 158]. Строго это выражается формулами для прямого
Как уже было отмечено выше, метод гауссовых волновых пакетов (Gaussian wave packet dynamics, GWPD [34]) это метод приближенного решения нестационарного уравнения Шрёдингера, который основан на приближении искомой волновой функции гауссовым волновым пакетом с центром, движущимся по классической траектории. Поэтому совершенно естественно, что этот метод может применяться лишь в условиях, когда квантовомеханическое описание мало отличается от классического. Последнее имеет место в квазиклассическом режиме при малых h. Метод гауссовых волновых пакетов является наиболее развитым квазиклассическим методом, существует много модификаций этого метода (см. обзор [33] и ссылки там), которые применяются, например, при моделировании ядерных реакций. Так как метод, предлагаемый в настоящей работе, является модификацией метода гауссовых волновых пакетов, мы представим некоторые результаты и построения, принятые в этом подходе. Наиболее близким к нашему является изложение [30, 4].
В этом случае гауссов волновой пакет (1. 13) с центром, движущимся по классической траектории, уже не будет удовлетворять уравнению (1.16) точно. Но, как видно из формул (1.14), при переходе к пределу при h — О ширина пакета стремится к нулю. Это позволяет при построении приближенного решения рассматривать лишь два первых члена в разложении потенциала вокруг классической траектории. При такой локальной "квадратичной" аппроксимации потенциала погрешность будет возникать за счет отброшенных в разложении членов и может быть сделана сколь угодно малой за счет сужения пакета при h — 0. Но на практике значение параметра ft хоть и мало, но все же фиксированно и ширина пакета "отлична от нуля". Это приводит к возникновению погрешностей, которые возрастают по мере того как пакет уширяется с течением времени. Таков механизм построения приближенного решения уравнения Шрёдингера (1.16) в квазиклассическом случае методом гауссовых волновых пакетов. Формально он заключается в следующем.
Таким образом, выполнено условие (1.11): Л (0)В(0) + В (0)А{0) = 2сг, а начальные данные для Q и Р будем считать произвольными: (ffoiPo) Є R2. Первые два уравнения системы (1.17) образуют систему уравнений Гамильтона классического движения в потенциале V(x). Эта система имеет решения Q(t) и P(t) для любых начальных данных. Локальная квадратическая аппроксимация потенциала вблизи фиксированной траектории Q(t) есть следующий отрезок ряда Тейлора.
Поэтому функция p(x,t) из (1.13), построенная на решении системы (1.17), удовлетворяет уравнению Шрёдингера с гамильтонианом (1.19) точно. При подстановке же в общее уравнение (1.16), мы получим отличную от нуля невязку, которая позволяет судить о погрешности приближенного решения. Величина невязки (,2-норма) определяется остаточным членом в разложении (1.18) и складывается из двух составляющих. Первая - определяется величиной пространственной области вокруг центра разложения Q, где сосредоточен пакет (1.13). Чем меньше эта область, тем "точнее" рассматриваемое в ней приближение (1.18). Поэтому эта часть невязки может быть Глава 1. Метод гауссовых волновых пакетов неограниченно уменьшена за счет сужения пакета при h —» 0. Вторая составляющая определяется величиной остатка во всем остальном пространстве, где погрешность разложения (1.18) становится существенной. Но эта составляющая также мала (при определенных условиях), так как функция tp практически равна нулю в этой области (быстро убывает при удалении от центра разложения: как с т ).
Таким образом, в квазиклассическом случае при малых h функция (1.13) является приближенным решением уравнения Шрёдингера с произвольным потенциалом. Оценка погрешности имеет асимптотический характер: 0(ft1/f2),& —» 0. Механизм аппроксимации заключается в том, что при ft —» 0 пространственная ширина пакета неограниченно уменьшается4. Но, как видно из формул (1.14), на пространственную ширину можно влиять и с помощь параметра а. Для этого надо рассмотреть
Здесь за С обозначены постоянные, не зависящие от t, а и /г, конкретные значениях которых неважны. Сразу отметим, что это уточнение "искусственно" в том смысле, что имеет смысл лишь когда приближенное решение уравнения (1.16) рассматривается на малых отрезках времени: 0 t ст. Обычно же требуется находить решение на максимально больших временных промежутках. В этом случае мы по-прежнему будем иметь лишь 0(Ь}12). Но тем не менее, оценка (1.23) послужит большей наглядности дальнейших построений.
В этом разделе мы опишем модификацию метода гауссовых волновых пакетов, предложенную в работах [4, 5] и которая будет использована нами в дальнейшем при построении численной процедуры.
Функция ip{x,t) является приближенным решением уравнения (1.16), удовлетворяющим начальному условию специального вида — гауссову волновому пакету (1.1). Но это ограничение можно обойти, использовав преобразование Фурье-Гаусса (1.3), т.е. полноту системы когерентных состояний в ,2(R). Действительно, Глава 1. Метод гауссовых волновых пакетов.
Подводя итог, подчеркнем особенности настоящей модификации метода гауссовых волновых пакетов. Прежде всего все полученные оценки (1.30),(1.33) являются равномерными по ft и начальным данным. Стандартный метод гауссовых волновых пакетов основывается на асимптотическом представлении решения (1.21) при ft —+ 0. На практике значение ft хоть и мало, но все же фиксировано и границы применимости полученного приближенного решения, вообще говоря, не ясны и требуют дополнительного исследования. В данном же случае, как следует из (1.35), мы можем получить решение задачи с любой точностью на произвольном фиксированном временном интервале за счет выбора параметров метода (достаточно малого шага по времени). И для этого совсем не требуется малость ft. Глава 1. Метод гауссовых волновых пакетов.
В заключение скажем несколько слов о практической стороне. Нахождение приближенного решения непосредственно по формуле (1.34) требует при каждом применении оператора /С вычислять образ Фурье-Гаусса аргумента. Эту процедуру можно сделать более эффективной, воспользовавшись тем, что FBI-образ функции W может быть вычислен явно (см. также [4]), используя известную формулу для интеграла Пуассона (см. параграф 1.3)
Описанная в этом разделе модификация метода гауссовых волновых пакетов и полученная в результате аппроксимация функции Грина G будут использованы нами далее для построения численной процедуры решения уравнения (1.16) (задачи (1)-(2)). Но уже на данном этапе можно интерпретировать эту процедуру как приближение континуального интеграла по когерентным состояниям (Coherent State Path Integral, CSPI), представляющего функцию Грина задачи, многократным Римановым интегралом (1.37). Отметим, что в настоящее время проявляется большой интерес к применению численного континуального интегрирования к решению широкого круга задач во Глава I. Метод гауссовых волновых пакетов многих областях, особенно в квантовой и статистической физике. Примером может служить обзор [10]; CSPI рассматривается в [24].
Численный алгоритм
Итак, согласно теореме 2.2.1, мы можем с любой степенью точности заменить на всех слоях по времени А: интегрирование в (2.5) по всей плоскости интегрированием по некоторому фиксированному (для данного начального условия) прямоугольнику D. Прямоугольник D будем выбирать со сторонами, параллельными координатным осям. После введения в D равномерной сетки Q с шагами hq и hp, и приближенного интегрирования в (2.5) по квадратурной формуле прямоугольников, получим следующую численную процедуру решения задачи вторая степень матрицы Gft(r). При этом, согласно теореме 1.2.5, матрица G\ является аппроксимацией функции Грина с эффективным шагом 2г, которая является более точной, чем Gft(2r). Таким образом, после n-кратного возведения в квадрат мы будем иметь матрицу пропагатора с эффективным шагом 2пг. Посчитанная таким образом матрица может быть сохранена и использована в дальнейшем, когда величина начального шага т много меньше рассматриваемых отрезков времени, и слишком мелкий шаг неудобен.
Численная процедура (2.21), соответствующая расчетам на фиксированной сетке, подробно описана в [I]. Параметрами алгоритма, определяющими точность вычислений, здесь являются: шаг по времени г, FBI-параметр и, размер прямоугольника D и его расположение на плоскости, шаги сетки П. Эти параметры должны выбираться эмпирически. Результаты численных расчетов, полученные в [1], продемонстрировали влияние параметров на точность метода и доказали его практическую применимость.
Итак, описанная только что численная процедура решения задачи (1)-(2) на фиксированной сетке оказывается неэффективной, когда решением является волновой пакет, перемешающийся в фазовой плоскости (q,p). В этом случае (см. Рис.2 на стр. 16) требуется выбор изначально большого прямоугольника D, который смог бы "покрыть" все возможные положения пакета на рассматриваемом отрезке времени [0,Т]. При этом в каждый отдельный момент времени использование такой большой области неоправданно, т.к. решение отлично от нуля лишь на малом подмножестве D. Поэтому применение процедуры (2,21) потребует большого количества лишних вычислений и памяти. В связи с этим возникает необходимость в модификации алгоритма за счет введения подвижной сетки. Перейдем к изложению алгоритма.
Вместо прямоугольника D рассмотрим на каждом слое по времени к наименьший прямоугольник D(k) со сторонами, параллельными координатным осям, и содержащий множество всех точек плоскости, где значения решения по абсолютной величине превышают некоторый малый порог 6. Другими словами, 6 - параметр метода, определяющий степень точности, с которой значения решения считаются равными нулю. В прямоугольнике D{k) вводится сетка Q(k), состоящая лишь из "существенных" на текущем слое узлов сетки П. Вместо (2.21) используется следующая схема вычисления по слоям.
С течением времени прямоугольник D{k) может перемещаться в плоскости (gtp) и изменять свой размер. Соответствующие преобразования сетки П(к), связанные с добавлением и исключением узлов при переходе на следующий слой по времени, предлагается выполнять согласно следующему алгоритму шаг 1) Пусть известно решение Ф (&) на сетке Q(k). Вычисляем по формуле (2.24) решение Ф +1 на следующем временном слое во всех узлах Q(k); шаг 2) Анализируем найденные значения и удаляем лишние узлы из сетки. Для этого значения 4 k+l в узлах сетки Q{k) просматриваются вглубь по строкам или столбцам, в зависимости от выбранной стороны прямоугольника D{k). Из сетки удаляются строки (столбцы) узлов, в которых все значения решения меньше порогового по модулю: Ф +1(&) . ДО тех пор пока не встретится строка Глава 2, Построение численной процедуры (столбец), где есть хотя бы одно значение, превышающее 5. После этого указанная процедура повторяется, начиная со следующей стороны. Таким образом рассматриваются все стороны прямоугольника в определенной последовательности (например, по часовой стрелке, начиная с левого края); шаг 3) Рассматривается возможность добавления в сетку новых узлов. Сетка может быть расширена на сторонах, где прямоугольник не был обрезан. Если такие стороны есть, то они рассматриваются последовательно. При этом значения решения Ф +1 в потенциально новых узлах, не принадлежащих сетке Q(k), рассчитываются, как и прежде, по формуле (2.24). Сетка расширяется за счет строк (столбцов) узлов, которые добавляются с соответствующих сторон, до тех пор пока в них есть хоть одно значение решения, превышающее по модулю порог 6.
Иллюстрация работы описанного алгоритма представлена на Рис. 2.1. Результатом является приближенное решение Фь+1 на следующем слое по времени, определенное, вообще говоря, на новой сетке.
Алгоритм расчета на подвижной сетке описан. В заключение скажем несколько слов, поясняющих механизм работы предложенного метода. Как уже было отмечено в [18], способ построения функции G делает процедуру вычисления по слоям (2.4)-(2.5) практически идентичной классическому методу частиц. Действительно, согласно (2.2), перепишем формулу (2.24) (как общий случай (2.21)) перехода на следующий временной слой в виде.
Таким образом переход на следующий слой по времени можно разбить на несколько этапов. Сначала приближенное решение задачи (1)-(2) фк(х), і R на А ом слое по времени представляется с помощью преобразования Фурье-Гаусса суммой гаус-сианов (1.1) с центрами в узлах сетки П(Аг) — частиц. Затем осуществляется малый шаг по времени т, в продолжении которого каждая частица-гауссиан перемещается по классической траектории. В результате мы получим функцию которую можно интерпретировать как приближенное решение па следующем слое по времени, полученное методом частиц.
Равноускоренное движение частицы
При практических вычислениях по формуле (3.1) следует обратить внимание на следующее обстоятельство. Квадратный корень в (3.1) извлекается из комплексной величины. Но было бы ошибкой понимать под ним лишь одну из ветвей (в данном случае — аналитическое продолжение арифметического квадратного корня на комплексную плоскость с разрезом). Дело в том, что функция ФС1, как решение дифференциального уравнения, должна быть аналитической (здесь этот термин понимается в смысле гладкости) функцией. Поэтому под символом у/ в (3.1) надо понимать аналитическое продолжение на всю двулистную риманову поверхность квадратного корня. Это означает, что фаза подкоренного выражения может с течением времени превысить 2тт, что соответствует тому, что точка совершит полный оборот вокруг начала координат в комплексной плоскости и перейдет на другой лист. При этом мы должны "взять" значение другой ветви, чтобы функция изменялась непрерывно. Стандартные процедуры вычисления квадратного корня в языках программирования, как правило, рассматривают только одну ветвь, не учитывая этот "набег фазы". Это приводит к тому, что выражение (3.1) делает скачок. Во избежании этой ошибки необходимо специально отслеживать фазу подкоренного выражения, для того чтобы выбрать правильную ветвь корня.Движение частицы в гармоническом поле соответствует задаче (1)-(2) с потенциалом
Формула (3.6) вместе с выражением (3.1) для точного решения использовались для контроля точности.
Гармонический осциллятор является примером финитного движения, когда решение является волновым пакетом, который с течением времени перемещается в ограниченном объеме фазового пространства. Особенностью гармонического потенциала также является то, что решение сохраняет гауссов профиль. Причем, как легко видеть из (3.5) ширина гауссиана остается неизменной, если только в начальный момент времени она совпадала с шириной основного состояния 70 = ft/w, В представлении Фурье-Гаусса решение (по абсолютной величине) качественно ведет себя как гауссиан, центр которого вращается в фазовой плоскости (qt р) по эллипсу вокруг точки равновесия (0,0). Это проиллюстрировано на Рис. 3.1, где изображены носители приближенного решения в различные моменты времени. Результат на
Результаты расчетов с этими параметрами для различных шагов по времени приведены в Таблица 3.1. Задача решалась на отрезке [0,10] и каждый раз вычислялось соответствующее число слоев N. О величине погрешности можно судить по наибольшему отклонению: приближенного решения (кол, 3), среднего значения координаты (кол, 4) и нормы (кол. 5) от точных значений в узлах сетки. Из таблицы видно, что с ростом т погрешность увеличивается. Используя связь с методом частиц можно сказать, что частицы начинают "разлетаться" слигиком далеко друг от друга. Использование же слишком мелкого шага т также неоправданно, так как число слоев увеличивается, а каждый слой вносит некоторую "неустранимую" ошибку (квадратуры).
Повысить точность можно за счет "измельчения" частиц. В Таблица 3.2 приведены результаты расчетов с шагом т = 0.25 для различных а. Остальные параметры задачи и метода оставлены без изменений. Шаг т = 0.25 оказался слишком велик для а = Глава 3. Результаты расчетов, уменьшая а до 0.5 удалось сделать погрешность приемлемой. Дальнейшее уменьшение вызывает увеличение ошибки. Это связано с тем, что пространственная сетка с шагами hq = hp = 0.9 оказывается слишком грубой для все более "мелких" частиц. Это подтверждает Таблица 3.3, где приведены результаты расчетов для различных шагов hq и hp при фиксированных т = a = 0.25. Замечание. Для применения предложенного метода к решению задачи {1)-(2) необходимо требовать выполнения условия (1.27) на потенциал. Для гармонического осциллятора и некоторых других задач в настоящей работе, в которых V{x) неограниченно возрастает, формально это условие не выполнено. Но так как задача решается приближенно на конечном отрезке времени, мы можем считать потенциал в уравнении (1) равным V{x) лишь в достаточно большой, но ограниченной области, где фактически перемешается сетка и решение. За пределами этой области, где решение практически обращается в ноль, потенциал можно взять равным нулю.
Рассеяние на потенциальном барьере
Задача (1)-(2) с потенциалом (3,11) иногда рассматривается как модель колебательного движения молекул. Можно показать, что рассматриваемая система имеет лишь конечное число связанных состояний, а также получить для соответствующих собственных значений явные выражения через параметры потенциала [12]. На примере этой задачи можно пронаблюдать эффект возвращения, когда хорошо локализованное начальное условие в процессе своей эволюции может сложным образом "размазываться" по фазовому пространству и затем, спустя довольно большой промежуток времени Trtv (много больший характерных времен), снова локализоваться, вернувшись в исходное состояние. В иностранной литературе такое поведение решения называется quantum wave packet revival. Основную роль в теоретическом исследованииt этого феномена играет автокорреляционная функция. При этом, знание собственных значений, как в случае потенциала Морзе, позволят оценить время Trev Более подробное рассмотрение эффекта возвращения, его демонстрацию в численных расчетах простых модельных систем, а также обширный обзор литературы, можно найти в [52]. Параметры потенциала были взяты нами из работы [24], где эта задача также рассматривалась в качестве тестовой.
В этом случае дискретный спектр содержит 20 собственных значений, и время возвращения имеет порядок [24] TTev PS 125, что полностью согласуется с численными результатами, полученными в настоящей работе. Это демонстрирует Рис. 3.8, где приведен график квадрата модуля автокорреляционной функции и нормы решения. Глава 3. Результаты расчетов
Результат на Рис. 3.8 был получен с помощь численной процедуры (2.21) на фиксированной сетке, состоящей из 26 х 26 узлов при г = 0.01, о — 0.2, Легко видеть, что такой мелкий шаг по времени, который диктуется необходимостью аппроксимации функции Грина2, не удобен при решении задачи на рассматриваемых временах порядка TTev w 125. Поэтому, фактически, для перехода на следующий слой использовалась матрица пропагатора с эффективным шагом по времени 25т = 0.32, которая была получена с помощью последовательного применения процедуры возведения в квадрат матрицы пропагатора, описанной в конце параграфа 2.3.1. Получить наглядное представление о матрице с эффективным шагом позволяет Рис. 3,9, где сопоставлены модули матричных элементов пропагаторов с элементарным и эффективным шагом.
Рассмотрим задачу о взаимодействии частицы с барьером в виде максвеллиана высоты Vo и характерной шириной d Такой профиль барьера выбран потому, что метод [4] требует гладкости потенциала. На этот барьер будем бросать гауссов волновой пакет с различной энергией. Центр начального условия q$ выбирается в зависимости от d и GQ так, чтобы пакет и барьер изначально "не пересекались" (см. Рис. 3.10) Конкретные расчеты были проведены для двух различных значений энергии пакета1 Е = Ро/2 : Е = 3.0 V0 и Е = 6.0 VQ. На Рис. 3.11 изображена зависимость среднего значения координаты частицы от времени (x)(t). В первом случае (Е VQ) мы имеем отражение среднего значения от барьера, а во втором (Е V0) — прохождение через барьер. Результат на Рис. 3.11 был получен при следующих параметрах метода
Рис. 3.12 дает представление о том, как ведет себя приближенное решение с течением времени для случая (Е VQ). Здесь изображены начальное условие и приближенное решение в финальный момент времени і = 35. Мы наблюдаем "преобладание" отраженной волны над прошедшей. "Преобладание" здесь понимается в смысле коэффициентов отражения и прохождения, которые при использовании нестационарного подхода определяются обычно как вероятности обнаружить частицу по эту или ту сторону барьера соответственно [8]
3 Напомним, что всюду в расчетах в настоящей работе мы считаем, что параметр т, соответствующий массе.Приведены результаты расчетов для двух различных шагов по времени г = 0.1 и т = 1.0. В последнем случае точность вычислений значительно хуже. На Рис. 3.14 сопоставлены приближенные решения в финальный момент времени t = 35.0, полученные для различных т. Результат соответствующий большему шагу по времени, содержит ложные осцилляции, которые возникают из-за того, что "частицы", составляющие пакет, успевают разлететься слишком далеко друг от друга.
