Manifestation of the hexatic phase in confined two-dimensional systems with circular symmetry

Full Text

ВВЕДЕНИЕ

Планарные системы одноименно заряженных частиц, взаимодействующих посредством кулоновского потенциала и запертых внешними потенциалами с высокой симметрией, играют важную роль в различных областях как экспериментальной, так и теоретической физики и химии. Примерами таких систем являются однокомпонентная плазма [1–4], диполярные коллоидные системы [5], самоорганизация коллоидных частиц на границе раздела двух жидкостей [6, 7], планарные квантовые точки [8]. В квазидвумерных системах изотропное отталкивание между одноименно заряженными частицами приводит к формированию гексагональной решетки, в которой шесть соседних частиц расположены симметрично относительно выбранной частицы. Эта решетка является аналогом трехмерного кристалла Вигнера [9] на плоскости [10]. Однако очень часто энергетически выгодным оказывается состояние с нарушенной симметрией, с избыточными плотностями дисклинационных дефектов. Напомним, что дисклинации представляют собой точечные дефекты вращения, ориентированные вдоль связи, а в гексагональной решетке они принимают форму точек, нарушающих шестикратную упаковку (например, [11]).

Очевидно, что большой интерес вызывают условия, при которых изменяется симметрия системы при изменении как внутренних, так и внешних параметров – температуры, плотности, магнитного поля. Изучение фазовых переходов является ключевым элементом в понимании поведения ферромагнетиков и сверхпроводников. Фазовые переходы в системах заряженных частиц, ограниченных потенциалом с круговой симметрией, представляют особый интерес в физике конденсированного состояния при создании наноструктурных объектов с конечным числом частиц (например, [8]).

Задача определения типа фазового перехода в планарных системах до сих пор не решена. Большое количество исследований посвящено плавлению гексагональной решетки при увеличении температуры системы и определению температуры соответствующего фазового перехода. Более того, исследователи описанных выше систем заряженных частиц не пришли к единому мнению о типе фазового перехода при увеличении температуры. Тем не менее одним из общепринятых сценариев является фазовый переход типа гексагональная –гексатическая фаза–изотропная квазижидкость (например, обзор [12]).

В соответствии с теорией Березинского– Костерлица–Таулеса (БКТ) [13–16] и результатами теоретических и экспериментальных исследований [10, 12, 17] при низких, близких к абсолютному нулю, температурах гексагональная решетка является энергетически наиболее выгодной в планарных системах классических, одноименно заряженных точечных частиц, связанных кулоновским взаимодействием. В то же время в системах с конечным числом таких частиц, запертых посредством круговых потенциалов с бесконечными стенками, равновесная конфигурация определяется необходимостью образования гексагональной решетки, с одной стороны, и кольцеобразной структуры, с другой [18–20]. Эта конкуренция приводит к образованию в таких системах внутренних дефектов в отличие от случая неограниченных областей, когда в основном состоянии системы нет дефектов. Симметрия гексагональной решетки нарушается при приближении к круговой границе. При числе частиц N >> 1000 случаю “изотропной квазижидкости” отвечает классическое равномерное распределение заряда по поверхности планарной системы.

Основной задачей настоящей работы был поиск фазовых характеристик и параметров порядка − предвестников фазового перехода типа гексагональная–гексатическая фаза в зависимости от числа частиц в системе при нулевой температуре. Ситуация существенно осложняется конечностью числа частиц, локализованных в ограниченной области вследствие наличия бесконечного, типа “твердая стенка”, потенциала запирания на границе и появления множества метастабильных состояний.

В последние годы компьютерное моделирование стало важным инструментом для изучения фазовых переходов. Оно позволяет изучать поведение систем в различных условиях, что часто невозможно в экспериментах. Кроме того, компьютерные исследования дают возможность уточнить теоретические модели и предсказания. Поэтому в качестве “экспериментальной базы” в работе использован численный подход [20], основанный на решении нелинейных вариационных уравнений для конечного числа одноименно заряженных частиц, запертых круговыми потенциалами с высокой симметрией [19, 20]. Отметим, что подобный подход позволяет значительно ускорить поиск равновесного состояния по сравнению с традиционными методами молекулярной динамики и подходами типа Монте-Карло.

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

ФИЗИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ

Дана система из N одноименно заряженных частиц с кулоновским взаимодействием и двумерным ограничивающим потенциалом радиуса R. Гамильтониан такой системы записывается следующим образом (1):

$H=\sum _{i=1}^{N}V\left(r_i\right)+\alpha \underset{i<j}{\overset{N}{\sum _{i,j=1}}}\frac{1}{\left|{\mathbf{r}}_i-{\mathbf{r}}_j\right|}+\sum _{i=1}^N{T}_i,$ (1)

где r_i = | r_i | – это расстояние до центра области, ограниченной потенциалом, α = e²/4πε₀ε_r – величина, характеризующая силу взаимодействия зарядов в среде, T_i – кинетическая энергия частицы. Ограничивающий потенциал V(r) определяется следующим образом:

$V\left(r\right)=\left\{\begin{array}{l}0, r<R\\ \infty , r\ge R\end{array}\right.$. (2)

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

$H=\sum _{i=1}^{N}V\left(r_i\right)+\alpha \underset{i<j}{\overset{N}{\sum _{i,j=1}}}\frac{1}{\left|{\mathbf{r}}_i-{\mathbf{r}}_j\right|}.$ (3)

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

МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ И СХЕМА ВЫЧИСЛЕНИЙ

Для поиска устойчивых конфигураций системы с минимумом энергии заряженных частиц, описываемых гамильтонианом (1), в настоящей работе использован подход, развитый в [21, 22]. Для моделирования эволюции описанной выше системы применяли метод молекулярной динамики [23] с учетом оценок, полученных при решении вариационных уравнений. В качестве начального приближения брали распределение частиц по кольцам, вычисленное с помощью уравнения (4), которое получено в результате решения уравнений модели устойчивых конфигураций [19, 20]. Термин “кольцо” обозначает кольцевой слой радиуса r_p (детали описаны в [19]), в котором располагаются частицы после достижения системой состояния с минимумом энергии:

N(p) = 4.1988 + 7.27p^1.66 + 0.48p^2.48. (4)

Здесь p – номер кольца от центра круговой области, N(p) – число частиц на p-м кольце. Полученная конфигурация использована для запуска молекулярной динамики с торможением. При достижении нулевой температуры системы расчеты, реализующие ее временну́ю эволюцию, заканчиваются, и рассчитывается энергия получившейся устойчивой конфигурации частиц в системе. Последний этап реализован аналогично традиционной схеме метода молекулярной динамики. Общее количество итераций, необходимое для достижения глобального минимума энергии, уменьшается на три–четыре порядка. Такой подход позволяет за разумное время проводить моделирование эволюции системы, состоящей из нескольких тысяч заряженных частиц. Например, для достижения основного состояния системы с минимумом энергии для 8000 частиц (рис. 1) методом молекулярной динамики с использованием технологии CUDA потребовалась 841 с на компьютерной системе, состоящей из процессора Intel(R) Xeon(R) Gold 6148 CPU @ 2.40 GHz и графических процессоров Tesla v100-sxm2-32gb.

 

Рис. 1. Распределение дефектов с ненулевым ориентационным параметром порядка для системы точечных электронов в круговой области с бесконечной жесткой границей. Темными точками обозначены электроны с параметром порядка Ψ₆(r_k) = 0, светлыми – электроны с ненулевым параметром порядка.

 

ТОПОЛОГИЧЕСКИЙ ФАЗОВЫЙ ПЕРЕХОД

В соответствии с теорией БКТ, фазовый переход гексагональная–гексатическая фаза в описанных выше двумерных системах является топологическим фазовым переходом, т.е. переходом между двумя топологическими фазами. Топологическая фаза материи – это состояние двумерной системы из большого числа сильно взаимодействующих частиц (конденсированной среды), характеризуемой определенным топологическим инвариантом. В рассматриваемом случае топологическим инвариантом является тип упорядочения частиц, описываемый определенными группами симметрии [11]. Гексагональная решетка характеризуется трансляционной и ориентационной дискретной симметрией шестого порядка (диэдральной группой симметрии шестого порядка). Как следствие, в гексагональной фазе реализуется так называемый дальний порядок. В гексатической фазе трансляционная симметрия нарушается за счет топологических дефектов, тогда как ориентационная симметрия шестого порядка может сохраняться. Поэтому в гексатической фазе реализуется так называемый квазидальний порядок. Таким образом, переход из одной фазы в другую реализуется вследствие появления топологических дефектов, нарушающих трансляционную симметрию периодической гексагональной структуры.

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

ОРИЕНТАЦИОННЫЙ ПАРАМЕТР ПОРЯДКА

Одним из возможных кандидатов на роль параметра порядка для фазового перехода гексагональная–гексатическая фаза может служить так называемый ориентационный порядок связи [24]:

${\mathrm{\psi }}_6\left(r_k\right)=\frac{1}{N_b}\sum _{l=1}^{N_b}\exp \left(i6{\mathrm{\theta }}_{kl}\right).$ (5)

Здесь N_b – число ближайших соседей, θ_kl – угол между направлением связи между частицами k и l и произвольно выбранным направлением. Параметр ψ₆ является важнейшей характеристикой двумерных систем частиц, образующих кристаллическую решетку с симметрией шестого порядка. Параметром порядка в соответствии с теорией фазовых переходов Л.Д. Ландау [25, 26] в этом случае может служить величина, которая может быть названа ориентационным параметром порядка:

Ψ₆₍r_k) = 1 – ψ₆(r_k) – ε₆. (6)

Здесь ε₆ – точность определения степени ориентационной симметрии шестого порядка, связанный с точностью вычисления функции ψ₆(r_k). В этом случае в узлах гексагональной решетки параметр порядка будет равен нулю (ψ₆(r_k) = 0). Нарушение гексагональной симметрии будет характеризоваться не равным нулю параметром порядка. Пример визуализации параметра порядка для числа точечных электронов N = 8000 с точностью ε₆ = 0.2 приведен на рис. 1.

Как видно из рисунка, нарушения гексагональной симметрии проявляются в виде дефектов решетки, имеющих форму цепочек и островков из частиц с ненулевым параметром порядка. Наиболее плотное скопление дефектов наблюдается вблизи границы. Эта картина существенно отличается от равномерного распределения дефектов в случае бесконечного числа одинаково заряженных частиц в неограниченной или даже достаточно большой области на плоскости. Отметим, что энергия рассмотренного состояния для N частиц с зарядом (–e) практически (с точностью до 2%) совпадает с энергией $E=\frac{\pi }{4}{N}^2$ однородно распределенного по поверхности диска заряда Q = eN. В случае гармонического запирания также наблюдается скопление дефектов вблизи границы, но гораздо менее плотное, чем в случае потенциала запирания (2).

ТОПОЛОГИЧЕСКИЙ ЗАРЯД

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

В соответствии с теорией Березинского– Костерлица–Таулеса–Хальперина–Нельсона–Янга [13–16, 24, 27, 28] существуют два типа топологических дефектов двумерных кристаллических решеток – дислокации и дисклинации, которые характеризуют нарушения трансляционной и вращательной симметрии. Дислокация в точке r определяется величиной приращения, которое получает контурный интеграл от поля смещений u(x, y) при обходе по замкнутому контуру, содержащему дислокацию [30]:

φ du = a₀b(r) = -n(r)a₀e₁ - m(r)a₀e₂, (7)

где b(r) – безразмерный вектор Бюргерса, a₀ – период рассматриваемой треугольной решетки, e₁ и e₂ – базисные векторы решетки, m и n – целые числа. Вектор Бюргерса является количественной характеристикой, описывающей искажения кристаллической решетки вокруг дислокации.

Дисклинация на гексагональной решетке характеризуется тем, что при обходе по замкнутому контуру, содержащему дисклинацию, интеграл от угла поворота

$\mathrm{\vartheta }(\mathrm{x},\mathrm{y})=\frac{1}{2}({\partial }_x{u}_y-{\partial }_y{u}_x),$

где u(x, y) = (u_x(x, y), u_y(x, y)) – вектор смещения некоторой малой области кристалла вследствие термодинамических флуктуаций при T ≠ 0, получает приращение, кратное 2π/6:

$\mathrm{\phi }d\mathrm{\vartheta }\left(\mathbf{r}\right)=-\frac{2\mathrm{\pi }}{6}s, s=\pm 1,\pm 2,...$ (8)

Величина s называется топологическим зарядом дисклинации. Дислокацию можно рассматривать как связанную пару дисклинаций с разными топологическими зарядами. Естественно предположить, что определение (8) справедливо и при нулевой температуре.

Для определения координат и формы дефектов был применен подход, основанный на подсчете количества ближайших соседей [18]. Для определения количества ближайших соседей использовали алгоритм Вороного [29], в основе которого лежит разбиение пространства на ячейки. Каждая ячейка состоит из точек, ближе всех расположенных к выбранной частице, чем к какой-либо другой частице из данного множества. При таком разбиении можно непосредственно наблюдать дефекты двух видов: дисклинации и дислокации. Дисклинацией называется ориентационный дефект с координационным числом C_d = 5 или C_d = 7. Координационное число определяется числом сторон выпуклого многоугольника, описанного вокруг частицы в разбиении Вороного [30, 31]. В случае шестиугольника координационное число будет равно C_d = 6.

Учитывая вышеизложенные факты, можно естественным образом связать топологический заряд (8) с координационным числом посредством соотношения s = 6 – C_d. Таким образом, дисклинацией называется топологический дефект, при котором топологический заряд частицы не равен нулю. Дислокацией называется дефект, состоящий из нескольких пар дисклинаций; соседние дисклинации в паре имеют равные по модулю топологические заряды с противоположными знаками. Топологический заряд узлов гексагональной решетки без дефектов в соответствии с определением равен нулю (рис. 2). Дислокации не нарушают гексагональную симметрию и, следовательно, такая пара дисклинаций не увеличивают энергию системы [18]. Общее количество дисклинаций N⁺ с топологическим зарядом s = 1 и N^– с топологическим зарядом s = –1 зависит от конкретной конфигурации.

 

Рис. 2. Распределение топологического заряда для числа частиц N: a – 92; б – 136; в – 187.

 

Напомним, что количество дисклинаций в плоской треугольной решетке, вписанной в круговую область, определяется теоремой Эйлера для треугольника и не может быть изменено. В этом случае чистый топологический заряд N⁺ – N^– всегда равен шести [2, 27] (строгое доказательство приведено в [2]).

ОРИЕНТАЦИОННАЯ КОРРЕЛЯЦИОННАЯ ФУНКЦИЯ

При описании фазовых переходов в двумерных системах одной из важнейших характеристик является поведение ориентационной корреляционной функции между двумя частицами i и j на расстоянии r_ij = r_i – r_j, которая определяется следующим образом:

 (9)

Здесь угловые скобки обозначают среднее статистическое значение для всех пар частиц на расстоянии | r_ij |.

Для поиска сигналов фазового перехода гексагональная–гексатическая фаза исследуем более детально распределение топологических зарядов и поведение корреляционной функции (9) для системы заряженных частиц с N = 92, 136, 187.

На рис. 2 целыми числами 0, ±1, ±2, … обозначен топологический заряд каждой частицы, который связан по определению с количеством сторон выпуклого многоугольника, внутри которого находится частица. Нулевому топологическому заряду соответствует шестиугольник, описанный вокруг частицы, вершинами которого являются ее ближайшие соседи. Однако это еще не означает, что частица находится в узле гексагональной (составленной из правильных шестиугольников) решетки. Частица принадлежит узлу гексагональной решетки только в том случае, когда ее ориентационный параметр порядка Ψ₆(r_k) (6) с абсолютной точностью равен нулю или ориентационный порядок связи ψ₆(r_k) (5) равен единице. Количество частиц, для которых ориентационный порядок связи строго равен единице, можно увидеть на гистограммах распределений ориентационного порядка связи (рис. 3).

 

Рис. 3. Распределение ориентационного порядка связи ψ₆ (5) для числа частиц N: a – 92; б – 136; в – 187.

 

Например, для N = 92 количество частиц со строго нулевым ориентационным параметром порядка равно семи. Это означает, что гексагональной решетке принадлежат только семь частиц. Аналогичная ситуация имеет место для N = 136, 187 (рис. 3б, в). Идеальной гексагональной решетке принадлежат только семь частиц в центре. Для остальных частиц гексагональная решетка начинает значительно искажаться по мере приближения к границе области. Стремление системы к наиболее выгодному энергетическому состоянию приводит к неоднородной локальной плотности. В результате формируется основное состояние, в котором гексагональная решетка изотропно расширяется и сжимается в области локализации взаимодействующих частиц (например, рис. 3а в [20]). Дальнейшее увеличение числа частиц приводит к искажению решетки и появлению дисклинаций, состоящих из частиц с ненулевым топологическим зарядом, что означает разрушение гексагональной решетки. Дисклинации образуются с равными и противоположными по знаку топологическими зарядами, которые, объединяясь в пары, образуют дислокации (комбинации точек с зарядами +1 и –1 на рис. 2а). Ориентационное искажение гексагональной решетки и образование дефектов в виде дисклинаций и дислокаций представляет собой проявление новой, так называемой гексатической фазы, которая может служить предвестником фазового перехода гексагональная–гексатическая фаза. Анализ распределения топологических зарядов указывает на доминирование дисклинаций с топологическим зарядом s = –1 c ростом числа частиц в рассматриваемой системе. Например, в случае N = 92 количество частиц с топологическим зарядом s = –1 равно 14, а с топологическим зарядом s = +1 – трем. В случае N = 187 количество частиц с топологическим зарядом s = –1 равно 34, а с топологическим зарядом s = +1 – восьми. Образование кольцевой структуры характеризуется значительным увеличением количества частиц с топологическим зарядом | s | > 2 (рис. 2). Так, в случае N = 92 количество частиц с топологическим зарядом | s | > 2 равно 59, а в случае N = 187 – 86.

Для более детального исследования фазового перехода проанализируем ориентационную корреляционную функцию G₆(r) (рис. 4). Цифрами 0 и 1 на графиках обозначены прямолинейные участки, которые отвечают гексагональной решетке. Соответствующие области на рис. 2 представлены точками с топологическим зарядом, равным нулю, характеризуемые параметром порядка связи, строго равным единице (последний бин на гистограммах рис. 3). Участки графиков между точками 1 и 2 представляют собой области гексатической фазы. Эти области соответствуют частицам, для которых начинается искажение гексагональной решетки. Как результат, это приводит к появлению дефектов в виде дисклинаций и дислокаций (комбинации точек с зарядами +1 и –1 на рис. 2). После точки 2 наблюдается резкое падение G₆, которое сменяется затухающей колебательной кривой, что свидетельствует о преобладании кольцеобразной фазы. В неограниченной области, в отличие от ограниченной круговой области с конечным числом частиц, после точки 2 функция G₆ обычно экспоненциально убывает, что свидетельствует о преобладании фазы “изотропной квазижидкости” (например, рис. 14 в [12]).

 

Рис. 4. Ориентационная корреляционная функция G₆(r) (9) для числа частиц N: a – 92; б – 136; в – 187.

 

ЗАКЛЮЧЕНИЕ

В работе предложен подход для поиска предвестников фазового перехода гексагональная–гексатическая фаза в зависимости от числа частиц при нулевой температуре в планарных системах одинаково заряженных частиц, взаимодействующих посредством кулоновского потенциала и запертых внешними потенциалами с круговой симметрией. Подход основан на совместном анализе зависимости введенного авторами ориентационного параметра порядка и ориентационной корреляционной функции от числа частиц в системе. Предложенный подход апробирован в рамках комбинированной модели [21, 22]. В этой модели с помощью решения нелинейных вариационных уравнений получены близкие к основному состоянию равновесные конфигурации, использование которых значительно ускоряет процесс поиска основного состояния в рамках молекулярной динамики.

Проведенный анализ продемонстрировал ряд характеристических особенностей предвестников фазового перехода гексагональная–гексатическая фаза в ограниченных конечных двумерных системах в зависимости от числа частиц при нулевой температуре. С ростом числа частиц гексагональная решетка, которая формируется в центральной области локализации частиц, начинает искажаться по направлению к границе. Введенный ориентационный параметр порядка позволяет установить появление дефектов решетки в результате конкуренции распределения частиц в узлах треугольной (гексагональной) решетки с кольцевой структурой, обусловленной симметрией потенциала. При дальнейшем приближении к границе решетка начинает разрушаться. В результате появляются дефекты кристаллической решетки в виде дисклинаций и дислокаций – пар дисклинаций с равными и противоположными по знаку топологическими зарядами. Дефекты такого типа и являются признаком гексатической фазы. При этом нарушается трансляционная симметрия, ответственная за дальний порядок, при сохранении ориентационного квазидальнего порядка. Более того, наличие круговой границы приводит к образованию кольцеобразных структур вблизи границы, что также влияет на характер перехода между гексагональной и гексатической фазами при росте числа частиц системы из-за конкуренции между кольцеобразными структурами и треугольной решеткой. Гексатическая фаза начинает проявляться в виде образовавшихся дефектов (дисклинаций и дислокаций) при гораздо меньшем числе частиц в случае бесконечного потенциала запирания на границе по сравнению с бесконечными неограниченными системами. С ростом числа частиц система стремится к стабильному состоянию за счет увеличения количества отрицательных топологических зарядов. Наличие потенциала с круговой симметрией приводит к образованию кольцевой структуры, которая характеризуется значительным увеличением количества топологических зарядов | s |>2.

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

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа выполнена в рамках проекта ОИЯИ № 05-6-1118-2014/2023 с использованием ресурсов гетерогенного кластера HybriLIT Лаборатории информационных технологий Объединенного института ядерных исследований. Исследования Э.Г. Никонова частично выполнены в рамках Программы фундаментальных исследований Национального исследовательского университета “Высшая школа экономики”.

Конфликт интересов. Авторы данной статьи заявляют, что у них нет конфликта интересов.

About the authors

E. G. Nikonov

Joint Institute for Nuclear Research; National Research University Higher School of Economics; Dubna State University

Author for correspondence.
Email: e.nikonov@jinr.ru
Russian Federation, 141980, Dubna; 119048, Moscow; 141980, Dubna

R. G. Nazmitdinov

Joint Institute for Nuclear Research; Dubna State University

Email: e.nikonov@jinr.ru
Russian Federation, 141980, Dubna; 141980, Dubna

P. I. Glukhovtsev

Joint Institute for Nuclear Research; Dubna State University

Email: e.nikonov@jinr.ru
Russian Federation, 141980, Dubna; 141980, Dubna

References

Supplementary files