ЖЭТФ, 2019, том 156, вып. 3 (9), стр. 545-556

СВОЙСТВА КРИСТАЛЛОВ И ЖИДКОСТИ ЮКАВЫ

В УСЛОВИЯХ ФАЗОВОГО РАВНОВЕСИЯ

В. В. Решетняк^*, А. В. Филиппов

ГНЦ РФ Троицкий институт инновационных и термоядерных исследований

108840, Троицк, Москва, Россия

Поступила в редакцию 28 февраля 2019 г.,

после переработки 28 февраля 2019 г.

Принята к публикации 12 марта 2019 г.

Изучены равновесные свойства системы частиц, взаимодействие которых описывается потенциалом Юка-

вы. Рассмотрены жидкость и кристаллы с ОЦК- и ГЦК-решетками вблизи линии плавления. Вычисле-

ны величины скрытой теплоты плавления кристаллов и изменения плотности при фазовом переходе.

Проанализированы температурные зависимости термодинамических величин вблизи линии плавления,

рассчитаны значения коэффициента теплового расширения и теплоемкости при постоянном давлении.

Показано, что состояние системы Юкавы не может быть однозначно определено парой безразмерных

величин: температурой и давлением, использующим в качестве масштабов параметры потенциала вза-

имодействия (произведение зарядов и длину экранирования). В частности, положение линии плавления

неоднозначно определено на двухпараметрической фазовой диаграмме, построенной в безразмерных ве-

личинах.

DOI: 10.1134/S0044451019090189

ность их изучения на наиболее подробном кинети-

ческом уровне, отслеживая траектории движения

1. ВВЕДЕНИЕ

«атомов», в роли которых выступают частицы КДФ

Плазма, содержащая погруженные в нее дис-

[5-14].

персные частицы конденсированной фазы (КДФ),

Полное теоретическое описание свойств пыле-

широко распространена в природе и технике и пред-

вой плазмы предполагает учет многокомпонентно-

ставляет значительный интерес как для фундамен-

сти, неравновесности, неидеальности, открытости

тальной науки, так и для ряда приложений [1-4].

системы и потому связано со значительными труд-

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

ностями. Однако задачу можно упростить, исполь-

КДФ приобретают электрический заряд, существен-

зуя приближение однокомпонентной системы час-

но превышающий по своей величине заряды элек-

тиц, взаимодействующих посредством псевдопотен-

тронов и ионов плазмы. Величина энергии электро-

циалов. Модель системы частиц Юкавы (дебаевская

статического взаимодействия макрочастиц при этом

модель) является одним из наиболее распростра-

зачастую достигает либо заметно превышает значе-

ненных приближений для теоретического изучения

ние энергии их теплового движения. Таким обра-

комплексной плазмы. Вычисление электростатиче-

зом, использование модели идеальной плазмы при

ского потенциала в плазме осуществляется в рам-

вычислении равновесных параметров системы ока-

ках теории Дебая - Хюккеля. Предполагается, что

зывается некорректным.

фоновая плазма равновесна, неподвижна, идеальна

Ансамбль сильно взаимодействующих частиц

и слабо взаимодействует с частицами КДФ. Разме-

КДФ, удерживаемых полем конфайнмента в огра-

ром и формой частиц пренебрегают. Тогда парное

ниченном объеме, может проявлять свойства кон-

взаимодействие в однокомпонентной системе частиц

денсированного вещества и испытывать фазовые пе-

КДФ описывается потенциалом

реходы. Это позволяет моделировать в лаборатор-

(

)

q_iq_j

rij

ных экспериментах различные процессы, протека-

u(r_ij ) =

exp

≡

ющие в конденсированных средах, и дает возмож-

4πε0 rij

λ_D

u₀

≡

exp(-k_Dr) .

(1)

^* E-mail: viktor.reshetnyak84@gmail.com

545

ЖЭТФ, вып. 3 (9)

В. В. Решетняк, А. В. Филиппов

ЖЭТФ, том 156, вып. 3 (9), 2019

Здесь r_ij — расстояние между частицами с индек-

статочно изученными остаются свойства жидкости

сами i и j, которое далее будет обозначаться просто

и кристалла на линии плавления. Сама же линия

как r, q_i,j — заряды частиц, которые часто принима-

плавления определена приближенно, о чем, в част-

ются равными между собой (q_i = q_j = q), u₀ — пара-

ности, свидетельствуют заметные расхождения ре-

метр, определяющий величину кулоновского взаи-

зультатов работ разных авторов, достигающие 10 %

модействия частиц, λ_D и k_D — длина дебаевского

[29]. Отчасти отмеченные расхождения можно объ-

экранирования и соответствующая ей обратная ве-

яснить тем фактом, что в ряде работ при исследо-

личина. Для изучения свойств систем Юкавы широ-

вании фазовой диаграммы из рассмотрения был ис-

ко используются приведенные параметры потенци-

ключен скачок плотности, возникающий при фазо-

ала (1): параметр неидеальности Γ = βu₀/δ и струк-

вом переходе. Например, в работах [24,26-30,36] для

турный параметр κ = δk_D. Здесь и далее β-1 =

определения фазового равновесия использовалось

= k_BT, k_B — постоянная Больцмана, T — темпера-

условие равенства свободной энергии Гельмгольца в

тура, n — концентрация частиц, δ = n-1/3 — среднее

различных фазах. Использование такого приближе-

межчастичное расстояние.

ния может быть оправдано малостью скачка плот-

Несмотря на простоту, дебаевская модель хоро-

ности на границах раздела фаз в системах Юкавы

шо зарекомендовала себя для описания структу-

[23, 25].

ры и свойств плазменно-пылевых систем. В рабо-

Настоящая работа посвящена исследованию рав-

тах [15-18] потенциал парного взаимодействия был

новесных свойств системы частиц Юкавы в ОЦК-,

восстановлен с использованием экспериментально

ГЦК- и жидкой фазах вблизи линии плавления. Рас-

измеренных траекторий частиц. Результаты свиде-

четы выполнены с использованием метода молеку-

тельствуют о том, что потенциал парного взаимо-

лярной динамики (МД). Сопоставление результатов

действия хорошо описывается выражением (1), хотя

расчетов с данными других авторов показало удо-

параметры u₀ и λ_D могут отличаться от дебаевских

влетворительное согласие с ними. Кроме того, для

значений. Аналитические и численные модели, учи-

жидкой фазы был применен метод интегральных

тывающие многокомпонентность системы, неравно-

уравнений (ИУ) теории жидкости [40], использую-

весность и неидеальность фоновой плазмы, также

щий приближения, отличные от стандартных допу-

указывают на экспоненциальную зависимость по-

щений МД. Результаты вычислений методами МД и

тенциальной энергии от расстояния между парой

ИУ совпадают с достаточно высокой точностью.

макроионов.

В отличие от работ [22, 24,26-30, 36], где изуче-

Для теоретического исследования фазовых пре-

ние фазовой диаграммы системы частиц Юкавы вы-

вращений традиционно используются два подхо-

полнено в переменных κ и Γ, в настоящей работе

да. Первый предполагает изучение равновесных

для исследования выбраны безразмерные давление

свойств различных фаз в разных точках фазовой

и температура p и

T, связанные с соответствующи-

диаграммы. Второй подход связан с изучением ре-

ми размерными величинами p и T формулами

лаксации неравновесной системы частиц при фа-

λ3D

k_BTλ_D

зовом переходе (например, кристаллизация пере-

p=p

T =

(2)

E₀

u₀

охлажденной жидкости). Фазовые переходы в трех-

мерных системах частиц Юкавы исследовались ра-

где E₀ — масштаб энергии, величина которого не

нее с использованием как первого, так и второго

изменялась в ходе расчетов и задавалась равной

подхода [12, 13, 19-39]. Результаты свидетельству-

k_BT_c, Tc = 300 K — постоянное значение темпе-

ют о том, что при определенных термодинамиче-

ратуры, примерно соответствующее условиям экс-

ских условиях для системы Юкавы возможен фа-

периментов. Также неизменной в настоящей работе

зовый переход первого рода между жидкой фазой

задавалась длина дебаевского экранирования λ_D =

и одной из двух кристаллических фаз: объемноцен-

= 5·10-3 см. Здесь и везде далее в статье безразмер-

трированной кубической (ОЦК) или гранецентриро-

ные величины обозначаются символом с тильдой, а

ванной кубической (ГЦК). Переход между ОЦК- и

отсутствие тильды указывает на наличие размерно-

ГЦК-фазами также изучался, в то время как о су-

сти.

ществовании стабильных трехмерных кристалличе-

Такой выбор переменных обусловлен прежде все-

ских фаз с другими структурами неизвестно. Заме-

го тем фактом, что на линии фазового равновесия

тим, однако, что, несмотря на обилие публикаций

плотность испытывает скачок [23], следовательно,

по данному вопросу, проблему фазовых переходов

параметры κ и Γ определены неоднозначно. Анало-

в системах Юкавы нельзя считать закрытой. Недо-

гичный нашему выбор масштабов при исследовании

546

ЖЭТФ, том 156, вып. 3 (9), 2019

Свойства кристаллов и жидкости Юкавы.. .

фазовой диаграммы системы Юкавы был использо-

нии (3) предполагается, что на границе раздела фаз

ван в работах [23,25].

выполняется равенство как температур, так и па-

Заметим, однако, что использование безразмер-

раметров потенциала (1). Использование уравнений

ных параметров, введенных в настоящей работе, без

(3) предполагает допущение о малости лапласового

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

давления на границе раздела фаз, вклад которого

бы привести к неоднозначности результатов. Так, к

не учитывается.

примеру, одному и тому же значению безразмерной

Анализируя температурную зависимость энтро-

температуры

T могут соответствовать различные T ,

пии и объема при постоянном давлении, можно

u₀ и λ_D. В отличие от обыкновенных конденсиро-

вычислить теплоемкость и коэффициент теплового

ванных веществ, где параметры потенциала парно-

расширения по формулам

го взаимодействия атомов неизменны, в ходе экс-

периментов с пылевой плазмой изменение величи-

⁽∂s⁾

⁽∂ñ⁾

cp =

α=-

(4)

ны

T может быть достигнуто управлением как тем-

∂T

T p,u0,λD

p,u0,λD

пературой частиц, так и параметрами потенциала

взаимодействия, причем эквивалентность перечис-

Индекс p, u₀, λ_D указывает на то, что дифференци-

ленных способов неочевидна. В настоящей работе

рование проводится при постоянных давлении и па-

во всех расчетах дебаевская длина предполагалась

раметрах потенциала.

неизменной, температура же и константа связи u₀

Кроме величин c_p и α в настоящей работе так-

варьировались.

же были исследованы параметры ĉ_p и α, определен-

Сначала для заданного безразмерного давления

ные аналогично (4) с тем лишь отличием, что при

выбиралось приближенное значение

T, причем зна-

расчете безразмерных температурных зависимостей

чение T задавалось равным T_c, а значение констан-

объема и энтропии постоянным считается значение

ты связи вычислялось таким образом, чтобы значе-

размерной температуры, а значение параметра u₀

ние эффективного параметра неидеальности,

варьируется:

(

)

Γ^∗ = Γ

1 + κ + 0.5κ²

e^-κ,

⁽∂s⁾

⁽∂ñ⁾

примерно соответствовало критерию плавления

ĉp =

α=-

(5)

∂T

[32, 34]. Согласно этому критерию плавление крис-

T p,T,λD

p,T,λD

таллов Юкавы осуществляется при Γ^∗

≈ 106.

Параметр κ определяется безразмерной плотностью

По аналогии с c_p и α далее в настоящей работе

κ = ñ-1/3, ñ = nλ3D, которая, в свою очередь,

будем называть эти параметры эффективной тепло-

вычислялась в программе LAMMPS [41, 42] для

емкостью и эффективным коэффициентом теплово-

заданного давления с использованием методов

го расширения. Еще раз подчеркнем, что поскольку

расчета NP T -ансамбля [43].

эффективные теплоемкость и коэффициент тепло-

Затем значение безразмерной температуры, со-

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

ответствующей фазовому равновесию при заданном

чении T и переменном u₀, их физический смысл от-

давлении, уточнялось путем исследования линей-

личен от смысла параметров c_p и α. Например, ве-

ного отклика безразмерных термодинамических па-

личина ĉ_p характеризует работу по изменению эн-

раметров на малые изменения

T. Рассматривались

тропии плазменно-пылевой системы при изменении

температурные зависимости химического потенциа-

заряда (а не температуры) макроионов при постоян-

ла μ = μ/E₀, плотности ñ, энтропии одной частицы

ном давлении. Постоянство давления при этом обес-

ансамбля s = s/k_B. При этом был проанализиро-

печивается за счет изменения плотности системы,

ван вопрос об эквивалентности двух разных спосо-

которое характеризует величина α. Таким образом,

бов изменения безразмерной температуры, для чего

численные значения теплоемкости и коэффициента

поочередно варьировались значения размерной тем-

теплового расширения также могут не совпадать с

пературы T и константы связи u₀. Положение точ-

соответствующими значениями эффективных пара-

ки фазового равновесия вычислялось путем реше-

метров.

ния уравнений

Кроме того, заметим, что в настоящей работе

расчеты выполнялись для однокомпонентной сис-

p₁ = p₂,

T₁ =

T₂,

μ₁ = μ₂,

(3)

темы, без учета самосогласованного характера из-

где индексы 1 и 2 обозначают различные фазы. При

менения параметров, имеющего место в реальной

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

плазменно-пылевой системе.

547

11*

В. В. Решетняк, А. В. Филиппов

ЖЭТФ, том 156, вып. 3 (9), 2019

2. ОПИСАНИЕ МОДЕЛИ

ме частиц путем введения параметра связи λ: u =

= u(λ). Параметр связи изменяют в пределах от 0

2.1. Метод молекулярной динамики

до 1 таким образом, чтобы при этом осуществлял-

Расчет термодинамических потенциалов осу-

ся переход из опорного состояния 1 в конечное со-

ществлялся с использованием метода МД, реализо-

стояние 2. Конечное состояние — исследуемый ан-

ванного в программе LAMMPS [41,42]. Траектории

самбль частиц, а в качестве опорного обычно вы-

вычислялись в кубической периодической ячейке

бирают состояние системы, свободная энергия ко-

для системы из 16000 атомов для ОЦК-кристалла

торого может быть вычислена аналитически. Тогда

и жидкости и 16384 атомов для ГЦК-кристалла.

избыточная свободная энергия может быть вычис-

Взаимодействие частиц определялось потенциа-

лена по изменению потенциальной энергии системы

лом

(1). Для вычисления радиуса обрезки r_cut

U согласно выражению

потенциала Юкавы использовалась оценка отно-

∫

сительного вклада хвостовой части интеграла для

⁷∂U⁸

F_ex = F₂ - F₁ =

dλ.

(6)

потенциальной энергии, детально обсуждавшаяся в

∂λ_λ

работе [44]. При этом значения r_cut в зависимости

от κ принимали значения от 13λ_D и выше. Размер

При расчете свободной энергии жидкости Юкавы в

вычислительной области во всех случаях превышал

качестве опорного состояния был выбран идеальный

удвоенное значение r_cut.

газ, свободная энергия которого известна [47]. Дета-

Численное интегрирование уравнений движения

ли вычислений свободной энергии жидкости описа-

выполнялось с использованием алгоритма Верле

ны в работе [44].

[45]. Временной шаг выбирался из условия Δt ≪

При расчете свободной энергии кристалла в ка-

≪ t₀ = δ/v_T, где t₀ — характерное время, v_T

честве опорного состояния использовался кристалл

= (k_B T/M)1/2 — средняя скорость теплового дви-

Эйнштейна [46], который представляет собой со-

жения частиц, а M = 4.19 · 10-12 г — масса частиц.

вокупность невзаимодействующих частиц, связан-

Значение Δt задавалось равным 4 · 10-4 c.

ных с узлами кристаллической решетки гармониче-

При расчете свойств кристаллов в начальный мо-

ским потенциалом. Жесткость связей a определя-

мент времени частицы располагались в узлах кри-

лась с использованием предварительно рассчитан-

сталлической решетки, а для жидкостей начальные

ного среднего квадратического смещения частиц в

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

состоянии 2 (кристалле Юкавы):

Объем системы подбирался из условия соответствия

параметров системы критерию плавления Γ^∗ = 106

(7)

[32, 34]. Затем координаты и импульсы частиц ре-

2β 〈(Δr)²〉

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

Свободная энергия системы не взаимодействующих

и температуре с использованием реализованного в

друг с другом осцилляторов известна [46, 48]:

программе LAMMPS алгоритма для NP T -ансамбля

[

]

[43]. На релаксацию отводилось 10000 расчетных

(

)3/2

шагов, что примерно на два порядка больше харак-

Fharm = -

ln Λ-3δ³

(8)

aδ²β

терного времени t₀. После релаксации системы вы-

числялся средний равновесный объем системы при

(

)1/2

Здесь Λ =

βh²/2πM

, h — постоянная План-

заданных условиях. Для определения среднего объ-

ка. Запишем выражение для приходящейся на один

ема выполнялось 10 статистически независимых вы-

атом свободной энергии эйнштейновского кристалла

числений, в каждом из которых усреднение прово-

f_harm, выраженной через соответствующее значение

дилось на временном интервале 1000 шагов. Далее,

для идеального газа f_ig:

при заданном объеме, вычисленном для выбранно-

го давления согласно описанной выше процедуре,

методом термодинамического интегрирования (ТИ)

f_harm = f_ig +

(9)

2β

aδ²β

была рассчитана свободная энергия Гельмгольца си-

стемы частиц. При этом использовался термостат

Для исследования фазовой диаграммы вблизи

Ланжевена для NV T-ансамбля.

линии плавления необходим учет ангармонизма ко-

Метод ТИ подробно описан в литературе (см.,

лебаний атомов в решетке. Вклад в свободную энер-

например, [46]). Идея метода состоит в непрерыв-

гию, связанный с ангармонизмом колебаний, вычис-

ном изменении потенциала взаимодействия в систе-

лялся методом ТИ. Параметр связи λ при этом вво-

548

ЖЭТФ, том 156, вып. 3 (9), 2019

Свойства кристаллов и жидкости Юкавы.. .

дился таким образом, что потенциальная энергия

Статический структурный фактор связан с пар-

системы частиц u определялась выражением

ной корреляционной функцией фурье-преобразова-

(

)

нием:

∫

u = λ^mu_p + (1 - λ^m)

u_b + u0p

(10)

S(k) = 1 + n h(r)eik·rdr.

(14)

Здесь u_p — энергия парного взаимодействия частиц,

Учитывая изотропность системы и сферическую

u_b — потенциальная энергия гармонического осцил-

симметрию парного потенциала, перепишем выра-

лятора, u0p — энергия парного взаимодействия час-

жение (14) в виде

тиц в кристалле, все атомы которого зафиксирова-

ны в узлах решетки. Показатель степени m — произ-

∫∞

4πn

вольное положительное число, которое в настоящей

S(k) = 1 +

h(r) sin(kr) r dr.

(15)

работе было выбрано равным 2.

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

Уравнение ОЦ решалось методом итераций со-

ратуры и объема были рассчитаны внутренняя энер-

гласно алгоритму, описанному в работах [50, 51].

гия e_int и свободная энергия Гельмгольца, энтропия

Для ускорения сходимости использовалась процеду-

и химический потенциал вычислялись по формулам

ра, предложенная в работе [52]. Контроль самосо-

s_ex = s - s_ig = (e_int - f_ex)/T,

гласованности решения осуществлялся путем срав-

(11)

нения вычисленных разными способами значений

μ_ex = μ - μ_ig = f_ex + pv - k_BT.

S(k) в точке k = 0. Первый способ предполагал ли-

Для каждого заданного значения давления описан-

нейную экстраполяцию данных по двум ближайшим

ная процедура расчета термодинамических пара-

к нулю точкам, а второй — использование уравнения

метров повторялась для пяти значений безразмер-

для изотермической сжимаемости K_T [47, 53]:

ной температуры

T, находящихся вблизи исходно-

lim

S(k) = nk_BT K_T ,

го значения. Варьирование безразмерной темпера-

k→0

туры осуществлялось двумя способами: при посто-

(16)

^[(∂P)]-1

янном значении T путем изменения u₀ и наоборот,

K_T =

∂n_T

после чего результаты сравнивались. Затем вычис-

лялись параметры линейной аппроксимации темпе-

При этом для определения производной давления P

ратурной зависимости термодинамических свойств

по концентрации пылевых частиц n уравнение ОЦ

системы вблизи исходной точки. Температура и хи-

решалось при двух близких значениях n.

мический потенциал в точке фазового равновесия

Давление вычислялось с использованием вири-

вычислялись с использованием полученных линей-

ального уравнения состояния [40]:

ных функций μ

T) для разных фаз путем решения

∫

n²

∂u(r)

системы уравнений (3). Значения остальных термо-

P = nk_BT -

g(r) dr,

(17)

∂r

динамических величин в этой точке определялись

по предварительно рассчитанным в линейном при-

которое для жидкости Юкавы принимает вид

ближении температурным зависимостям.

∞

∫

2πu₀n²

P = nk_BT +

(1 + k_Dr)e-kDrg(r) r dr. (18)

2.2. Метод интегральных уравнений

Использование метода ИУ связано с решением

Свободная энергия вычислялась с использовани-

уравнения Орнштейна - Цернике (ОЦ), которое име-

ем метода термодинамического интегрирования, со-

ет вид [40]

∫

гласно описанному в [44] алгоритму.

h(r) = c(r) + n h(r₁)c(|r - r₁|) dr₁,

(12)

3. РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

где h(r) — парная корреляционная функция, c(r) —

прямая корреляционная функция. Для замыкания

Значения безразмерных термодинамических па-

уравнения ОЦ использовалось гиперцепное (ГЦ)

раметров систем Юкавы в ОЦК-, ГЦК- и жидкой

приближение [49]:

фазах на линии плавления-кристаллизации пред-

ставлены в табл. 1. Для жидкости отдельно приведе-

c(r) = exp[-βu(r) + γ(r)] - γ(r) - 1,

(13)

ны параметры в точках, соответствующих кристал-

где γ(r) = h(r) - c(r).

лизации с образованием ГЦК- и ОЦК-фазы.

549

В. В. Решетняк, А. В. Филиппов

ЖЭТФ, том 156, вып. 3 (9), 2019

Таблица 1. Параметры системы Юкавы на линии плавления ОЦК- и ГЦК-кристаллов

ОЦК

T = const = T0, варьируется u0

u0 = const, варьируется T

T · 10⁴

μex

sex

ñ · 10²

Γ^∗

T · 10⁴

μex

sex

ñ · 10²

Γ^∗

0.025

0.034

30.45

-4.53

0.117

9.49

130.1

0.035

29.98

-4.63

0.119

9.44

132.7

0.05

0.155

36.06

-4.51

0.200

7.94

117.4

0.145

37.00

-4.52

0.195

8.00

118.1

0.10

0.552

42.86

-4.46

0.343

6.63

106.7

0.484

45.15

-4.50

0.325

6.75

109.3

0.20

1.365

53.88

-4.50

0.558

5.64

104.2

1.338

54.68

-4.49

0.549

5.67

103.5

0.40

2.919

67.97

-4.52

0.906

4.80

102.0

2.749

70.57

-4.52

0.870

4.86

102.3

0.80

5.766

84.95

-4.48

1.488

4.07

97.5

5.146

90.69

-4.53

1.389

4.16

100.7

1.60

10.03

106.8

-4.44

2.419

3.46

94.8

7.951

123.3

-4.60

2.090

3.63

103.0

3.20

13.66

146.8

-5.29

3.624

3.02

101.3

12.83

154.4

-4.63

3.435

3.08

102.9

6.40

22.60

177.3

-4.84

6.086

2.54

92.8

17.26

214.4

-4.68

5.058

2.70

105.6

12.8

24.93

260.7

-4.72

8.563

2.27

106.8

26.36

269.7

-4.56

8.261

2.30

98.6

ГЦК

0.025

0.037

29.93

-4.61

0.120

9.41

128.6

0.036

30.08

-4.65

0.119

9.44

129.0

0.05

0.161

35.73

-4.60

0.203

7.90

117.1

0.146

37.02

-4.62

0.196

7.99

117.3

0.10

0.534

43.29

-4.58

0.341

6.64

109.1

0.480

45.13

-4.61

0.326

6.74

110.2

0.20

1.292

54.67

-4.62

0.549

5.67

107.2

1.310

54.63

-4.61

0.550

5.67

105.7

0.40

2.785

68.98

-4.63

0.892

4.82

104.6

2.658

70.47

-4.65

0.871

4.86

105.8

0.80

5.051

89.00

-4.67

1.417

4.13

105.1

4.927

90.56

-4.67

1.390

4.16

105.2

1.60

8.333

115.7

-4.70

2.238

3.55

105.5

7.641

123.2

-4.73

2.092

3.63

107.2

3.20

13.26

148.6

-5.37

3.580

3.03

103.3

12.29

154.3

-4.74

3.437

3.08

107.5

6.40

20.01

189.6

-5.00

5.741

2.60

100.3

16.31

214.2

-4.82

5.060

2.70

111.7

12.8

26.23

255.0

-4.67

8.763

2.25

103.1

23.95

269.4

-4.76

8.267

2.30

108.6

Жидкость-ОЦК

0.025

0.034

30.45

-3.82

0.114

9.57

120.6

0.035

29.98

-3.85

0.116

9.51

123.3

0.05

0.155

36.06

-3.80

0.197

7.98

113.2

0.145

37.00

-3.81

0.192

8.04

113.8

0.10

0.552

42.86

-3.75

0.338

6.66

103.7

0.484

45.15

-3.79

0.321

6.78

106.6

0.20

1.365

53.88

-3.79

0.552

5.66

102.4

1.338

54.68

-3.77

0.543

5.69

101.6

0.40

2.919

67.97

-3.80

0.898

4.81

100.8

2.749

70.57

-3.80

0.863

4.88

101.1

0.80

5.766

84.95

-3.76

1.479

4.07

96.8

5.146

90.69

-3.82

1.381

4.17

100.0

1.60

10.03

106.8

-3.70

2.408

3.46

94.4

7.951

123.3

-3.89

2.081

3.64

102.6

3.20

13.66

146.8

-3.86

3.614

3.02

101.1

12.83

154.4

-3.90

3.425

3.08

102.7

6.40

22.60

177.3

-3.62

6.071

2.54

92.6

17.26

214.4

-3.97

5.047

2.71

105.4

12.8

24.93

260.7

-3.99

8.552

2.27

106.8

26.36

269.7

-3.83

8.247

2.30

98.5

Жидкость-ГЦК

0.025

0.037

29.93

-3.77

0.116

9.52

116.6

0.036

30.08

-3.79

0.116

9.52

119.8

0.05

0.161

35.73

-3.77

0.199

7.95

111.7

0.146

37.02

-3.79

0.191

8.06

111.6

0.10

0.534

43.29

-3.78

0.335

6.68

105.3

0.480

45.13

-3.80

0.321

6.78

107.5

0.20

1.292

54.67

-3.83

0.543

5.69

105.2

1.310

54.63

-3.81

0.544

5.69

104.1

0.40

2.785

68.98

-3.85

0.884

4.84

103.3

2.658

70.47

-3.87

0.864

4.87

104.8

0.80

5.051

89.00

-3.90

1.408

4.14

104.3

4.927

90.56

-3.91

1.382

4.17

104.5

1.60

8.333

115.7

-3.93

2.229

3.55

105.1

7.641

123.2

-3.96

2.082

3.64

106.8

3.20

13.26

148.6

-3.90

3.569

3.04

103.0

12.29

154.3

-3.99

3.427

3.08

107.3

6.40

20.01

189.6

-3.82

5.727

2.59

100.1

16.31

214.2

-4.07

5.050

2.71

111.6

12.8

26.23

255.0

-3.92

8.750

2.25

103.0

23.95

269.4

-4.02

8.253

2.30

108.5

550

ЖЭТФ, том 156, вып. 3 (9), 2019

Свойства кристаллов и жидкости Юкавы.. .

Таблица 2. Параметры жидкости Юкавы на линии плавления ОЦК-фазы, определенные методами молекулярной

динамики (МД) и интегральных уравнений жидкости (ИУ)

МД

ИУ

ñ · 10²

T · 10⁴

sex

8.55

24.93

12.8

107.9

-3.99

12.9

109.2

-3.89

6.07

22.60

6.40

71.2

-3.62

6.59

74.7

-3.68

3.61

13.66

3.20

55.4

-3.86

3.24

56.6

-3.87

2.41

10.03

1.60

38.2

-3.70

1.65

39.9

-3.75

1.48

5.766

0.80

28.1

-3.76

0.820

29.2

-3.80

0.898

2.919

0.40

20.6

-3.80

0.413

21.6

-3.87

0.552

1.365

0.20

14.8

-3.79

0.209

15.7

-3.87

0.338

0.552

0.10

10.6

-3.75

0.107

11.5

-3.85

0.197

0.155

0.05

7.84

-3.80

0.054

8.64

-3.93

0.114

0.034

0.025

5.76

-3.82

0.028

6.54

-4.00

Вычисление значения параметра неидеальности

затем сравнивались с результатами, полученными

Γ, структурного параметра κ, а также эффективно-

методом МД. Результаты сравнения представлены

го параметра неидеальности Γ^∗ осуществлялось по

в табл. 2.

рассчитанным значениям безразмерных температу-

Как видно из табл. 2, данные, полученные с ис-

ры и плотности согласно формулам κ = ñ-1/3, Γ =

пользованием методов МД и ИУ, совпадают доста-

= (

T)-1. Значения параметров Γ, κ и Γ^∗ для ОЦК-

точно хорошо, при этом наилучшее согласие резуль-

и ГЦК-структур в точках с одинаковыми значения-

татов имеет место при больших значениях

T и ñ:

ми давлений, температур и параметров потенциала

максимальным значениям параметров соответству-

различаются незначительно: отличия соответствую-

ет разница вычисленных величин, не превышающая

щих значений Γ и κ не превышают 0.2 %, а значе-

1.5 %. При уменьшении

T и ñ разница становит-

ний Γ^∗ — 1 %. В то же время для ОЦК-кристалла

ся более существенной и достигает 13 % для мини-

и жидкости различия более существенны: соответ-

мальных значений параметров в выбранном диапа-

ствующие параметры Γ и κ отличаются примерно

зоне. Объяснением данного результата может слу-

на 1%, а параметры Γ^∗ — примерно на 8 %.

жить ограниченная применимость метода ИУ для

Неопределенность, связанная с неоднозначнос-

сильно неидеальных систем. Степень неидеальности

тью выбора способа изменения безразмерной темпе-

системы определяется параметром Γ, который мо-

ратуры, оказывается заметно больше: для одного и

нотонно убывает с ростом

T и ñ. Так, для парамет-

того же давления ОЦК-кристалла значения безраз-

ров

T = 2.493 · 10-3 и ñ = 8.55 · 10-2, при которых

мерной температуры в точке плавления могут от-

совпадение результатов наилучшее, значение пара-

личаться более чем на 30 %, а вычисленная разница

метра неидеальности Γ ≈ 176, в то время как наи-

значений κ, Γ и Γ^∗ достигает 6 %, 19 % и 12 % соот-

большие различия наблюдаются при

T = 3.37·10-6,

ветственно.

ñ = 1.14 · 10-3 и Γ ≈ 3.094 · 10⁴.

Для тестирования результатов вычислений ис-

Результаты вычислений также находятся в

пользовались результаты, полученные в настоящей

неплохом согласии с данными других авторов.

работе методом ИУ, а также данные, известные из

Среднее значение Γ^∗ на линии кристаллизации

работ других авторов. Уравнение ОЦ решалось для

жидкости

〈Γ^∗〉_T

≈ 103.2, 〈Γ∗〉u0

≈ 103.7 (см.

жидкости Юкавы при температуре T = 300 K, пара-

табл. 3), в то время как по данным [32], где были

метры которой (плотность и безразмерная темпера-

обобщены и проанализированы результаты работы

тура) задавались с использованием представленных

[29], 〈Γ^∗〉 = 106.9.

в табл. 1 данных в условиях фазового равновесия

Сравнение результатов с данными [29] также

с ОЦК-кристаллом. Вычислялись давление, потен-

указывает на разумное соответствие положений

циальная энергия и избыточная энтропия, которые

кривых плавления. Для сравнения результатов нас-

551

В. В. Решетняк, А. В. Филиппов

ЖЭТФ, том 156, вып. 3 (9), 2019

Таблица 3. Линия плавления в переменных κ, Γ и Γ^∗

Γ^∗

u₀ = const

T = const

Данные [29]

u₀ = const

T = const

Данные [32]

2.26

161.3

169.2

166.8

97.8

103.0

101.3

3.22

273.3

261.0

273.0

102.7

98.1

102.3

4.19

472.8

460.6

470.8

100.0

97.4

99.1

4.84

728.1

726.7

735.1

101.1

100.9

102.3

5.80

1428

1434

1475

102.1

102.6

105.5

6.45

2353

2316

2380

105.1

103.4

106.6

7.42

5067

5033

5341

109.1

108.4

115.1

8.06

8511

8537

9343

111.7

112.0

122.6

тоящей работы с данными [29] была использована

несколько выше: L_fcc = 0.85. Этот результат непло-

линейная интерполяция зависимостей Γ(κ), рассчи-

хо согласуется с данными работы [23]. Положитель-

танных для жидкости Юкавы на линии кристалли-

ное значение теплоты плавления кристаллов гово-

зации. При этом считалось, что если безразмерная

рит о том, что плавление осуществляется при увели-

температура плавления кристалла с ОЦК-решеткой

чении температуры. Заметим, что разность энтро-

в данной точке выше, чем для ГЦК-кристалла, то

пий на линиях плавления ОЦК- и ГЦК-кристаллов

кристаллизация жидкости происходит с образова-

Δs мала и соизмерима со значением дисперсии. Та-

нием ОЦК-фазы, а в обратном случае — с образо-

ким образом, исследование фазовых переходов меж-

ванием ГЦК-фазы. Результаты сопоставления сви-

ду ОЦК- и ГЦК-фазами требует повышения точ-

детельствуют о совпадении с результатами [29] для

ности вычислений. В частности, может понадобит-

различных точек фазовой диаграммы в пределах

ся увеличение выборки для снижения статистиче-

10 %. Имеющиеся расхождения могут быть связаны

ской ошибки вычислений, увеличение радиуса об-

как с неоднозначностью определения линии плавле-

резки потенциала, уточнение положения точек фа-

ния в безразмерных величинах, так и с использован-

зового равновесия.

ной в настоящей работе линейной аппроксимацией.

Другим важным параметром, входящим в урав-

Анализ данных табл. 1 свидетельствует о сла-

нение Клапейрона - Клаузиуса, является разность

бой зависимости энтропии от давления. Кроме того,

удельных объемов для разных фаз — величина, ха-

неоднозначность выбора способа изменения безраз-

рактеризующая скачок плотности при фазовом пе-

мерной температуры мало сказывается на значении

реходе. Анализ табл. 1 свидетельствует о том, что

энтропии (средние значения различаются примерно

плотности ГЦК- и ОЦК-кристаллов на линии плав-

на 2 %). Среднее значение энтропии ОЦК-кристал-

ления мало различаются. При плавлении плотность

ла составляет 〈s^bccex〉 ≈ -4.57. Для ГЦК-кристалла

кристаллов изменяется скачком, причем величина

〈s^fccex〉 ≈ -4.69. Энтропия жидкости выше, ее значе-

этого изменения зависит от давления и темпера-

ния изменяются от -3.99 до -3.70, а среднее 〈s^liqex 〉 ≈

туры. Зависимость отношения плотностей кристал-

≈ -3.84.

лических и жидкой фаз от температуры на линии

плавления представлена на рисунке.

Разность энтропий на линии фазового равнове-

сия позволяет вычислить скрытую теплоту перехо-

Из рисунка видно, что температурная зависи-

да между фазами 1 и 2: L₂₁ = T (s₂ - s₁). Скрытая

мость отношения плотностей жидкой и кристалли-

теплота входит в качестве параметра в уравнение

ческих фаз неплохо аппроксимируется логарифми-

Клапейрона - Клаузиуса, решение которого широко

ческой функцией. С ростом температуры (соответ-

используется для описания фазовых переходов пер-

ственно, на линии плавления растет и давление) во

вого рода. Вычисленная в единицах E₀ с использо-

всем рассмотренном диапазоне параметров величи-

ванием средних значений энтропии скрытая тепло-

на ñ_cr/ñ_liq медленно убывает, оставаясь при этом

та плавления ОЦК-кристалла L_bcc = 0.73, в то вре-

больше 1. Индексом «cr» здесь обозначен кристалл

мя как скрытая теплота плавления ГЦК-кристалла

ОЦК-фазы (bcc) либо ГЦК-фазы (fcc). Видно так-

552

ЖЭТФ, том 156, вып. 3 (9), 2019

Свойства кристаллов и жидкости Юкавы.. .

n_bcc/nliq

nfcc /nliq

1.030

1.025

1.020

1.015

1.010

1.005

1.000

10^-1

10⁰

10¹

10²

10^-1

10⁰

10¹

10²

T.104

Температурные зависимости отношения плотностей кристаллов и жидкости на линии плавления: a — ОЦК, б — ГЦК,

1 — численный расчет при T = const, 2 — аппроксимация результатов, полученных при T = const, 3 — численный расчет

при u0 = const, 4 — аппроксимация результатов, полученных при u0 = const

же, что соотношение плотностей слабо зависит от

стемы считалась неизменной и равной 300 К. Таким

выбора способа изменения безразмерной температу-

образом, поскольку Γ^∗ ∝

T-1, следует ожидать, что

ры

T. Заметим, что при слабой зависимости отно-

при постоянных температуре и объеме рост

T приве-

шения плотностей от температуры и давления абсо-

дет к уменьшению давления, а при постоянном дав-

лютное значение разности объемов Δv = v_cr - v_liq

лении — к уменьшению объема, т. е. α < 0.

меняется в выбранном диапазоне параметров весь-

Аналогичный результат может быть получен

ма существенно: давлению p = 0.025 соответствует

аналитически, при рассмотрении почти идеального

величина Δv ≈ 20, а для давления p = 12.8 харак-

газа частиц Юкавы. Заметим, что хотя такая систе-

терно значение Δv ≈ 0.02, что примерно в 1000 раз

ма не является моделью слабонеидеальной больц-

меньше.

мановской плазмы [54], ее рассмотрение важно для

По результатам исследования температурной за-

понимания термодинамики юкавских систем.

висимости энтропии и объема по формулам (4) были

рассчитаны значения коэффициента теплового рас-

Функция радиального распределения почти иде-

ширения и теплоемкости, а также по формулам (5)

ального газа может быть записана в виде [54]

соответствующие значения эффективных парамет-

ров. Результаты представлены в табл. 4.

(

)

Из табл. 4 видно, что значения параметров α и

g(r) ≈ exp

≈1-

(19)

-k_BT

k_BT

cp, вычисленные при постоянном значении констан-

ты связи u₀ по формуле (4), существенно отличают-

ся от значений параметров α и ĉ_p, рассчитанных при

Давление вычисляется с использованием вириаль-

постоянной температуре по формуле (5). Более то-

ного уравнения состояния (18). Подставляя выра-

го, значения коэффициента теплового расширения

жение (19) в (18) и интегрируя, получим

α положительны, в то время как α < 0 во всем рас-

смотренном диапазоне параметров. Результат для

[

(

)]

α находится в качественном соответствии с резуль-

2πu₀n

u₀

p = nk_BT

(20)

татами [44], где исследовалась зависимость термо-

k_BTk2D

4k_BTλ_D

динамических свойств жидкости Юкавы от Γ^∗ при

постоянном объеме и было показано, что давление

жидкости — монотонно возрастающая функция Γ^∗.

Дифференциал выражения (20) при постоянном

Заметим, что в ходе исследований температура си-

давлении равен нулю:

553

В. В. Решетняк, А. В. Филиппов

ЖЭТФ, том 156, вып. 3 (9), 2019

Таблица 4. Параметры системы Юкавы на линии

Из уравнения (21) следует, что

плавления

)

κ³ +

T-2

⁽∂n

ОЦК

∂T_p

T κ³ +4

T-1 - π^˜-2

0,u0

(

)

ĉp

(22)

(

)

πnk_D

T-1

∂n

0.025

3.043 · 10⁴

-5.453 · 10⁴

2.956

0.609

(

∂u₀

-2

p0,T

k_BT κ³ + 4

T-1 -

0.05

4.743 · 10³

-1.583 · 10⁴

2.184

0.631

0.10

1.043 · 10³

-5.625 · 10³

1.870

0.682

Выражения (22) свидетельствуют о том, что

при постоянном объеме в почти идеальной системе

0.20

3.371 · 10²

-2.156 · 10³

2.320

0.841

Юкавы плотность является монотонно убывающей

0.40

1.126 · 10²

-1.179 · 10³

1.998

0.923

функцией как температуры, так и константы связи.

0.80

4.632 · 10¹

-6.979 · 10²

2.010

0.983

Подставляя (22) в (4) и (5), получим

1.60

2.281 · 10¹

-4.847 · 10²

1.947

1.042

-2

κ³ +

3.20

1.021 · 10¹

-3.144 · 10²

1.944

1.166

α=

> 0,

6.40

5.652

-2.508 · 10²

1.874

1.190

T κ³ +4

T-1 -

T-2

(23)

12.8

2.890

-1.762 · 10²

1.891

1.272

T-1

α=-

< 0.

T² κ³ + 4

T-1 -

T-2

ГЦК

Неравенства (23) справедливы для почти идеально-

0.025

2.283 · 10⁴

-5.390 · 10⁴

2.067

0.612

го газа в силу малости

T-1. Заметим, что поскольку

0.05

4.392 · 10³

-1.589 · 10⁴

1.995

0.643

значение параметра неидеальности также определе-

0.10

1.033 · 10³

-5.619 · 10³

1.938

0.707

но с использованием соотношения константы свя-

0.20

3.045 · 10²

-2.172 · 10³

1.975

0.832

зи и температуры, выводы настоящей работы, по-

лученные для

T, справедливы и для Γ.

0.40

1.128 · 10²

-1.184 · 10³

1.908

0.920

0.80

4.550 · 10¹

-6.923 · 10²

1.936

1.074

1.60

2.240 · 10¹

-4.818 · 10²

1.878

1.086

4. ЗАКЛЮЧЕНИЕ

3.20

1.027 · 10¹

-3.149 · 10²

1.922

1.215

6.40

5.735

-2.509 · 10²

1.925

1.214

В работе исследованы свойства жидкости, ГЦК-

12.8

5.411

-1.768 · 10²

2.216

1.356

и ОЦК-кристаллов Юкавы вблизи линии плавле-

Жидкость

ния. Вычислены давления, температуры и химичес-

кие потенциалы систем на линии фазового равнове-

0.025

2.114 · 10⁴

-5.490 · 10⁴

1.920

0.580

сия. Определены энтропии и плотности разных фаз.

0.05

4.296 · 10³

-1.590 · 10⁴

1.939

0.764

Показано, что вдоль линии плавления энтропия

0.10

1.055 · 10³

-5.604 · 10³

1.882

0.697

кристаллов и жидкостей слабо меняется и флук-

туирует относительно некоторых средних значений,

0.20

2.943 · 10²

-2.180 · 10³

2.010

0.815

причем дисперсия этих флуктуаций для кристаллов

0.40

1.112 · 10²

-1.183 · 10³

1.864

1.041

не превышает 3 %, а для жидкости примерно равна

0.80

4.503 · 10¹

-6.971 · 10²

1.966

1.012

4.2 %. При этом энтропия кристаллов ниже, чем у

1.60

2.229 · 10¹

-4.833 · 10²

1.908

1.083

жидкости, примерно на 30 % (относительно энтро-

3.20

1.011 · 10¹

-3.164 · 10²

1.986

1.191

пии жидкости). По разности энтропий была вычис-

лена скрытая теплота плавления кристаллов, кото-

6.40

5.645

-2.497 · 10²

1.871

1.372

рая положительна как для ОЦК-фазы, так и для

12.8

2.865

-1.773 · 10²

1.909

1.353

ГЦК-фазы.

Плотности разных фаз на линии плавления сла-

(

)

πn²

u²⁰

бо различаются: не более чем на 3 % для выбранного

dp ≡

0 = nk_B +

dT +

диапазона параметров. Вдоль линии плавления наб-

2k_D k_BT²

)

людается монотонно убывающая логарифмическая

⁽2πn

πn²

u₀

du₀ +

температурная зависимость отношения плотностей

k2D

k_D k_BT

(

)

n_cr/n_liq, для которой предложены аппроксимирую-

4πnu₀

πnu²⁰

щие выражения.

+ k_BT +

dn.

(21)

k2D

k_Dk_B

554

ЖЭТФ, том 156, вып. 3 (9), 2019

Свойства кристаллов и жидкости Юкавы.. .

Проведенный в работе анализ свидетельствует о

В. Н. Цытович, УФН 167, 57 (1997).

неоднозначной зависимости термодинамических па-

В. Е. Фортов, А. Г. Храпак, С. А. Храпак,

раметров системы Юкавы от безразмерной темпера-

В. И. Молотков, О. Ф. Петров, УФН 174, 495

туры

T. Эта неоднозначность обусловлена возмож-

(2004).

ностью изменения

T разными способами, при ко-

торых размерная температура и параметры потен-

I. Mann, N. Meyer-Vernet, and A. Czechowski, Phys.

Rep. 536, 1 (2014).

циала взаимодействия могут либо изменяться, либо

фиксироваться. Физический смысл отмеченных спо-

D. Samsonov, S. K. Zhdanov, R. A. Quinn, S. I. Po-

собов изменения безразмерной температуры разли-

pel, and G. E. Morfill, Phys. Rev. Lett. 92, 255004

чен: в то время как изменение T при постоянных

(2004).

значениях параметров потенциала означает нагрев

S. Nunomura, D. Samsonov, S. Zhdanov, and G. Mor-

вещества, изменяя параметры потенциала при по-

fill, Phys. Rev. Lett. 95, 025003 (2005).

стоянной температуре, мы получаем разные веще-

ства при одинаковой температуре. В экспериментах

A. I. Momot, A. G. Zagorodny, and O. V. Momot,

с пылевой плазмой оба сценария могут быть легко

Phys. Plasmas 25, 073706 (2018).

реализованы.

B. Liu and J. Goree, Phys. Rev. Lett. 100, 055003

Существенные различия, связанные с разными

(2008).

способами изменения безразмерной температуры

S. Ratynskaia et al., Phys. Rev. Lett. 96, 105010

при исследовании температурных зависимостей тер-

(2006).

модинамических величин, наблюдаются при вычис-

лении положения точек фазового равновесия на

10.

G. E. Morfill et al., Phys. Rev. Lett. 83, 1598 (1999).

фазовой диаграмме, а также плотности и энтропии

11.

A. V. Filippov, V. N. Babichev, A. F. Pal’, A. N. Sta-

в них. При исследовании линейного отклика плот-

rostin, and V. E. Cherkovets, Contrib. Plasma Phys.

ности и энтропии на малое изменение безразмерной

56, 286 (2016).

температуры при постоянных значениях парамет-

ров потенциала нами были вычислены коэффициент

12.

B. Klumov et al., Europhys. Lett. 92, 15003 (2010).

теплового расширения α и теплоемкость при посто-

13.

B. Klumov et al., Plasma Phys. Control. Fusion 51,

янном объеме c_p. Кроме того, нами были проанали-

124028 (2009).

зированы соответствующие эффективные значения

параметров α и ĉ_p, определенные при постоянной

14.

M. Rubin-Zuzic et al., Nature Phys. 2, 181 (2006).

температуре и переменном значении константы u₀

15.

U. Konopka, L. Ratke, and H. M. Thomas, Phys. Rev.

в потенциале (1). Вычисленные таким образом па-

Lett. 79, 1269 (1997).

раметры различны. Показано, что величина α по-

ложительна и вдоль линии плавления монотонно

16.

U. Konopka, G. E. Morfill, and L. Ratke, Phys. Rev.

убывает с ростом давления и температуры, в то

Lett. 84, 891 (2000).

время как величина α отрицательна и ее темпера-

17.

O. S. Vaulina, O. F. Petrov, A. V. Gavrikov, and

турная зависимость монотонно возрастает. Для сла-

V. E. Fortov, Plasma Phys. Rep. 33, 278 (2007).

бонеидеального газа Юкавы аналогичный результат

18.

О. С. Ваулина, Е. А. Лисин, А. В. Гавриков,

получен аналитически. Выводы настоящей статьи

О. Ф. Петров, В. Е. Фортов, ЖЭТФ 137, 751

о неоднозначности определения состояния системы

(2010) [O. S. Vaulina, E. A. Lisin, A. V. Gavrikov,

Юкавы безразмерной температурой

T справедливы

O. F. Petrov, and V. E. Fortov, J. Exp. Theor. Phys.

также и для параметра неидеальности Γ.

110, 662 (2010)].

Финансирование. Работа выполнена при

19.

J. M. e Silva and B. J. Mokross, Phys. Rev. B 21,

2972 (1980).

поддержке Российского научного фонда (проект

№16-12-10424).

20.

D. Hone, S. Alexander, P. M. Chaikin, and P. Pincus,

J. Chem. Phys. 79, 1474 (1983).

21.

K. Kremer, M. O. Robbins, and G. S. Grest, Phys.

ЛИТЕРАТУРА

Rev. Lett. 57, 2694 (1986).

1. V. E. Fortov, A. V. Ivlev, S. A. Khrapak, A. G. Khra-

22.

M. O. Robbins, K. Kremer, and G. S. Grest, J. Chem.

pak, and G. E. Morfill, Phys. Rep. 421, 1 (2005).

Phys. 88, 3286 (1988).

555

В. В. Решетняк, А. В. Филиппов

ЖЭТФ, том 156, вып. 3 (9), 2019

23.

E. J. Meijer and D. Frenkel, J. Chem. Phys. 94, 2269

39.

O. S. Vaulina and X. G. Koss, Phys. Rev. E 92,

(1991).

042155 (2015).

24.

R. T. Farouki and S. Hamaguchi, Appl. Phys. Lett.

40.

Н. П. Коваленко, И. З. Фишер, УФН 108, 209

61, 2973 (1992).

(1972) [N. P. Kovalenko and I. Z. Fisher, Phys. Usp.

15, 592 (1973)].

25.

C. F. Tejero, J. F. Lutsko, J. L. Colot, and M. Baus,

Phys. Rev. A 46, 3373 (1992).

41.

S. Plimpton, J. Comput. Phys. 117, 1 (1995).

26.

G. Dupont, S. Moulinasse, J. P. Ryckaert, and M. Ba-

42.

LAMMPS http://lammps.sandia.gov.

us, Mol. Phys. 79, 453 (1993).

43.

W. Shinoda, M. Shiga, and M. Mikami, Phys. Rev.

27.

R. T. Farouki and S. Hamaguchi, J. Chem. Phys. 101,

B 69, 134103 (2004).

9885 (1994).

44.

В. В. Решетняк, А. Н. Старостин, А. В. Филиппов,

28.

S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, J.

ЖЭТФ 154, 1258 (2018).

Chem. Phys. 105, 7641 (1996).

45.

L. Verlet, Phys. Rev. 159, 98 (1967).

29.

S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin,

46.

D. Frenkel and B. Smit, Understanding Molecular Si-

Phys. Rev. E 56, 4671 (1997).

mulation: from Algorithms to Applications, Vol. 1,

30.

S. Hamaguchi, Plasmas Ions 2, 57 (1999).

Elsevier (2001).

31.

G. P. Hoffmann and H. Löwen, J. Phys.: Condens.

47.

Л. Д. Ландау, Е. М. Лифшиц, Статистическая

Matter 12, 7359 (2000).

физика, ч. 1, т. 5, Наука, Москва (1995).

32.

O. S. Vaulina and S. A. Khrapak, ЖЭТФ 117, 326

48.

T. J. H. Vlugt, J. P. J. M. Van der Eerden,

(2000) [O. S. Vaulina and S. A. Khrapak, J. Exp.

M. Dijkstra, B. Smit, and D. Frenkel, Introduc-

Theor. Phys. 90, 287 (2000)].

tion to Molecular Simulation and Statistical Thermo-

dynamics, ebook (2008).

33.

О. С. Ваулина, С. А. Храпак, ЖЭТФ 119, 264

(2001).

49.

T. Morita, Prog. Theor. Phys. 23, 829 (1960).

34.

O. S. Vaulina, S. V. Vladimirov, O. F. Petrov, and

50.

A. V. Filippov, A. N. Starostin, I. M. Tkachenko, and

V. E. Fortov, Phys. Rev. Lett. 88, 245002 (2002).

V. E. Fortov, Phys. Lett. A 376, 31 (2011).

35.

M. J. Stevens and M. O. Robbins, J. Chem. Phys.

51.

A. V. Filippov, A. N. Starostin, I. M. Tkachenko, and

98, 2319 (1993).

V. E. Fortov, Contrib. Plasma Phys. 53, 442 (2013).

36.

R. S. Hoy and M. O. Robbins, Phys. Rev. E 69,

52.

K.-C. Ng, J. Chem. Phys. 61, 2680 (1974).

056103 (2004).

53.

N. H. March and M. P. Tosi, Introduction to Liquid

37.

B. A. Klumov, M. Rubin-Zuzic, and G. E. Morfill,

State Physics, World Sci., New Jersey-London-Sin-

JETP Lett. 84, 542 (2007).

gapore-Hong Kong (2002).

38.

Б. А. Клумов, УФН 180, 1095 (2010) [B. A. Klumov,

54.

Л. П. Кудрин, Статистическая физика плазмы,

Phys. Usp. 53, 1053 (2010)].

Атомиздат, Москва (1974).

556