Диссертация (Исследование релятивистских и нерелятивистских квантовых систем с помощью вычисления континуальных интегралов методом Монте-Карло), страница 2

Для описания актуальности данной темыбыло бы достаточно сослаться на Гинзбурга, включившего металлический водород в список важнейших проблем современной физики. Но будет правильнее указать два основныхнаправления, в которых экспериментальная и наблюдательная наука в настоящее времясталкивается с необходимостью теоретического описания данной системы. Во-первых, влабораторных экспериментах в самое недавнее время наконец была получена фаза, которая претендует на роль металлическго водорода. Тут надо отметить, что при высокихплотностях у водорода есть достаточно много фазовых состояний. В частности металлможет быть реализован как на основе решётки ионов H+, так и на основе решётке атомов(мы будем исследовать второй случай).

Другой сферой, которая сталкивается с даннымвеществом является астрофизика. Есть существенные аргументы в пользу того, что ядрапланет-гигантов в значительной части состоят из металлического водорода. Таким образом, для описания строения и эволюции таких планет, которых открыто большое количество и за пределами солнечной системы, необходимо исследовать его свойства.1.2Цель работыПроведение данной работы преследовало следующие цели:- Изучение существующих подходов к численному моделированию в физике конденсированного состояния, исследование их производительности, создание максимально эффективных методов и алгоритмов их реализации.- Построение и исследование численными методами моделей конденсированных сред,в первую очередь металлического водорода.

В данных моделях исследуются уравнениясостояния, а также фазовые состояния и переходы.1.3Научная новизна исследованияВ данной работе достигнуты существенно новые резултаты:- На основе теории численных методов и результатов тестовых вычислений построенывысокопроизводительные методы Монте-Карло моделирования в задачах физики конденсированного состояния. Также составлено подробное описание алгоритмов, позволяющеедругим исследователям быстро реализовывать их, не затрачивая силы на изучение теоретических основ и тонкостей их построения.

- Для вычисления физичесчки важных величинпостроены где это возможно соответствующие "решёточные наблюдаемые"и алгоритмыих расчёта. Также разработаны и реализованы методы расчёта некоторых физическихпараметров, на выражающихся черех решёточные наблюдаемые напрямую. - С прменением вышеуказанных методов исследованы модели физики конденсированного состояниявещества. В частности, для металлического водорода при высоких плотностях получены уравнения состояния, обнаружении изучен фазовый переход "кристалл-жидкость".Данное исследование позволило получить описание данной среды в широком диапазонетемператур и плотностей.2Вычисление интеграла по траекториям методом МонтеКарло (PIMC)Этот раздел содержит краткое введение в формализм и методику PIMC, более подроб-но см.

[17] и [3].2.1Метод Монте-Карло для интеграла по траекториямРассмотрим систему, описываемую координатами x. В общем случае: в трёхмерном пространстве для системы, состоящей из Np частиц – таких координат будет 3Np штук. Временна́я эволюция системы определяется гамильтонианом H. При конечной температуре,которую для удобства будем описывать параметром β :β = 1/kB T,(1)где kB – постоянная Больцмана, T – температура в кельвинах.Матрица плотности такой системы выглядит следующим образом (здесь и далее мыработаем в координатном представлении):ρx0 →xNt = hx0 |e−βH |xNt i.Соответствующая (2) функция рапределения имеет вид:ZZ = trρ = dx0 hx0 |e−βH |x0 i.Среднее значение наблюдаемой A вычисляется по формулеZ11hAi = tr(Aρ) =dx0 hx0 |Ae−βH |x0 i.ZZ(2)(3)(4)Отметим, что в квантовой теории как правило используется нормированная матрицаплотностиρn =ρρ= .trρZ(5)В таком случае очевидно trρn = 1, а формула для среднего значения наблюдаемой принимает видhAi = tr(Aρn )(6)Перейдём к формулировке теории на языке интеграла по траекториям.

Для этоговведём "временно́й шаг"τ и (целое) число шагов Nt , связанные соотношениемβ = Nt τ.(7)Разложим матрицу плотности в произведение Nt матриц плотностиρx0 →xNt = ρx0 →x1 . . . ρxt−1 →xt . . . ρxNt −1 →xNt ,(8)или, в координатном представленииhx0 |e−Nt τ H |xNt i = hx0 |e−τ H |x1 i × . . . × hxt−1 |e−τ H |xt i × . . . × hxNt −1 |e−τ H |xNt i.(9)Через эти промежуточные матрицы определяется "действие шага"St , по определению равноеρxt−1 →xt = hxt−1 |e−τ H |xt i ≡ e−St .(10)"Решёточное действие"S по определению равноS=NtXSt .(11)t=1Это разложение, получаемое при помощи формулы Троттера, полностью справедливотолько при Nt → ∞, но для практических целей численного моделирования Nt можетбыть выбрано достаточно большим, чтобы от него не зависел результат.

Далее введёмобозначениеDx =NtYdxt .(12)t=1Тогда, подставив (9) в формулы (3) и (4), можно их представить в следующем виде:ZZ = Dxe−S ,(13)RZDxAe−SDxe−SRA=hAi = RDxe−SDxe−S(14)Формула (14) раскрывает основную идею метода PIMC. Пусть у нас есть (достаточно большой) набор траекторий x = x0 . . . xt . . . xNt , где вероятность того, что некотораятраектория будет включена в данный набор пропорциональна её "статвесу"π(x) ∼ e−S(x) .(15)Среднее значение любой наблюдаемой может быть получено простым (арифметическим)усреднением по данному набору.2.2Интеграл по траекториям при нулевой температуреВ случае нулевой температуры (обычная квантовая механика) можно построить абсолютно аналогичный формализм. При этом некоторые исходные положения и их физическийсмысл отличаются, но используются те же основные идеи.

В фейнмановской интерпретации квантовой механики частица описывается интегралом по траекториямZ−iH(tN −t0 )/~iSm (x)/~ hx0 |e|xN i =Dxe x(t0 )=x0 , .(16)x(tN )=xNЗдесь Dx = lim dx1 . . . dxN −1 ,n→∞"Zmdt2Sm (x) =dxdt#2− V (x)(17)– действие в пространстве Минковского. Средние значения наблюдаемых величин вычисляются следующим образом:RhAi =DxA(x)eiSm (x)/~R.DxeiSm (x)/~(18)По причинам, которые будут пояснены в следующем параграфе, в этих формулах необходимо, сделав виков поворот t → −iτ, перейти в Евклидово пространство.

После этогоуже становится очевидно совпадение формализма континуального интеграла в квантовоймеханике при нулевой температуре с формализмом при конечной температуре:ZZ = Dxe−S(x)/~ ,1hAi =ZZDxA(x)e−S(x)/~ ,(19)(20)гдеZS(x) ="mdτ2dxdτ2#+ V (x)(21)– Евклидово действие.Следующим естественным шагом является переход от бесконечномерных интегралов кконечномерным, проще говоря, дискретизация. Будем пользоваться равномерной решеткой по времени τn = na, (a – шаг решетки, n = 0, . . .

, N − 1.). Набор {x(τn ) ≡ xn } ≡ xназывается траекторией или конфигурацией. Граничные условия как правило (и всюдув данной работе) выбираются периодические:xN = x0 . Решеточное действие выглядитследующим образом:Nstep −1 S(x) =Xn=0m (xn+1 − xn )2+ V (xn ) a.2a2(22)Перегруппировкой членов или интегрированием по частям в (21) его можно привести квидуNstep −1 S(x) =Xn=0Nstep −1Xm (−xn )(xn+1 − 2xn + xn−1 )+V(x)a=[K(xn ) + V (xn )] .n2a2n=0(23)Выражения (19), (20) после дискретизации переходят соответственно вZ=Xe−S(x)/~ ,(24)confhAi =Pгде1 XA(x)e−S(x)/~ ,Z conf(25)обозначает сумму по всем конфигурациям (которых бесконечное число).confТеоретически, набор формул (23), (24), (25) позволяет производить вычисление средних: для этого надо взять достаточно большой (чтобы смоделировать бесконечную суммуP) набор конфигураций и применить к нему эти формулы.

В реальности действие наconfразных конфигурациях различается во много раз, а их статистический весπ(x) ≡1 −S(x)/~eZ(26)– на много порядков. Поэтому для получения сколько-нибудь точного ответа требуетсяогромное количество конфигураций и, следовательно, огромное время. Эту проблему решает метод существенной выборки.2.3Метод существенной выборки. Термализация.Формулы для вычисления средних (20),(25) можно записать так:RZDxA(x)e−S(x)/~Dxe−S(x)/~R=dµA(x),dµ=,hAi = RDxe−S(x)/~Dxe−S(x)/~A(x)e−S(x)/~Xe−S(x)/~P −S(x)/~ =A(x)P(x), P(x) = P −S(x)/~ .eeconf(27)PhAi =confconf(28)confВидно, что с математической точки зрения вычисление средних сводится к интегрированию с некоторой мерой (27) или к взятию среднего с некоторой вероятностью (28).Надо отметить, что величину P(x) можно интерпретировать как вероятность именно благодаря тому, что в излагаемой теории был совершен переход к евклидовому пространству.Формула (28) раскрывает идею метода существенной выборки: для вычисления среднегозначения некой величины достаточно взять среднее арифметическое по набору конфигураций, в котором вероятность появления конфигурации в наборе пропорциональна еестатистическому весу:hAi =1NconfNconfXA(xk ),(29)k=1P(xk ) = π(xk ) ∼ e−S(xk )/~ .(30)Метод получения таких специфических наборов конфигураций основан на свойствахмарковских случайных процессов.

В марковском процессе вероятность P (x → x0 ) переходаиз одного состояния x в другое x0 определяется только начальным и конечным состоянием.Как известно, такой процесс обладает предельным распределением (так как нас интересуетполучение распределения (30), будем обозначать его той же буквой), которое не зависито начального состояния:lim Pn→∞(n)Z0(x → x ) ≡ limn→∞Zdx1 . .

.dxn−1 P (x → x1 ) . . . P (xn−1 → x0 ) = P(x0 ), ∀x.(31)Чтобы получить заданный вид предельного распределения P(x0 ), надо соответствующим образом подобрать вероятность перехода P (x → x0 ). Достаточным (но не необходимым) требованием к ней является условие детального баланса:P (x → x0 )P(x) = P (x0 → x)P(x0 ).(32)Данное условие не фиксирует вероятность перехода однозначно. Конкретный вид P (x →x0 ), хотя и не влияет на предельное распределение, но радикально влияет на скорость иточность вычислений. Кроме того, зависит время достижения предельного распределениязависит как от вероятности перехода, так и от начального состояния.Таким образом, необходимые для вычислений наборы конфигураций, удовлетворяющие условию (30), получаются в результате термализации – марковского случайного процесса, удовлетворяющего требованию детального баланса:P (x → x0 )π(x) = P (x0 → x)π(x0 ).(33)В следующих параграфах будут рассмотрены некоторые наиболее распространенные иэффективные алгоритмы, основанные на различных способах выбора P (x → x0 ).3Вычислительные алгоритмы3.1Алгоритмы генерации конфигурацийВсе нижеизложенные алгоритмы можно считать обобщением или частным случаем методаМетрополиса-Гастингса.

В его основу положено разложение вероятности перехода[3]Z0000P (x → x ) = T (x → x )A(x → x ) + δ(x − x ) 1 − dyT (x → y)A(x → y) ,(34)гдеT (x0 → x)π(x0 ).A(x → x ) = min 1,T (x → x0 )π(x)0(35)Легко проверить, что (34) удовлетворяет условию (33) для любого T (x → x0 ). Смысл формулы (34) следующий: генерируется новая пробная конфигурация при помощи пробнойвероятности T (x → x0 ), затем с вероятностью принятия A(x → x0 ) эта конфигурацияпринимается (x0 добавляется в набор), с вероятностью 1 − A(x → x0 ) конфигурация отклоняется (в набор добавляется еще одна копия x).Конкретный вид алгоритма типа Метрополиса-Гастингса определяется выбором пробной вероятности T (x → x0 ).