CC BY
21
Монте-Карло моделирование металлического водорода: фазовый переход и уравнение состояния
А. А. Новоселов1,0, О. В. Павловский2,6, М. В. Улыбышев2
1 Московский государственный университет имени М. В. Ломоносова, физический факультет, кафедра квантовой теории и физики высоких энергий.
2 Институт теоретических проблем микромира МГУ.
Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.
E-mail: а novoselov@goa.bog.msu.ru, 6 ovp@goa.bog.msu.ru Статья поступила 02.10.2013, подписана в печать 06.11.2013.
Проведено численное моделированию атомарного (металлического) водорода методом PIMC (Монте-Карло для интеграла по траекториям). Исследован диапазон температур и плотностей, в котором поведение электронов определяется квантовой статистикой, а протонов — классической. В этой области получены уравнения состояния в виде зависимостей внутренней энергии и давления от температуры и плотности. Эти зависимости позволяют выявить и исследовать фазовый переход между кристаллом и жидкостью.
Ключевые слова: метод Монте-Карло, металлический водород, уравнение состояния.
УДК: 538.9. PACS: 62.50.-p, 61.20.Ja, 64.30.-t, 67.80.F-.
Введение
Статья посвящена моделированию металлического водорода при помощи вычисления интеграла по траекториям методом Монте-Карло (Path integral Monte Carlo — PIMC) [1]. В исследуемом диапазоне температур и плотностей электроны образуют вырожденный ферми-газ, в то время как ядра могут быть описаны классической статистикой. Эти два факты в совокупности позволяют на самом раннем этапе избежать известной «проблемы знака» для фермионов. Цель исследования — получение зависимостей внутренней энергии и давления от температуры и плотности. Также изучается фазовый переход кристалл-жидкость в широком диапазоне плотностей.
Изучение металлического водорода представляет интерес как с точки зрения астрофизики, так и в связи с прогрессом в экспериментах на алмазных наковальнях (diamond anvil cell), в которых недавно был получен металлический водород [2, 3].
Мы изучаем зависимость свойств атомарного водорода от температуры и давления, описываемых соответственно параметрами ß
ß = 1/kßT (1)
и rs (радиус Вигнера-Зейтца):
Р = А • (2)
з nrS3
Моделируется область конечного объема, содержащая Np = 128 частиц. В процессе моделирования вычисляются внутренняя энергия (как сумма кинетической и потенциальной энергии частиц)
E = K + V (3)
и давление P. Как известно, функции E(p, T) и P(p, T) обеспечивают полное термодинамическое описание системы. В качестве естественной меры упорядоченности системы мы вычисляем отношение Линдеманна
Л =
Здесь (х2) — смещение частицы относительно ее положения в кристаллической решетке, а — расстояние между ближайшими соседними частицами. Отношение Линдеманна позволяет явно различить кристаллическую и жидкую фазы.
В работе часто используются «ядерные» единицы измерения: ке = Н = е = тр = 1; здесь ке — постоянная Кулона, а тр — масса протона. Соответствующие единицы основных физических величин следующие: ядерный боровский радиус а0м = = 2.88199 • 10~14 м для расстояния; ядерный хартри На = Ем = 8.00516 • • 10-15 Дж для энергии; рм = 3.34420 • 1026 Па для давления; рм = 6.98747 • 1013 кг/м3 для плотности и Тм = 5.79811 • 108 К для температуры. В ядерных единицах масса электрона составляет те = 5.44617 • 10~4, а (электронный) боровский радиус а0е = 1/те = 1.83615 • 103.
1. Модель
Гамильтониан модели металлического водорода запишем в виде
Hfull = Kn + Ke + Vo + V + V,t.
(5)
Здесь Км и Ке — кинетическая энергия ядер (протонов) и электронов соответственно; У0, Уе и — потенциальная энергия взаимодействия межъядерного, электрон-электронного и ядер с электронами соответственно; все они представляют собой сумму парных кулоновских взаимодействий, например:
np ¿1-1
vo=Е£
¿1=1 ¿2=1
1
¿1 ¿2
(6)
Rn
(4)
Существует достаточно широкий диапазон температур и плотностей, в котором, с одной стороны, электроны можно рассматривать как вырожденный ферми-газ и использовать модель Томаса-Ферми [4] (т. е. квантовый характер их статистики имеет определяющий
характер), но, с другой стороны, протоны заведомо невырождены и их статистика не важна. В таких предположениях мы можем рассматривать только протоны, причем применять для них классическую (больцма-новскую) статистику. Учет электронов приводит к то-мас-фермиевской экранировке взаимодействия [5]. Таким образом, эффективный гамильтониан [6]
Н\и\\ = Км + Ум.
(7)
Здесь Ум — потенциальная энергия протонов с экранированным взаимодействием:
N ¿1-1
Ум=ЕЕ
¿1 = 1 ¿2=1
ехр{-гМ2/йтр }
(8)
Томас-фермиевский радиус экранировки ЛтР определяется по формуле
= V 12 .
(9)
Гамильтониан (5) может быть быть сведен к (7) при следующих условиях. Во-первых, ядерные (сильные) взаимодействия между протонами должны быть пренебрежимо малы, т. е. расстояние между ними (приблизительно равное ) должно быть много больше их размера Яр. Во-вторых, теория Томаса-Ферми должна быть применима к электронам: для этого внутри радиуса экранировки должно находиться большое количество электронов. Отсюда очевидно следует, что расстояние между ядрами должно быть меньше (электронного) боровского радиуса. Таким образом, наши приближения применимы для плотностей, соответствующих
К < Г < аое.
(10)
Соответствующие ограничения в ядерных единицах
и в СИ следующие: « 3 • 10
9 - 10
-16
ми
а0е « 2 • 103 « 5 • 10 11 м. Приведем также оценки для предельных плотностей (2): ~ 3 • 103 кг/м3 и Ртах « 6 • 1017 кг/м3.
Кроме этого наше приближение корректно, если электроны вырождены, а протоны — нет. Температуру вырождения можно оценить как «тг2. Таким образом, допустимый диапазон температур зависит от давления и определяется соотношениями
тег2 < в < г2.
(11)
Предельные температуры в ядерных единицах следующие: втт ~ 5 • 10-4 т1 и втах ~ Г2. Отсюда получаются следующие оценки при заданных плотностях (в единицах СИ): Ттп « 0.9р2/3 кг-2/3м2К и Ттах « 2 • 105р2/3 кг-2/3м2К.
2. Вычисление интеграла по траекториям методом Монте-Карло
2.1. Формализм и методика PIMC
Этот раздел содержит краткое введение в формализм и методику Р1МС; более подробно см. работы [7, 8].
Рассмотрим систему с гамильтонианом Н и соответствующей ему матрицей плотности рх0^щ = = {х0\е-вН\хМ1}, описываемую координатами х, при
температуре в. Перейдем к мнимому времени и, введя «временной шаг» т, по определению равный
1/Т = в = N т, (12)
разложим матрицу плотности в произведение N матриц плотности:
N-1 N-1
{Х0\е-вН\хм,} = П х\е-тН\х+1} = Д е-5' = е-5
t=0
'=0
(13)
где 5 — «решеточное действие». Это разложение, получаемое при помощи формулы Троттера, полностью справедливо только при N — то, но для практических целей численного моделирования N может быть выбрано достаточно большим, чтобы от него не зависел результат.
Далее вводя обозначение Т>х = 1 ^х', запишем формулы для функции распределения и среднего значения наблюдаемой А:
Z = 1г р =
йхо{х0\е вН\х0} =
хе
5
{А} = 1 \х(Ар) = 1 = / 'ВхАе-5
йх0ЫАе-вН\х0} =
¡Vxe
5
А
хе
5
5
(14)
(15)
Формула (15) раскрывает основную идею метода Р1МС. Пусть у нас есть (достаточно большой) набор траекторий х = х0,..., х',..., хN, где вероятность того, что некоторая траектория будет включена в данный набор, пропорциональна ее «статвесу»
п(х) - е-5(х\ (16)
Среднее значение любой наблюдаемой может быть получено простым (арифметическим) усреднением по данному набору.
Способ получения набора траекторий с правильным (16) распределением основан на свойстве марковских цепей сходиться к некоторому предельному распределению. Для того чтобы марковская цепь с вероятностью перехода V(х — х') сходилась к предельному распределению п(х), достаточно выполнения условия детального баланса
V(х — х')п(х) = V(х' — х)п(х'). (17)
Обобщенный алгоритм Метрополиса-Гастингса основан на разложении
V(х — х')= Т(х — х')А(х — х')+ (18)
+ д(х - х'){ 1 -
йуТ (х — у)А(х — у)
А(х — х') = тт
1,
Т (х' — х)п(х')
(19)
Т (х — х')п(х)
Такая вероятность перехода удовлетворяет условию детального баланса для любого Т(х — х'). Формула (19) означает следующее: генерируется пробная траектория с вероятностью Т; она добавляется в набор с вероятностью А или происходит ее возвращение к предыдущей траектории и добавление в набор еще одной ее копии с вероятностью 1 - А. Во-первых, распределение Т(х — х') должно быть близко к п(х');
г
1 2
во-вторых, должен существовать алгоритм быстрой генерации соответствующих случайных чисел.
В примитивном алгоритме переход от старой траектории к новой состоит в попытке изменения координат на одном слое по мнимому времени I. Для больших систем и большого числа временных слоев такой алгоритм приводит к гигантской автокорреляции, в результате требуется много вычислительного времени для получения действительно статистически независимых траекторий. Эта проблема решается применением многоуровневого алгоритма [8]. Он основан на быстрой генерации траектории грубого приближения, которая увеличивает вероятность принятия более точной траектории в дальнейшем.
Рассмотрим многоуровневый алгоритм с делением пополам. Возьмем участок траектории длиной 2Ыеуе1 слоев, например 5 = (х0,..., х2«1еуе1). Эта часть траектории делится на уровни: нулевой 50 = (х0, х2«1еуе1) (не будет изменяться в течение данного многоуровневого хода), первый = (х2ы еуе1 -1), второй 52 = (х2«1еуе-2,х2«1 еуе1 -1 + 2«1 еуе1 -2) и т.д. Введем «действие уровня» пк($к) = пк($0,..., 5к—1, 5к), которое является функцией 5к, а координаты на предыдущих уровнях являются параметрами. Действия на промежуточных уровнях могут быть выбраны произвольно, но действие последнего уровня должно быть равно настоящему решеточному действию:
ПЫеуе1(5Ы еуе^ = ^^
(20)
Алгоритм типа Метрополиса-Гастингса с распределениями вероятностей
7к($к) = 7к(..., «к— 1; $к; $¿+1,... ^ ..., ; $к; ...),
(21)
Ли(5'к) = Л(..., 1; $¿+1,... ^ ..., 1; $к; 5*+ь ...) =
= Ш1П
Тк(5к)пи(5'к )пк—1(5к) Ти(5'к )пк(5к)пк—1(5'к )]
(22)
удовлетворяет уровневому условию детального баланса на каждом уровне
Ри($'к )-
пи($и)
= Гк(5к)~
пк(5к )
Пк— 1 ($к— 1) Пк_1(5к—1)'
что ведет к полному детальному балансу (17).
(23)
нием между ними
\
¿>?1(0 - х«2(0 + Ьпа)2. (25)
а=1
Решеточное действие, соответствующее гамильтониану (7) с потенциальной энергией (8), имеет вид
- 1 п п = 5 = Бт + ,
(26)
^ = £ £ £ (х?(*) - хП*-[)+(°)2, (27)
t=1 ¿=1 а= 1 2Т
ДД^ 1 ехр-"^3 (О/ад ^=£££ £ .п^пзт. (28)
t=1 ¿1 = 1 ¿2 = 1 п1,2'3(Л = _ 1 ¿1 ¿2 УЧ
11 ¿2
Будем использовать многоуровневый алгоритм. Действие уровня можно выбирать произвольно, но существует оптимальный выбор — интеграл полного решеточного действия по координатам следующих уровней:
Пк($к) =
d5k+1 ... d5Nlеуе1 п($).
(29)
Этому условию удовлетворяет следующий алгоритм. В качестве пробной вероятности возьмем кинетическую часть действия
Т(х^(1)\п^ + 1) - n(t)) ~ ехр |- (х°° - Х0^ ,
(30)
где введены обозначения х0 = x(t+1)+x(t-1)+¿ИШЬ^», xто(t) = х! ({) + Ln/(t). Соответственно требуется лишь генерация распределения Гаусса (преобразованием Бокса-Мюллера) и определение п'(0, п+ 1) и х'({) за счет условий -1/2 < х'({) < 1/2 и п'^ + 1) - п'(t) = п^ + 1) - n(t). На каждом уровне используется временной шаг т ^ тк = 2Ыеуе1—кт, а число намоток сохраняется пк^ + 1) - пк({) = пк— 1(t). Соответствующая вероятность принятия
Лк ($к) = т1п
е
(5')
е—5у ($к)
(31)
($к) определяется соотношением (28), в котором первая сумма берется только по временным слоям уровня 5к, т в данном случае — настоящий временной шаг.
2.2. Подробности моделирования
Для моделирования используется кубическая ячейка с периодическими пространственными граничными условиями. Размер ячейки
Ь = 3/^ /V (24)
а = х, у, г координаты I -й частицы на t -м временном слое обозначим х®(0. Для полного описания системы также требуются «числа намоток» па(0 = -1,0, 1, которые показывают «перескок» траектории с одной границы ячейки на противоположную через периодические пространственные граничные условия. Потенциальная энергия взаимодействия частиц определяется расстоя-
3. Результаты
Вычисления проводились для Ыр = 128 частиц при следующих значениях параметров: Ыа = 1/в от 0.5 • 10—5 до 4.75 • 10—5 с шагом 0.25 • 10—5 и дополнительными точками 0.57 • 10—5 и 0.66 • 10—5; г5 изменялся от 200 до 450 с шагом 50. Эти значения параметров в ядерных единицах соответствуют значениям температуры от 2.9 • 103 до 27.7 •103 К и плотности от 183 • 103 до 2085 • 103 кг/м3. Сетка точек, в которых производились вычисления, показана на рис. 1. Во всех отмеченных точках получены значения давления, кинетической и потенциальной энергии и отношения Линдеманна. В табл. 1, 2 представлены полученные значения при максимальной и минимальной исследованной плотности.
Рис. 1. Фазовая плоскость
Таблица 2
Параметры системы при р = 183 • 103 кг/м3, Ур = 20 686.0 кДж/моль
Т, К, У - Ур, Е - Ур, Р,
кК МДж/моль МДж/моль МДж/моль ППа
11.6 311.5(9) 392.0(1) 703.4(9) 4.3908(1)
10.2 304.1(8) 385.1(1) 689.1(8) 4.3894(1)
8.7 296.4(8) 378.5(1) 674.8(8) 4.3879(1)
7.2 292.9(7) 371.4(1) 664.2(7) 4.3869(1)
5.8 292.8(6) 327.5(1) 620.3(6) 4.3835(1)
4.3 290.8(5) 324.1(1) 614.9(5) 4.3831(1)
3.9 290.6(5) 323.4(1) 614.1(5) 4.3830(1)
3.3 290.4(5) 322.9(1) 613.3(5) 4.3829(1)
2.9 289.8(4) 322.9(1) 612.7(4) 4.3829(1)
Т а б л и ц а 1
Параметры системы при р = 2085 • 103 кг/м3,
уР = 117.678 МДж/моль
Т, К, У - Ур, Е - Ур, Р,
кК МДж/моль МДж/моль МДж/моль ППа
27.7 1.035(1) 2.089(1) 3.124(1) 262.010(2
26.1 1.028(1) 2.085(1) 3.113(1) 261.997(2
24.7 1.026(1) 2.080(1) 3.105(1) 261.988(2
23.2 1.020(1) 2.076(1) 3.096(1) 261.975(2
21.7 1.014(1) 2.072(1) 3.086(1) 261.964(2
20.3 1.008(1) 2.068(1) 3.076(1) 261.950(2
18.8 1.005(1) 2.065(1) 3.070(1) 261.942(2
17.4 1.001(1) 2.061(1) 3.062(1) 261.932(2
15.9 0.999(1) 2.057(1) 3.057(1) 261.925(1
14.5 0.993(1) 2.054(1) 3.048(1) 261.912(1
13.1 1.018(1) 1.967(1) 2.985(1) 261.899(1
11.6 1.019(1) 1.964(1) 2.983(1) 261.899(1
10.2 1.019(1) 1.961(1) 2.980(1) 261.894(1
8.7 1.021(1) 1.959(1) 2.980(1) 261.895(1
7.2 1.019(1) 1.958(1) 2.977(1) 261.891(1
Мы видим, что свойства системы зависят от плотности значительно сильнее, чем от температуры. Внутренняя энергия почти полностью состоит из потенциальной энергии, определяемой расстоянием между протонами, т. е. плотностью. Чтобы выделить значимую часть потенциальной энергии, введем Ур — аналог энергии Маделунга. Она равна потенциальной энергии идеального кристалла при абсолютном нуле с юкавовским взаимодействием между узлами
^ ¿1 1 ехш — г~. / ктс у
(32)
Ур = ЕЕ
ехрЬ^/ВД
1 =1 2 =1 1 2
где г0 — положения частиц в узлах объемно-центрированной кубической решетки.
3.1. Энергия, давление, уравнение состояния
Средняя внутренняя энергия {Е} вычисляется по формуле (3) с учетом: (12), (27), (28):
Е = К + {У}-(^ - в
М I).
(33)
На рис. 2 показана потенциальная энергия У как функция температуры при плотности 2085 • 103 кг/м3 (г = 200). Мы наблюдаем небольшой рост с ростом температуры и резкий скачок при определенной температуре, указывающий на фазовый переход. На рис. 3 изображена кинетическая энергия при данной плотности. Она также растет с температурой, но при фазовом переходе имеет скачок вниз. При более низких плотностях этот скачок исчезает, и остается только излом (скачок производной). Таким образом, зависимость выглядит как фазовый переход первого рода при плотностях 261 • 103 кг/м3 и меньше и как фазовый переход второго рода при плотностях 618 • 103 кг/м3 и больше.
V, МДж/моль
119.76
10 15 20 25 Т, кК
Рис. 2. У(Т) при р = 2085 • 103 кг/м3 (г = 200)
Рис. 3. К(Т) при р = 2085 • 103 кг/м3 (г$ = 200)
Наблюдаемая для давления имеет вид
р > = 31? (<К>- г( £ щ ) (34)
На рис. 4 изображена зависимость давления от температуры при плотности 2085 • 103 кг/м3 (г5 = 200). Она испытывает скачок при той же температуре, что и энергия, что подтверждает существование фазового перехода. Как и энергия, давление в основном зависит от плотности и незначительно меняется с температурой, именно этого следовало ожидать в конденсированной среде. Недостатком этого свойства является сложность численных термодинамических расчетов, так как частные производные по температуре и плотности, встречающиеся в уравнениях, различаются на порядки величин. Формально функция Р(р, Т) позволяет построить изобары, но разрешения численных данных для этого недостаточно, несмотря на большое количество точек на фазовой плоскости.
Р, ППа
Рис. 4. Р(Т) при р = 2085 • 103 кг/м3 (г* = 200)
3.2. Фазовый переход
Отношение Линдеманна (4) — хорошая мера неупорядоченности кристаллической решетки, поэтому его можно использовать для обнаружения фазового перехода. На рис. 5 показано отношение Линдеманна для всего исследованного диапазона плотностей и температур: ясно виден фазовый переход. Пологий участок этих зависимостей в твердой фазе (слева внизу) содержит некоторую физическую информацию о неупорядоченности кристаллической решетки. Пологий участок в жидкой фазе (справа вверху) — артефакт конечного объема, его уровень зависит только от размера ячейки моделирования.
Л
Рис. 5. Отношение Линдеманна при различных плотностях
Определив положение фазового перехода, построим фазовую диаграмму для металлического водорода (см. рис. 1). Тут надо сделать существенное замечание. Наблюдаемые в методе Р1МС являются средними по набору траекторий, получаемых в результате термализа-ции, начатой с некоторой траектории. С некоторого момента начинают получаться равновесно распределенные конфигурации, но неизвестно, как скоро это произойдет. Также известно, что модели систем в окрестностях фазового перехода обычно плохо термализуются, особенно через этот переход. Например, мы начинаем с идеального кристалла при абсолютном нуле. Он быстро термализуется в некоторое твердое состояние, которое кажется устойчивым, но может впоследствии дотермализоваться в жидкое. Пример показан на рис. 6. Получение физических траекторий занимает непредсказуемо много вычислительного времени. Таким образом, полученное положение фазового перехода формально является его верхним пределом. Нижний предел можно получить моделированием, начиная с «жидкой фазы», но не стоит ожидать, что он будет сильно отличаться от полученного нами верхнего предела. Кроме того, оказывается, что термализация из жидкого в твердое состояние происходит так медленно, что ее вряд ли можно получить за приемлемое вычислительное время.
Другой способ определения положения фазового перехода заключается в термодинамическом интегрировании свободной энергии [9]. Это позволяет избежать
nconf
Рис. 6. Термализация отношения Линдеманна при р = 188.4 • 106 кг/м3, T = 13.1 кК (кривая 1) and T = 14.5 кК (кривая 2)
долгой термализации из жидкого состояния в твердое. Но этот подход основан на интегрировании данных моделирования, содержащих некоторую погрешность, на большом промежутке.
Наконец, надо пояснить, почему мы утверждаем, что фазовый переход происходит именно между фазами жидкости и кристалла с объемно-центрированной кубической решеткой. Это можно легко увидеть, просто нарисовав соответствующие траектории. Примеры траекторий в обеих фазах приведены на рис. 7, 8. Существование жидкое (не кристаллического) состояния над переходом также подтверждается отношением Линдеманна, которое в этой области много больше единицы (см. рис. 5). Объемно-центрированная кубическая кристаллическая структура под переходом видна по характерным траекториям (рис. 8).
10001
500-
-500
-1000 -800
800
Рис. 7. Траектория в жидкой фазе
800
Рис. 8. Траектория в фазе кристалла с объемно-центрированной кубической решеткой
Заключение
Было проведено моделирование металлического атомарного водорода методом Р1МС. Были исследованы его термодинамические характеристики в широком диапазоне температур и плотностей, получены численные уравнения состояния. Обнаружен и исследован фазовый переход между жидкой и кристаллической фазами.
Для рассматриваемой системы известны основные термодинамические параметры: температура, плотность, давление и энергия. Однако, не измерена ее энтропия, что будет предметом дальнейших исследований. Методы расчета энтропии и построения адиабат несколько сложнее, чем, например, для изотерм, потому что для энтропии нет Р1МС-наблюдаемой. Таким образом, приходится решать термодинамические дифференциальные уравнения. Формально они дают полную информацию о системе, так как нам известны функции Е(р, Т) и Р(р, Т), но получить эту информацию численно непросто. Для достаточно точного вычисления производных по температуре и давлению и последующего интегрирования сетка расчетных точек должна включать близкие точки на большом промежутке, что требует большого объема вычислений. Также мы планируем применить альтернативный способ вычисления производных, основанный на построении наблюдаемых непосредственно для них.
Еще одна цель дальнейшей работы — осуществление термализации из жидкого состояния в твердое для определения нижнего предела для фазового перехода. Ожидается, что это может быть сделано относительно быстро, если начинать термализацию с двухфазной системы. Это позволит не только более точно определить положение фазового перехода, что может быть сделано также при помощи термодинамического интегрирования, но и определить уравнения состояния
и другие свойства системы около перехода, которое термодинамическим интегрированием получены быть не могут.
Работа выполнена с использованием ресурсов суперкомпьютерного комплекса МГУ имени М. В. Ломоносова [10] при частичной финансовой поддержке Минобрнауки РФ (грант № 8376).
Список литературы
1. McMahon J.M., Morales M.A., Pierleoni C., Ceperley D.M.
// Rev. Mod. Phys. 2012. 84. P. 1607.
2. Mao H, Hemley R.J. // Rev. Mod. Phys. 1994. 66. P. 671.
3. Eremets M.I., Troyan I.A. // Nature materials. 2011. 10. P. 927.
4. Thomas L.H. // Proc. Cambridge Phil. Soc. 1927. 23 (5). P. 542.
5. Ashcroft N.W., Mermin N.D. Solid State Physics. Toronto, 1976.
6. Militzer B., Graham R.L. // J. of Phys. and Chem. of Solids. 2006. 67. P. 2136.
7. Feynman R.P., Hibbs A. // Quantum Mechanics and Path Integrals. N.Y., 1965.
8. Ceperley D.M. // Rev. Mod. Phys. 1995. 67. P. 279.
9. Militzer B. // Phys. Rev. B. 2009. 79. P. 155105.
10. Воеводин Вл.В., Жуматий С.А., Соболев С.И. и др. Практика суперкомпьютера «Ломоносов» // Открытые системы. М., 2012.
Monte Carlo modeling of metallic hydrogen: the phase transition and the equation of state
A. A. Novoselov1,a, O. V. Pavlovsky2, M. V. Ulybyshev2
1 Department of Quantum Theory and High-Energy Physics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.
2 Institute for Theoretical Problems of Microworld, M. V. Lomonosov Moscow State University, Moscow 119991, Russia
E-mail: a novoselov@goa.bog.msu.ru.
We conducted numerical modeling of atomic (metallic) hydrogen using the PIMC (path integral Monte Carlo) method. The temperature and density range in which the electron (proton) behavior is governed by quantum (classical) statistics was studied. The equations of state in the form of dependences of the internal energy and pressure on temperature and density were obtained in that region. These dependences allow one to reveal and study the phase transition between crystal and liquid phases.
Keywords: Monte Carlo method, metallic hydrogen, equation of state. PACS: 62.50.-p, 61.20.Ja, 64.30.-t, 67.80.F-.
Received 2 October 2013.
English version: Moscow University Physics Bulletin 1(2014). Сведения об авторах
1. Новоселов Александр Андреевич — аспирант; e-mail: novoselov@goa.bog.msu.ru.
2. Павловский Олег Владимирович — канд. физ.-мат. наук, ст. науч. сотрудник; тел.: (495) 939-26-96, e-mail: ovp@goa.bog.msu.ru.
3. Улыбышев Максим Владимирович — канд. физ.-мат. наук, ст. науч. сотрудник; тел.: (495) 939-26-96.
