Исследование фазовых переходов методом отжига популяции тема диссертации и автореферата по ВАК РФ 00.00.00, кандидат наук Мозоленко Вячеслав Константинович

• Мозоленко Вячеслав Константинович
• кандидат наук
• 2025

Оглавление диссертации кандидат наук Мозоленко Вячеслав Константинович

1.1 Статсумма и плотность состояний

1.2 Выжимки из теории фазовых переходов и критического поведения

1.2.1 Фазовые переходы первого и второго рода

1.2.2 Критические точки и критическое поведение

1.2.3 Конечномерное поведение

1.2.4 Намагниченность и концентрация

1.2.5 Кумулянт Биндера

1.2.6 Критическое замедление

1.2.6.1 ФП второго рода

1.2.6.2 ФП первого рода

1.3 Обзор вычислительных методов в статистической физике

1.3.1 Метод трансфер-матриц

1.3.2 Классический метод Монте-Карло

1.3.3 Метод отжига популяции

1.3.4 Метод Ванга-Ландау

2 Метод МСРА

2.1 Оригинальный алгоритм Розе-Махты

2.1.1 Проблема 1: Генерация состояния с максимальной энергией

2.1.2 Проблема 2: Верхняя часть спектра

2.1.3 Проблема 3: Знание спектра системы

2.2 Алгоритм МСРА

2.2.1 Охлаждение

2.2.2 Нагревание

2.2.3 Пересчет е(Е) в 5(Е)

2.2.4 Сшивка ветвей 5(Е)

2.2.5 Выбор нормировки энтропии

2.2.6 Расчет термодинамических функций

2.2.7 Расчет намагниченности

2.2.8 Расчет концентрации и дырочной восприимчивости

2.2.9 Поиск точки фазового перехода

2.2.10 Технические детали реализации алгоритма

2.2.11 Генерация нескореллированных конфигурация для тернарной классификации фаз в области фазового перехода первого рода

2.3 Исследование точности алгоритма

2.3.1 Ш модель Изинга

2.3.2 2Э модель Изинга

3 Сравнение метода MCPA с методом Ванга-Ландау

3.1 Модель Поттса

3.2 Детали вычислений

3.2.1 Вычислительная сложность алгоритмов

3.2.2 Выбор параметров MCPA

3.3 Результаты сравнения

3.3.1 Плотность состояний

3.3.2 Удельная теплоемкость

3.3.3 Биндеровский кумулянт

3.3.4 Распределения вероятностей энергии P(E, T)

3.3.5 Энергии упорядоченной и разупорядоченной фаз

3.4 Выводы

4 Исследование модели Блюма-Капеля

4.1 История исследования модели Блюма-Капеля

4.1.1 Оригинальные работы Блюма и Капеля

4.1.2 Трикритическая точка в смеси He3 — He4

4.1.3 Современное состояние модели

4.2 Объяснение трикритической точки с помощью теории фазовых переходов Ландау

4.3 Общие замечания

4.3.1 Общий вид симулированных параметров

4.3.2 Первые возбужденные состояния

4.3.2.1 Основное состояние

4.3.2.2 Первое возбужденное состояние

4.3.2.3 Второе возбужденное состояние

4.3.2.4 Третье возбужденное состояние

4.3.2.5 Четвертое возбужденное состояние

4.3.2.6 Пятое возбужденное состояние

4.3.3 Сравнение с решением перебором для малых размеров системы

4.3.4 Периодическая структура S при дробных значениях D

4.3.5 Число уровней энергии и время симуляции

4.3.6 Общий вид вычисляемых в MCPA параметров

4.3.7 Общий вид основных термодинамических величин

4.4 Конечномерный анализ для области фазового перехода II рода

4.4.1 Удельная теплоемкость

4.4.1.1 Детали аппроксимации удельной теплоемкости

4.4.1.2 Детали аппроксимации критической температуры

4.4.2 Биндеровский кумулянт энергии

4.4.3 Магнитная восприимчивость

4.4.4 Биндеровский кумулянт намагниченности

4.4.5 Дырочная восприимчивость

4.4.6 Максимальный наклон Биндеровского кумулянта намагниченно-

сти и логарифма степеней намагниченности

4.4.7 Сравнение оценок критических температур

4.5 Конечномерный анализ для области фазового перехода I рода

4.5.1 Энергии упорядоченной и разупорядоченной фаз

4.5.2 Латентное тепло

5 Заключение

Список литературы

Введение диссертации (часть автореферата) на тему «Исследование фазовых переходов методом отжига популяции»

1 Введение

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

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

• (полу)точные решения: метод трансфер-матриц,

• Монте-Карло (МК): классический метод МК,

• прямая оценка энтропии: метод отжига популяции, гистограммный метод, метод Ванга-Ландау.

Прямое вычисление плотности состояний (density of states, DoS, g(E)) является одним из самых удобных способов исследования физических систем, которое служит важным элементом в понимании термодинамических свойств систем с множеством состояний. Такие традиционные методы, как метод Ванга-Ландау, продемонстрировали свою эффективность, однако они часто сталкиваются с ограничениями. Например, при необходимости рассчитывать дополнительные переменные: намагниченность m и концентрацию нулевых спинов (дырок) X.

В диссертации мы представляем новый алгоритм прямой оценки энтропии, основанный на методах Монте-Карло и популяционном отжиге (microcanonical population annealing, MCPA). Этот метод сочетает в себе принципы метода Монте-Карло — статистического подхода, который использует случайные блуждания в фазовом пространстве — и концепцию популяционного отжига, которая позволяет эффективно исследовать конфигурации системы, избегая застревания в локальных минимумах энергии. Популяционный отжиг представляет собой процесс, в котором множество реплик (или конфигураций) системы одновременно "отжигаются".

В рамках диссертации мы исследуем эффективность MCPA на примерах одномерной и двумерной моделей Изинга, которые представляют собой "плодовые мушки" -тестовые системы, для которых ответ известен аналитически. Мы также рассматриваем применение данного метода к модели Блюма-Капеля, которая была нами выбрана за ее богатую фазовую диаграмму и внешнюю простоту, что позволяет глубже понять механизмы фазовых переходов и критического поведения. Кроме того, мы проводим сравнение с методом Ванга-Ландау на примере модели Поттса, что поможет выявить как сильные, так и слабые стороны нового подхода, а также его потенциальные области применения.

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

1.1 Статсумма и плотность состояний

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

• Микроканоническое состояние (или микросостояние) определяется полным набором информации про все частицы системы: координаты, скорость и т.д.

• Множество всех микросостояний системы называется фазовым пространством системы.

• В противовес микросостоянию мы также можем рассмотреть макросостояние, характеризующееся определенными значениями макроскопических параметров (температура, давление, объем и т. д.). Одному макросостоянию соответствует множество микросостояний.

Статистическая физика рассматривает положение системы в ее фазовом пространстве как случайную величину. Соответственно, каждому микросостоянию в приписывается его вероятность т = т(в), зависящая от положения в фазовом пространстве (состояния) в.

В равновесии плотность вероятности т подчиняется распределению Гиббса:

Цв) а е-вЕ(в),

где в - обратная температура, а Е(в) - энергия системы в состоянии в. Зная т, мы можем выразить среднее значение (А) любой термодинамической величины А как:

(А) = | Е т(в)А(в), (1)

2 = Е «>М, (2)

в

где 2 - нормирующая константа.

Если наша термодинамическая величина А является функцией энергии системы Е или же допускает усреднение вида А(Е) = 9Щ Е(в)=Е А(в), то мы можем выразить А: 9

(А) = 2 Е 9(Е)А(Е)е-вЕ, (3)

2

Е

2 = £ 9(Е )е-вЕ. (4)

Е

Величину g(E) = ^e(s)=e 1, которая представляет собой число состояний с энергией E, называют плотностью состояний (density of states, DoS). Знание g(E) дает нам возможность рассчитать Z для любой температуры. Также в этой работе мы будем часто упоминать величину S(E) = ln g(E), которую мы будем называть энтропией.

Величину Z также называют статсуммой (нем. Zustandssumme). Она чрезвычайно полезна тем, что ее знание дает нам полную информацию о термодинамических свойствах системы [2]. Приведем несколько примеров: среднюю энергию E, теплоемкость системы cv, энтропию S и свободную энергию F можно вычислить как производные статсуммы. Именно на этапе вычисления статсуммы в наших вычислениях наконец-то появляется температура системы и начинается изучение термодинамики.

(E > =

д ln Z

дв '

(5)

cv =

д(E> _ 1 дТ = кв Т2

(SE 2>,

(6)

д дF

S Ws ln Ws = кв (ln Z + в<E>) = дТ (квТ ln Z) =

дТ

(7)

^ = (Е) - Т£ = -квТ 1п (8)

Вычислительные алгоритмы, направленные на непосредственный расчет значений термодинамических величин (обыкновенно это методы Монте-Карло), и алгоритмы, направленные на оценку плотности состояний (алгоритм Ванга-Ландау, метод гистограмм), значительно отличаются в итоговой производительности около точек фазового перехода. Дело в феномене критического замедления для фазовых переходов второго рода [6] и в феномене образования доменной стенки для фазовых переходов первого рода[2, 6]. Эти феномены различны по своей природе, однако оба ведут к катастрофичному замедлению традиционных Монте-Карло алгоритмов вблизи точки фазового перехода.

Методы оценки плотности состояний (они же микроканонические методы) как будто бы не испытывают подобного замедления либо же испытывают его в значительно меньшей степени [7, 8]. Именно таким методом является представленный в работе метод МСРА.

Также отметим, что метод МСРА позволяет естественным образом включать расчет таких параметров, как намагниченность и концентрация дырок в алгоритм как функции энергии т = т(Е), X = X(Е), что в свою очередь позволяет очень удобно вычислять термодинамику этих параметров. Это преимущество значительно отличает МСРА от, например, алгоритма Ванга-Ландау.

1.2 Выжимки из теории фазовых переходов и критического поведения

Фазовые переходы (ФП) представляют собой явления, при которых система изменяет свое состояние (или фазу) при изменении внешних условий, таких как температура, давление или магнитное поле. Эти переходы могут быть как первого, так и второго рода, в зависимости от характера изменений в системе.

1.2.1 Фазовые переходы первого и второго рода

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

В отличие от них, фазовые переходы второго рода происходят без скрытой теплоты, и их характеристиками являются непрерывные изменения термодинамических свойств, таких как намагниченность или энергия. Фазовые переходы второго рода сопровождаются изменением симметрии вещества, и это изменение может быть количественно характеризовано параметром порядка. Параметр порядка для разных систем определяется по-разному, например: в переходе ферромагнетик-парамагнетик, у ферромагнетика имеется векторный параметр порядка - спонтанная намагниченность т. В отличие от фазового перехода первого рода, в точке ФП II рода состояния обеих фаз совпадают [2].

1.2.2 Критические точки и критическое поведение

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

Размер флуктуаций можно описать с помощью корреляционной функции О(г), которая описывает скоррелированность элементов. Для спиновой системы корреляционная функция будет выглядеть следующим образом:

ад = адмго + г)), (9)

где в(г) - значения спина в узле с координатой г. Вдали от точки фазового перехода и на больших г О ведет себя как:

О(г) а г-^ехр(-г/£), (10)

где ё - размерность системы, а £ - корреляционная длина, которая и характеризует масштаб флуктуаций в системе [6]. Вблизи же фазового перехода второго рода корреляционная длина, а также некоторые термодинамические величины становятся

бесконечными. Законы, по которому они расходятся, оказываются следующими [6]:

Корреляционная длина: £ а |т-,

Намагниченность: т а |т|в,

дт

Магнитная восприимчивость: х = —— а |т|

дН

Теплоемкость: С = а |т|"

где т =1 - Т/Тсгй - приведенная температура, Н - внешнее магнитное поле, а коэффициенты называют критическими показателями. Оказывается, что критические показатели являются взаимосвязанными и принимают некоторый дискретный набор значений, которые называют классами универсальности.

1.2.3 Конечномерное поведение

Все вышеописанное применимо к системам в термодинамическом пределе Ь ^ то. Однако в вычислительной физике мы оперируем системами значительно более скромных размеров. И в системах малых размеров возникает "псевдокритическое" (конечномерное критическое) поведение.

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

Ь « £ а |т|т| а Ь-1/^, (11)

где т =1 - Т/Тсгй - приведенная температура.

Отсюда же следует, что при приближении к точке фазового перехода максимальные значения термодинамических величин будут также достигнуты раньше точки фазового перехода. Их максимальные значения мы тогда можем оценить как [6]:

Мтах а Ь-^

Хтах а С пг Та/"

Стах а Ь

Ь^ (12)

1.2.4 Намагниченность и концентрация

В работе нам также понадобятся определения для намагниченности т и концентрации дырок X.

Рассмотрим спиновую систему размера N = Ь2 со спином 1. Каждый конкретный спин на решетке может принимать одно из трех значений: а € [-1, 0,1].

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

1

т = N

Е

а»

(13)

и концентрацию X как:

11

X = ¿.,0 = 1 - ¿Е (14)

1.2.5 Кумулянт Биндера

При определении фазовых переходов также часто используется кумулянт Биндера намагниченности (БКМ) ВСт:

ВС =1 (т4 (Т)) (15)

ВСт = 1 - 3(т2(Т))2 ■ (15)

В пределе Ь ^ ж, ВСт ^ 0 при Т > Тс и ВСт ^ 2/3 при Т < Тс. Для "достаточно больших" размеров систем функции температуры ВСт (Т) пересекаются в критической точке, которая не зависит от размера системы Ь [6, 9]. Исследуя пересечения линий БКМ, мы сможем получить еще одну оценку положения критических точек, причем она не должна иметь конечномерной зависимости от Ь.

Мы также можем определить кумулянт Биндера энергии (БКЭ) ВСе:

ВСЕ = 1 - ' (16)

который имеет минимум в критической точке. Однако, в отличие от БКМ, БКЭ имеет конечномерную зависимость от Ь. Для области фазового перехода первого рода это [6, 10]:

ТВСЕ (Ь) = ТВСЕ (ж) + ВвсЕ/Ь2, (17)

а для области фазового перехода второго рода, на общих основаниях:

ТВСЕ (Ь) = ТВСЕ (ж) + ВвсЕ/Ь1^ ■ (18)

1.2.6 Критическое замедление

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

1.2.6.1 ФП второго рода Для окрестности ФП второго рода проблемой становится расходящаяся корреляционная длина £ а |т -, которая растет вплоть до размера системы (до бесконечности в термодинамическом пределе). Из-за этого динамика системы становится "глобальной", а время, за которое система "забывает" свои начальное состояние, увеличивается как £^ ~ И-^, где г - динамический критический показатель [6, 11].

Примером такого поведения может быть модель Изинга при Т ^ Тс, в которой магнитные домены растут до размера системы, а их переструктурирование требует огромного времени.

1.2.6.2 ФП первого рода В окрестности ФП первого рода проблемой становится энергетический барьер между фазами. Предположим, у нас есть метастабильное состояние из двух фаз, разделенных фазовым барьером (интерфейсом). Этот интерфейс обладает собственной энергией Д^, который система должна преодолеть, чтобы окончательно оказаться в одной из фаз. Энергия интерфейса зависит линейно от размера

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

ехр (тёг1) [6].

Примером такой ситуации может служить модель Поттса с д > 4 (с ФП первого рода), где интерфейс между разными фазами имеет высокую энергию, и система застревает в метастабильных состояниях, пока не произойдет гигантский переход через барьер.

1.3 Обзор вычислительных методов в статистической физике

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

1.3.1 Метод трансфер-матриц

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

Для модели Блюма-Капеля с гамильтонианом в форме:

Н = -3 ^ ^ + Д ^ а2,

{¡-л > *

где 3 - константа ферромагнитной связи, Д - анизотропное поле, а - спины, равные 1, на двумерной решетке М х Ь с периодическими граничными условиями, которую мы станем рассматривать как набор М замкнутых цепочек длины Ь с периодическими

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

м

я=Е V++1+2V;+0,

3=1 4 J

где V - элемент энергии внутри цепочки (г - индекс внутри цепочки, ] - индекс цепочек), а W - между цепочек:

ь

V' = а»,3а»+1,3 + Аа?,з

»=1

ь

Щ,3+1 = а»,3+1

»=1

Тогда статсумму 2 можно записать как

М Г в„ _ в "

2 = ЕПехр

М 3=1

- 2 V;- - вWз,з+1 - 2 V+1

где суммирование по {а} - это суммирование по всем возможным конфигурациям спинов. Оно может быть выражено как:

Ем = Е ■■■ Е ■

о-1=±1,0 оь=±1,0

Введем трансфер-матрицу Т с элементами:

Т(-1, -2) = ехр [^а^ - Ав(-2 + а|)/2].

Тогда мы можем рассмотреть суммирование по спинам а2... ам как умножение матриц Т, а суммирование по спину а1 как взятие следа матрицы. Тогда статсумма может быть записана как:

2 = Тг(Тм)

Зная собственные числа матрицы Т, мы можем написать:

2 = ЛМ + ■■■

В термодинамическом пределе наибольшие собственные числа матрицы Т будут определять свойства системы, так как именно они будут вносить наибольший вклад в 2. Например, свободная энергия системы будет выражаться как:

^ = -Т 1п 2 = -Т 1п(ЛМ + ■■■)

Вычислительная сложность этого алгоритма связана с построением, диагонализа-цией и поиском собственных значений матрицы Т.

1.3.2 Классический метод Монте-Карло

Идея классического метода Монте-Карло (Monte-Carlo, MC), используемая, например, в работах [13, 14] - симулировать интегрирование по фазовому пространству физической системы в состоянии термодинамического равновесия. Это, пожалуй, самый фундаментальный алгоритм вычислительной физики [6].

Пусть мы исследуем среднее значение некоторой термодинамической функции A

(A) = i £ A(s)e-eE <">

s

по всем состояниям системы s. Если вместо того, чтобы пройти по всем состояниям системы (обычно это очень большое число, например, для модели Поттса с q компонен-

Т2\ ~

тами оно составляет qT ), мы предъявим вероятностный процесс Si, который посещает состояния системы с вероятностью, пропорциональной их больцмановскому весу, то для него можно будет записать:

1 ^ ' 'Si)e"eE(si)

(A)« )

2

г

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

P P

1 m 1 n

где Pm и Pn — вероятности того, что система находится в состояниях m и n, wmn — вероятность прямого перехода системы из состояния n в состояние m, а wnm — вероятность обратного перехода системы из состояния m в состояние n. Примером алгоритма, строящего si в соответствии с принципом детального равновесия, может служить алгоритм Метрополиса-Гастингса, который предлагает вероятность перехода между состояниями si и s2 как:

p(s1 ^ s2) = min(1, exp -в(ES1 - Es2)).

Точность метода Монте-Карло оценивается как к 1/VN [6]. Метод в классической интерпретации требует нового расчета для каждого интересного нам значения температуры T = 1/в и параметров системы, что может быть частично нивелировано методами перевзвешивания [6].

1.3.3 Метод отжига популяции

Метод отжига популяции (Population Annealing, PA), описанный, например, в работах [15, 16], представляет собой эволюционный подход, разработанный для преодоления

ограничений классического метода Монте-Карло, особенно в условиях критического замедления. В отличие от традиционного Монте-Карло, где система исследуется в рамках одной траектории, PA использует ансамбль независимых копий (реплик) системы, чтобы эффективно изучать фазовое пространство. Этот подход позволяет избегать застревания в локальных минимумах и улучшает исследование сложных энергетических ландшафтов, характерных для систем с фазовыми переходами.

Рассмотрим вначале большое число R реплик системы при обратной температуре во. На каждом шаге алгоритма популяция состояний S = Si, s2,..., sn эволюционирует через три этапа: понижение температуры T = 1/в, перевзвешивание и релаксация.

• Перевзвешивание: каждое состояние Si с вероятностью, пропорциональной его Больцмановскому весу e-ßE(si), может быть включено в новую выборку (а может и быть отринуто и выброшено из популяции). Это обеспечивает динамическое увеличение числа представителей низкоэнергетических состояний при понижении температуры T.

• Релаксация: после перевзвешивания состояния подвергаются случайным изменениям, например, с помощью алгоритма Метрополиса-Гастингса при температуре T. Это позволяет системам избегать застревания в локальных минимумах и исследовать новые области фазового пространства.

• Понижение температуры T = 1/в, что имитирует физический процесс отжига. Важно отметить, что PA сохраняет детальное равновесие на каждом этапе, обеспечивая корректность статистических оценок. На этом шаге также рассчитываются интересные нам термодинамические параметры.

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

1.3.4 Метод Ванга-Ландау

Гистограммный метод Ванга-Ландау (Wang-Landau, WL), разработанный Фугао Ван-гом и Дэвидом Ландау в 2001 году [17, 18], представляет собой итерационный алгоритм, направленный на прямую оценку плотности состояний ) системы. Этот метод позволяет обойти ограничения классического Монте-Карло и метода отжига популяции, обеспечивая эффективное вычисление термодинамических функций для

широкого диапазона температур в рамках единого эксперимента. Использование алгоритма для исследования модели Блюма-Капеля можно найти, например, в работе

[19].

Идея метода заключается в том, чтобы накапливать гистограмму посещений энергетических уровней E, итеративно корректируя веса переходов между состояниями так, чтобы гистограмма посещений оставалась плоской. Это достигается с помощью модифицированной вероятности перехода p, которая постепенно устраняет предпочтение в пользу высоко- или низкоэнергетических состояний. В результате, оригинальный алгоритм оценивает g(E) с точностью до нескольких процентов [17], что позволяет вычислять любые термодинамические величины через статсумму Z.

Оригинальный алгоритм Ванга-Ландау (WL) выглядит как описано в Алгоритме

1.

Algorithm 1 Ванга-Ландау

0. Инициализация. Введем начальное случайное состояние системы s0, приближение плотности состояний g0(E) = 1 для всех E. Определим гистограмму посещений H(E) = 0 и начальное значение модификатора f0, например, e ~ 2.71828 и конечное значение fend, например, fend = exp(10-8) [17]. while f > fend do

while H(E) не достаточно плоская (относительное отклонение > 5%) do

1. Произведем переход от текущего состояния si к следующему состоянию Sj по правилу Метрополиса-Гастингса с вероятностью Ванга-Ландау:

p(si ^ Sj) = mm(——( F ^, 1)

•in g(Esi) •In g(Es.)

2. Обновим гистограмму посещений: H(ESj) ^ H(ESj) + 1

3. Обновим плотность состояний на фактор f: ing(ESj) ^ ing(ESj) + in f end while

4. Обновляем значение фактора f по некоторому правилу: fk+1 ^ y(fk), например, y(f) = v7, очищаем гистограмму H(E) = 0. end while

Алгоритм оказывается чувствителен к выбору начального фактора f и правилу его изменения у по ходу алгоритма. Например, алгоритм с ) = может не сходиться к верному решению. Оптимальный выбор у оставался открытым до работ [20—22], в которых аргументируется и исследуется выбор у вида у = 1/^ на поздних стадиях вычисления. Также исследуется и ослабление условия для Н для запуска нового цикла: вместо плоскости гистограммы предлагается менее строгое, но все еще достаточное условие: посетить все доступные энергии так, чтобы Н = 0.

Итого, оптимальный алгоритм выглядит следующим образом: до выполнения условия 1п f > NE/£, где NE - число уровней энергии системы, а I - число элементарных переворотов спинов исследуемой реплики, используется оригинальный алгоритм WL.

Однако после этой точки мы меняем правило изменения / на /к+1 = Nе/£ и более не используем гистограмму Н: / обновляется при каждом элементарном перевороте спина. Теоретические основания сходимости этого алгоритма представлены в работе [21].

Именно эту модификацию алгоритма Ванга-Ландау мы станем использовать для сравнения с исследуемым в работе методом МСРА в секции 1.3.4.

Одним из недостатков алгоритма является его неудобность для расчета, например, намагниченности, ведь для этого нам придется строить уже гистограммы Н(Е,т) или д(Е,Х), что значительно удорожит процесс вычислений.

2 Метод MCPA

Представленный в диссертации метод отжига популяции (MCPA) является методом прямой оценки энтропии - это значит, что для систем, характеризующихся неким функционалом (в нашем случае энергией E), результатом работы алгоритма будет функция S (E ) = ln g(E ), с помощью которой можно рассчитать статсумму Z (T ) для любой температуры. Метод также годится для нахождения экстремума такого функционала (E), то есть основного и максимально возбужденного состояния системы. В диссертации рассматриваются двумерные модели Поттса, Изинга и Блюма-Капеля и, минимально, одномерная модель Изинга.

Метод MCPA не оперирует внутри себя температурными вероятностями, а потому, по задумке, не испытывает критического замедления в привычном смысле.

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

2.1 Оригинальный алгоритм Розе-Махты

В своей работе [7] Джон Махта и Натан Розе представили новую технику моделирования равновесных систем в микроканоническом ансамбле с помощью отжига популяции под энергетическим порогом. По построению эта техника предназначена для получения численных значений плотности состояний системы (density of states, DOS, g(E)) и конфигураций системы для заданных наперед значений энергии.

Алгоритм Розе-Махты 2 может выглядеть необычно, поэтому мы рассмотрим его в деталях, с иллюстрациями (См. Рис.1 и Рис.2).

Рассмотрим R реплик интересующей нас, например, спиновой системы, для которой мы знаем шаг спектра (либо сам спектр).

Время релаксации T, несмотря на обманчивое именование, не является настоящим физическим временем - так мы называем некоторым образом нормированное число переворотов спинов. Эта величина стала ассоциироваться с настоящим физическим временем, так как именно эта операция занимает большую часть времени исполнения алгоритма.

Величину T станем изменять в Монте-Карло шагах (Monte-Carlo Steps, MCS): 1 MCS = L2 единичных переворотов спина. Если не оговорено другое, T =10 MCS = 10 L2 единичных переворотов спина. Конкретный выбор числа T зависит от параметров задачи и обсуждается в секции 2.3.

Алгоритм 2 Розе-Махты

0. Инициализация (Initialization): создадим в памяти R реплик системы с максимальной энергией E = Emax. Введем начальный верхний порог энергии U > Emax (иллюстрация на Рис. 1а). while U > Emin do

Список литературы диссертационного исследования кандидат наук Мозоленко Вячеслав Константинович, 2025 год

// // V -

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

E/L2

(e) D = l.8

(f) D = l.9

E/L2 E/L2

(g) D = l.95 (h) D = l.97

- т — X

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

Е/Ь2

(1) Б = 1.98

Ш Б = 1.99

Рис. 35: Вид параметров 5, т и X как функций энергии Е для разных значений Б и Ь = 32.

4.3.7 Общий вид основных термодинамических величин

Для дальнейшей наработки интуиции разберем общий вид основных термодинамических величин, которые мы будем использовать в этой работе: среднюю энергию Е(Т), термодинамическую энтропию 5(Т) = 1п ^ + Е(Т)/Т, вероятность распределения Р(Е,Т), удельную теплоемкость с(Т) = С(Т)Ь—2, Биндеровский кумулянт энергии ВСЕ, намагниченность т(Т), магнитную восприимчивость х(Т), Биндеровский кумулянт намагниченности ВСт, концентрацию X(Т), восприимчивость концентрации У(Т). Все они рассчитаны согласно описанию в секции 2.2.6 и проиллюстрированы на Рис. 36.

(а) Средняя энергия Е(Т)

(Ь) Термодинамическая энтропия 5(Т)

(с) Удельная теплоемкость с(Т)

(а) Плотность вероятности Р(Е, Т) для демонстрационной температуры Т =1.

(е) Биндеровский кумулянт энергии ВСЕ(Т). СредНЯЯ намагниченность т(Е)

(Ь) Биндеровский кумулянт намагниченно-(g) Магнитная восприимчивость х(Т). сти ВСт(Т).

(1) Концентрация дырок X (Т).

П) Дырочная восприимчивость У (Т).

Рис. 36: Вид термодинамических величин Е(Т), 5(Т), т(Т) и X(Т) как функций температуры Т для разных значений Б и Ь = 32.

4.4 Конечномерный анализ для области фазового перехода II рода

Определить (конечномерную) критическую точку Рс можно различными способами. Для области фазового перехода второго рода мы определим и исследуем следующие:

• PC - положение максимума удельной теплоёмкости c

bc

• Pc E - положение минимума Биндеровского кумулянта энергии BCe

• PX - положение максимума удельной магнитной восприимчивости х

• PcCm - положение минимума Биндеровского кумулянта намагниченности BCm

Критическая точка, соответствующая термодинамической величине A, характеризуется температурой T^, значением величины Aс при температуре TA и по умолчанию является функцией размера системы L и параметра модели D.

В области значения параметра —то < D < 0, когда дырочные состояния спинов подавлены, мы ожидаем, что модель Блюма-Капеля будет находится в классе эквивалентности 2D модели Изинга, а в области около трикритической точки Ttr ~ 1.966 в некотором другом классе эквивалентности, который мы будет называть трикрити-ческим. В области значений 0 < D < 1.966 мы ожидаем увидеть смешанный режим. Ожидаемые значения критических показателей, согласно работе [37], представлены в таблице 6.

Класс эквивалентности а в Y V

2D Изинга 0 1/8 7/4 1

Трикритический 8/9 1/24 37/36 5/9

Таблица 6: Ожидаемые значения критических показателей, согласно работе [37]

4.4.1 Удельная теплоемкость

Напомним нотацию: удельная теплоемкость о - это полная теплоемкость О системы, отнесенная к размеру решетки: о(Т) = О(Т)Ь-2.

Как мы уже обсуждали в секциях 1.2.2 и 1.2.3, мы ожидаем, что конечномерная зависимость значения критической удельной теплоемкости ос(Ь) (отоах(Ь)) и соответствующей температуры ТС(Ь) от размера системы Ь выглядит следующим образом:

ос (Ь) = , (48)

ТС (Ь) = ТС (то) + £СЬ-1^. (49)

Однако для случая класса эквивалентности 2Э модели Изинга показатель а = 0, а значит, для удельной теплоемкости нам потребуется следующий порядок точности. Согласно работам [6, 41, 42], в следующем порядке точности мы ожидаем, что

сс (Ь) = со(1п Ь + С1) (50)

Константы не являются универсальными между разными моделями, и все же аналитический расчет констант со и с1 из уравнения (50) для модели Изинга на торе (квадратной решетке с периодическими граничными условиями) может быть найден в работах [41—44], и оказывается, что с0 = 8/п(3 /Тс)2 « 0.494.

Иллюстрация максимума удельной теплоемкости изображена на Рис. 37 и 36с. Для поиска максимума мы используем метод ГштЬои^ библиотеки эс1ру.

0 1 2 3 4 5 6

Т

Рис. 37: Иллюстрация вида зависимости теплоемкости системы для О = 0 и Ь = 40. Красные пунктирная линии обозначают положение максимума теплоемкости (Т,СЬ-2) = (1.704,3.217).

4.4.1.1 Детали аппроксимации удельной теплоемкости Рассмотрим уравнения (48) и (50) и иллюстрацию с соответствующими аппроксимациями на Рис. 38, где мы качественно видим переход от зависимости типа (48) к (50). Если мы хотим количественно оценить этот переход, мы можем исследовать также разброс параметров аппроксимации.

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

диссертации аппроксимации выведены в пределе больших Ь, и чем меньше Ь, тем больше должен быть вклад неучтенных порядков разложения. Первую причину мы оценим с помощью статистической погрешности параметров аппроксимации, второй же требуется больше внимания. Типичный набор размеров решеток в работе: Ь Е [4, 6, 8,10,12,14,16, 20, 24, 28, 32,40]. Мы предлагаем следующую процедуру: построим не одну аппроксимацию, учитывающую все точки, а к (в нашем случае типичное число к = 4), последовательно удаляя точки с наименьшим Ь (процесс проиллюстрирован на Рис. 39). В качестве истинного значения будем брать аппроксимацию с минимальным числом точек (то есть последнюю), а систематический разброс будем оценивать по этим к аппроксимациям. В предположении о независимости ошибок суммарная ошибка некого параметра А будет выглядеть как о\ = а\ + а\ 8у8.

Если не сказано иного, то все погрешности параметров конечномерного анализа далее рассчитываются именно так, с разбивкой на статистическую и систематическую части.

Результаты аппроксимаций и разброс их параметров представлены в таблице 7 и проиллюстрированы на Рис. 40.

(с) В = 1.9 (а) В = 1.95

ю1

L

(е) В = 1.966

Рис. 38: Результаты аппроксимации удельной теплоемкости с помощью уравнений (48) и (50). Ось Ь представлена в лог-масштабе. Легко заметить, что при В < 1.9 логарифмическая аппроксимация (50) работает лучше, чем степенная (48), а после В >= 1.95 ситуация меняется на обратную.

data •

without 0 first points

without 1 first points

— without 2 first points

without 3 first points ^ - ' m

■■■■ ✓ * s • ,,-•* y y y-' s*

101

L

Рис. 39: Иллюстрация к процессу оценки систематической ошибки. Аппроксимации проводятся последовательно, на каждом шаге удаляется все больше точек с наименьшим Ь. В = 1.953.

ce (L) = co (ln L + ci ) ce (L) = AeLa/v

D co ci Ac a/v

-2.0 0.73(1) 0.59(4) 1.14(5) 0.28(1)

-1.5 0.753(9) 0.56(4) 1.14(5) 0.28(1)

-1.0 0.75(2) 0.58(9) 1.16(6) 0.28(2)

0.0 0.78(1) 0.35(6) 1.06(5) 0.30(1)

0.1 0.7б(1) 0.47(б) 1.09(б) 0.29(2)

0.2 0.7б8(8) 0.37(3) 1.0б(б) 0.30(1)

0.3 0.778(7) 0.30(3) 1.02(б) 0.30(1)

0.4 0.7б9(б) 0.34(2) 1.02(б) 0.30(1)

0.б 0.7б(1) 0.29(4) 0.99(4) 0.30(1)

0.б 0.740(б) 0.32(2) 0.99(б) 0.30(1)

0.7 0.722(8) 0.3б(3) 0.98(б) 0.30(1)

0.8 0.72(1) 0.31(4) 0.9б(б) 0.30(1)

0.9 0.б89(4) 0.37(2) 0.9б(б) 0.30(1)

1.0 0.б9б(9) 0.28(4) 0.91(4) 0.30(1)

1.1 0.бб8(8) 0.31(4) 0.89(4) 0.30(1)

1.2 0.ббб(б) 0.2б(3) 0.84(3) 0.31(1)

1.3 0.б3б(7) 0.2б(3) 0.82(3) 0.31(1)

1.4 0.б0(1) 0.33(б) 0.80(3) 0.30(1)

1.б 0.ббб(7) 0.39(4) 0.77(3) 0.29(1)

1.б 0.б4(1) 0.3б(9) 0.72(2) 0.299(9)

1.7 0.48(1) 0.б(1) 0.72(1) 0.283(б)

1.7б 0.47(1) 0.б(1) 0.717(8) 0.280(3)

1.8 0.41(1) 1.0(1) 0.77(1) 0.2б1(б)

1.8б 0.37(1) 1.б(1) 0.88(2) 0.22(1)

1.87 0.3б(1) 1.9(2) 0.94(4) 0.21(1)

1.9 0.41(2) 1.б(2) 0.98(б) 0.22(2)

1.91 0.47(2) 1.2(2) 0.97(б) 0.24(2)

1.92 0.б2(3) 1.0(2) 1.00(8) 0.2б(3)

1.93 0.71(2) 0.24(9) 0.90(7) 0.31(2)

1.94 1.01(1) -0.47(3) 0.77(б) 0.39(2)

1.9б 1.б4(9) -1.2(1) 0.бб(2) 0.бб(1)

1.9б1 1.7(1) -1.3(1) 0.б2(2) 0.б8(1)

1.9б2 1.8(1) -1.3(1) 0.49(2) 0.б0(1)

1.9б3 2.0(1) -1.4(1) 0.44(1) 0.б4(1)

1.9б4 2.1(1) -1.б(1) 0.42(1) 0.бб4(8)

1.9бб 2.2(2) -1.б(1) 0.390(7)

1.9бб 2.4(2) -1.б(1) O^l^) 0.738(б)

1.9б7 2.б(2) -1.б(1) 0.332(б) 0.7б3(б)

1.9б8 2.7(3) -1.7(2) 0.30б(б) 0.797(б)

1.9б9 2.9(3) -1.8(2) 0.272(8) 0.840(9)

1.9б 3.2(3) -1.8(2) 0.24(1) 0.89(1)

1.9б1 3.4(4) -1.9(2) 0.22(1) 0.92(1)

1.9б2 3.б(б) -1.9(2) 0.20(1) 0.97(1)

1.9б3 3.9(б) -2.0(2) 0.17(1) 1.02(2)

1.9б4 4.2(б) -2.0(2) 0.14(1) 1.09(3)

1.965 4.5(6) -2.0(2) 0.13(1) 1.13(3)

1.966 4.8(7) -2.1(2) 0.11(1) 1.17(4)

1.967 5.1(8) -2.1(2) 0.10(1) 1.24(4)

1.968 5.3(9) -2.1(2) 0.08(1) 1.28(4)

1.969 5.7(10) -2.1(2) 0.08(1) 1.32(5)

1.97 5.8(10) -2.2(2) 0.07(1) 1.37(5)

1.971 6.1(11) -2.2(2) 0.06(1) 1.39(6)

1.972 6.1(11) -2.2(2) 0.05(1) 1.44(6)

1.973 6.2(11) -2.2(2) 0.05(1) 1.46(6)

1.974 6.4(12) -2.2(2) 0.05(1) 1.47(6)

1.975 6.3(12) -2.2(2) 0.04(1) 1.50(6)

1.976 6.3(11) -2.2(2) 0.04(1) 1.50(7)

1.977 6.3(12) -2.2(2) 0.05(1) 1.48(7)

1.978 6.2(11) -2.2(2) 0.05(1) 1.47(7)

1.979 6.0(11) -2.2(2) 0.04(1) 1.48(7)

1.98 5.7(11) -2.2(2) 0.04(1) 1.49(7)

1.981 5.5(10) -2.2(2) 0.04(1) 1.47(7)

1.982 4.1(8) -2.1(2) 0.07(2) 1.30(8)

1.983 3.9(7) -2.0(2) 0.07(2) 1.27(8)

1.984 3.6(6) -2.0(2) 0.07(2) 1.27(9)

1.985 2.8(5) -1.9(3) 0.10(4) 1.1(1)

1.986 2.6(5) -1.9(3) 0.10(4) 1.1(1)

1.987 2.4(5) -1.8(3) 0.10(5) 1.1(1)

1.988 1.2(2) -1.2(4) 0.28(8) 0.70(8)

1.989 1.1(2) -1.1(4) 0.30(8) 0.65(8)

1.99 0.9(1) -0.8(5) 0.3(1) 0.58(9)

1.995 0.01(4) - 1.5(1) -0.09(4)

2.0 0.02(2) - 1.20(9) -0.03(2)

Таблица 7: Результаты конечномерного анализа удельной теплоемкости сс с помощью аппроксимаций (48) и (50), иллюстрированные на Рис. 40

(я) с0 для всех исследуемых значений Б. (Ь) с0 в окрестности трикритической точки.

(с) сх для всех исследуемых значений Б. (^ сх в окрестности трикритической точки.

(е) Ас для всех исследуемых значений Б. (Г) Ас в окрестности трикритической точки.

а/у для всех исследуемых значений Б. (Ь) а/у в окрестности трикритической точки.

Рис. 40: Иллюстрация полученных значений параметров с0 и с1 из аппроксимации уравнением (50) и Ас и а/у из аппроксимации уравнением (48). Зеленая пунктирная линия вверху обозначает ожидаемое значение а/у = 8/5, а большой красной звездочкой обозначено значение величины для Б = 1.966.

На Рис. 40 мы можем наблюдать несколько занимательных тенденций. Во-первых, логарифмическая аппроксимация (50) обладает отличной точностью при Б < 1.95 и начинает ухудшаться при приближении к трикритической точке. Во-вторых, степенная аппроксимация ведет себя наоборот: улучшается при приближении к трикритической точке, как и ожидалось. В-третьих, значение отношения критических показателей а/у оказывается значительно ниже ожидаемого значения 8/5.

4.4.1.2 Детали аппроксимации критической температуры Рассмотрим внимательнее уравнение (49) и иллюстрацию на Рис. 41. Как проиллюстрировано, когда

мы пытаемся приблизить данные функцией из уравнения (49) с двумя свободными переменными V и ТС(то), мы получаем довольно абсурдные результаты. Это распространенная ситуация в задачах оптимизации, и, чтобы ее решить, у нас есть несколько способов.

Во-первых, мы можем зафиксировать Тс, но нужно его откуда-то взять, а аналитического решения модели на момент написания работы не известно. Во-вторых, мы можем зафиксировать V - его ожидаемые аналитические значения нам как раз известны. В-третьих, мы могли бы попытаться ввести некоторую регуляризацию на значения параметров. Однако этот путь не будет иметь физического смысла и будет значительно зависеть от выбранной регуляризации, то есть от желания исследователя. Поэтому мы остановимся на втором варианте: зафиксируем V =1 и V = 5/9 и исследуем ошибки полученных аппроксимаций. Соответствующие упрощенные выражения для температуры приведены в уравнениях (51) и (52). Результаты этого упражнения представлены в таблице и проиллюстрированы на Рис. 42.

ТС (Ь) = ТС М + БСЬ-1 (51)

ТС (Ь) = ТС (то) + БсЬ-9/ъ (52)

1.68 1.66

0.00 0.05 0.10 0.15 0.20 0.25

1/Ь

Рис. 41: Иллюстрация для аппроксимации температуры критической точки удельной теплоемкости. Зеленой пунктирной линией изображена аппроксимация согласно уравнению (49) с неправдоподобными оптимальными параметрами V = -147, БС = 148 и V = 20188. Сплошная оранжевая линия - аппроксимация с помощью упрощенной формулы (51).

Tc = Tc (то) + BcL-1 tc = tc (то) + BcL-9/5

D tc (то) Bc Tc (то) Bc

-2.0 2.001(1) 0.б3(4) 2.011(1) 3.4(5)

-1.б 1.947(1) 0.4б(3) 1.956(1) 2.9(4)

-1.0 1.879(2) 0.43(б) 1.887(1) 2.7(5)

0.0 1.б9б(1) 0.33(3) 1.702(1) 2.0(4)

0.1 1.б74(1) 0.30(4) 1.6793(9) 1.9(4)

0.2 1.бб1б(б) 0.24(2) 1.6559(6) 1.5(2)

0.3 1.б21(1) 0.32(4) 1.627(1) 2.0(4)

0.4 1.б98б(9) 0.22(3) 1.6025(6) 1.4(2)

0.б 1.бб97(9) 0.23(2) 1.5738(7) 1.5(2)

0.б 1.б44(1) 0.14(4) 1.5467(8) 0.9(3)

0.7 1.б07(1) 0.2б(3) 1.5112(9) 1.5(3)

0.8 1.4778(9) 0.14(2) 1.4802(5) 0.9(2)

0.9 1.43б(1) 0.2б(б) 1.441(1) 1.6(4)

1.0 1.402(1) 0.1б(4) 1.4047(8) 1.0(3)

1.1 1.3б19(9) 0.11(2) 1.3640(6) 0.6(2)

1.2 1.319(1) 0.0б(3) 1.3198(6) 0.4(2)

1.3 1.2б8(1) 0.0б(3) 1.2697(6) 0.4(2)

1.4 1.21б(1) 0.04(б) 1.2166(9) 0.1(3)

1.б 1.1б7(2) 0.00(б) 1.1568(9) -0.1(3)

1.б 1.089(2) -0.0б(б) 1.0880(8) -0.6(3)

1.7 1.012(1) -0.23(б) 1.0082(5) -1.6(2)

1.7б 0.9бб(2) -0.28(б) 0.9605(6) -2.0(2)

1.8 0.91б(1) -0.43(б) 0.9083(4) -3.0(1)

1.8б 0.8б7(1) -0.б9(4) 0.8472(7) -4.0(3)

1.87 0.8279(б) -0.б3(2) 0.8169(9) -4.1(3)

1.9 0.773б(8) -0.б1(2) 0.7647(7) -3.4(2)

1.91 0.7ббб(4) -0.б1(1) 0.7466(7) -3.3(3)

1.92 0.73б(1) -0.48(3) 0.7275(9) -3.1(4)

1.93 0.712б(б) -0.39(1) 0.7059(5) -2.6(2)

1.94 0.б880(4) -0.31(1) 0.6826(4) -2.1(1)

1.9б 0.бб11(4) -0.22(1) 0.6573(2) -1.5(1)

1.9б1 0.бб84(4) -0.21(1) 0.6548(2) -1.4(1)

1.9б2 0^2(2) -0.220(7) 0.6524(3) -1.4(1)

1.9б3 0.бб29(4) -0.20(1) 0.6495(2) -1.3(1)

1.9б4 0.бб01(3) -0.193(9) 0.6468(2) -1.3(1)

1.9бб 0.б4б7(4) -0.17(1) 0.6437(2) -1.2(1)

1.9бб 0.б439(3) -0.1б8(9) 0.6410(2) -1.12(9)

1.9б7 0.б410(3) -0.1б(1) 0.6382(2) -1.08(9)

1.9б8 0.б379(3) -0.1б0(9) 0.6353(1) -1.01(8)

1.959 0.6343(3) -0.134(9) 0.6320(1) -0.90(7)

1.96 0.6307(4) -0.11(1) 0.6288(1) -0.79(7)

1.961 0.6279(3) -0.12(1) 0.6259(1) -0.77(7)

1.962 0.6244(2) -0.101(7) 0.6226(1) -0.67(5)

1.963 0.6208(4) -0.09(1) 0.6194(1) -0.58(6)

1.964 0.6174(2) -0.077(8) 0.6161(1) -0.53(5)

1.965 0.6139(2) -0.069(7) 0.6128(1) -0.48(5)

1.966 0.6099(2) -0.051(7) 0.6091(1) -0.34(4)

1.967 0.6060(3) -0.03(1) 0.6055(1) -0.24(5)

1.968 0.6020(2) -0.018(6) 0.6017(1) -0.13(4)

1.969 0.5980(2) -0.002(7) 0.5980(1) -0.03(4)

1.97 0.5939(1) 0.015(5) 0.5942(1) 0.08(3)

1.971 0.5897(1) 0.026(5) 0.59017(9) 0.15(3)

1.972 0.5852(2) 0.046(6) 0.5860(1) 0.28(5)

1.973 0.5805(1) 0.063(3) 0.58164(9) 0.41(3)

1.974 0.5761(1) 0.077(5) 0.57743(9) 0.49(3)

1.975 0.5715(3) 0.09(1) 0.5731(1) 0.59(5)

1.976 0.5658(1) 0.127(5) 0.5680(1) 0.82(7)

1.977 0.5608(2) 0.145(8) 0.5634(1) 0.94(6)

1.978 0.5564(4) 0.14(1) 0.5589(1) 0.95(6)

1.979 0.5487(4) 0.22(1) 0.5526(3) 1.4(1)

1.98 0.5413(7) 0.28(2) 0.5462(6) 1.8(2)

1.981 0.5405(9) 0.19(3) 0.5435(2) 1.3(1)

1.982 0.531(1) 0.28(3) 0.5368(6) 1.7(1)

1.983 0.527(3) 0.26(8) 0.532(1) 1.6(3)

1.984 0.521(1) 0.31(3) 0.5270(5) 1.9(1)

1.985 0.502(1) 0.52(2) 0.514(1) 2.9(3)

1.986 0.509(5) 0.3(1) 0.516(2) 2.0(4)

1.987 0.498(4) 0.43(9) 0.507(1) 2.5(4)

1.988 0.458(3) 0.92(6) 0.485(4) 4.4(6)

1.989 0.460(5) 0.85(8) 0.485(4) 4.1(6)

1.99 0.449(6) 0.9(1) 0.475(4) 4.5(7)

1.995 0.60(2) -1.4(7) 0.58(1) -8.3(46)

2.0 0.548(5) 0.3(1) 0.555(3) 1.9(11)

Таблица 8: Результаты конечномерного анализа критической температуры Т^ согласно уравнениям (51) и (52), иллюстрированные на Рис. 42

(Ь) ТС (го) в окрестности трикритической (а) ТС (го) для всех исследуемых значений Б. точки.

(с) ВС для всех исследуемых значений Б. (^ ВС в окрестности трикритической точки.

Рис. 42: Иллюстрация полученных значений параметров Т^(го) и Вс из аппроксимации уравнениями (51) и (52). Большая красная и желтая звездочки соответственно обозначают значения величин для Б = 1.966.

Из Рис. 42 можно наблюдать, что температуры, полученные с помощью двух аппроксимаций, довольно хорошо совпадают. Погрешности для аппроксимации (51) оказываются исключительно маленькими во всей области исследуемых Б, а погрешности аппроксимации (52), будучи значительно больше, уменьшаются в окрестности трикритической точки. И все же сделать итоговый вывод о классах эквивалентности на примере критической температуры удельной теплоемкости не выходит: критический показатель V в целом оказывается довольно сложно оценить [6].

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

4.4.2 Биндеровский кумулянт энергии

Биндеровский кумулянт энергии (БКЭ) рассчитывается согласно уравнению (25) и в нем проявляется несколько сложных эффектов. Общий вид БКЭ уже изображен на Рис. 36е, однако нам понадобится взглянуть на него поближе на Рис. 43. Как мы можем видеть, положение минимума БКЭ совершенно не совпадает с ожидаемым положением критической точки, в котором, однако, можно наблюдать локальный минимум. При разных значениях Б он может находится по разные стороны от большого минимума.

(а) В = 0.5 (Ь) В = 1.95

Рис. 43: Иллюстрация поведения БКЭ для разных значений В. Зеленой пунктирной линией изображено ожидаемое положение критической точки (согласно таблице 8), красной точко-пунктирной линией изображено положение "фальшивого" (большого) минимума БКЭ.

Для того, чтобы разобраться в природе происходящего, напомним, как мы считаем БКЭ:

(Е 4)

ВСЕ (Т) = 1 - о/ .

3(Е

Рассмотрим распределение вероятности Р(Е, Т) и величины ЕР(Е, Т) при температуре фальшивого пика на рис. 44.

(а) Р(Е, Т) для параметров В = 0.5, Т = 2.8 (Ь) Р(Е, Т) для параметров В = 1.95, Т = 0.5

(с) ЕР(Е, Т), В = 0.5, Т = 2.8 (а) ЕР(Е, Т), В = 1.95, Т = 0.5

Рис. 44: Иллюстрация поведения плотности вероятности Р(Е, Т) и величины ЕР(Е, Т) для разных значений В при температуре Т, соответствующей фальшивому пику.

Мы предполагаем, что дело в пересечении плотностью вероятности линии E = 0, а фальшивый пик - результат математического артефакта. Избавиться от этого артефакта можно несколькими разными способами. Например, мы можем разделить энергию на части: Изинга и Блюма-Капеля:

H = - ^ aaj +D ^ а2 = Ei + EB, (53)

(ij ^ V J у

Ising Part Blume-Capel Part

и построим соответствующие БКЭ:

E)

BCi (T) = 1 -BCb (T) = 1 -

3<Ej)

(E4B) 3<EB)

Заметим, что усреднение (*) = ^Е (Е)е-Е/т - это все еще усреднение по Гиббсовско-му распределению с полной энергией, НЕ по Гиббсовскому распределению с частями энергии Ехв. Все эти способы дадут довольно похожие результаты, проиллюстрированные на Рис. 45.

2

2

(а) 2/3 - ВС, В = 0.5, Т = 2.8 (Ь) 2/3 - ВС, В = 1.95, Т = 0.5

Рис. 45: Иллюстрация поведения плотности вероятности Р(Е, Т) и величины ЕР(Е, Т) для разных значений В при температуре Т, соответствующей фальшивому пику. Для наглядности вычтен фактор 2/3, домножено на —1 и представлено в лог-масштабе.

Исследуем именно пики ВС/ и ВСВ и их конечномерное поведение. Рассмотрим их типовой вид на Рис. 46. Как обсуждалось в секции 1.2.3, мы ожидаем увидеть конечномерную зависимость для критической температуры БКЭ вида:

ТВ с (Ь) = ТВ с (то) + Ввс, (54)

но по причинам, описанным в секции 4.4.1.2, на практике мы будем использовать упрощенную версию зависимости с V =1:

ТВс (Ь) = ТВс (то) + Ввс/Ь. (55)

На Рис. 46 мы наблюдаем, что выбранная аппроксимация хорошо ложится на данные моделирования.

.........

[ • НС п. с

■ ВС1:С

--- ВС^сЫ+В^Ь

1 /ь 1/ь

(а) Значения критических ВС, В = 1.5 (Ь) Температуры критических Твс, В = 1.5 Рис. 46: Иллюстрация для конечномерного анализа БКЭ для частей энергии.

Результаты такого конечномерного анализа представлены в таблице 9 и проиллюстрированы на Рис. 47.

В ТВ св (то) Вв Тв с' (то) В/

1.0 1.394(1) 0.625(1) 1.380(7) 1.7(2)

1.1 1.355(2) 0.573(2) 1.320(8) 2.1(1)

1.3 1.262(1) 0.535(1) 1.23(3) 1.9(7)

1.5 1.150(1) 0.480(1) 1.13(2) 1.7(5)

1.9 0.764(1) 0.327(1) 0.75(1) 0.9(2)

1.95 0.654(3) 0.376(3) 0.65(1) 0.7(2)

1.96 0.623(5) 0.444(5) 0.62(1) 0.8(3)

1.97 0.581(4) 0.715(4) 0.57(2) 1.1(3)

1.98 0.537(3) 0.842(3) 0.53(1) 1.1(2)

1.99 0.476(3) 1.193(3) 0.48(1) 1.5(2)

Таблица 9: Результаты конечномерного анализа критических температур БКЭ. Малое число различных О связано с двумя факторами: минимум БК не всегда существует (из-за аномалии) и расчет частей энергии требует дополнительных прогонов алгоритма МСРЛ в условиях ограниченного времени.

Рис. 47: Результаты конечномерного анализа БКЭ для частей энергии.

4.4.3 Магнитная восприимчивость

Как мы уже обсуждали в секциях 1.2.2 и 1.2.3, мы ожидаем, что конечномерная зависимость значения критической магнитной восприимчивости хс(Ь) и соответствующей температуры Т£(Ь) от размера системы Ь выглядит следующим образом:

Хс (Ь) = , (56)

ТЭД = Т*(<х>) + Вх/Ь^. (57)

Однако, как мы описывали в секции 4.4.1.2, нам приходится упростить асимптотику температуры до:

Т*(Ь) = Т*М + Вх/Ь. (58)

В отличие от ситуации с удельной теплоемкостью в секции 4.4.1, мы ожидаем, что 7 будет меняться от 7/4 вдали от трикритической точки до 37/36 вблизи трикрити-ческой точки, и нам не понадобятся логарифмические поправки к скейлингу и чересчур сложный анализ. Результаты аппроксимации с помощью уравнений (56) и (58) представлены в таблице 10 и проиллюстрированы на Рис. 48. Мы можем наблюдать замечательное соответствие полученного значения 7/^ ожидаемому аналитическому значению.

Хс (Ь) = ЛХЬ Т*(Ь) = Т*М + Вх/Ь

О Лх 1/у т*м Вх

-2.0 0.0466(6) 1.770(4) 2.0019(7) 1.46(1)

-1.5 0.0471(6) 1.769(4) 1.9461(8) 1.40(1)

-1.0 0.0478(4) 1.768(3) 1.879(1) 1.34(1)

0.0 0.0502(1) 1 . 763(1) 1.696(1) 1.18(1)

0.1 0.0507(2) 1.761(1) 1.672(1) 1.170(9)

0.2 0.0506(4) 1.764(2) 1.648(1) 1.158(9)

0.3 0.0512(4) 1.762(2) 1.622(1) 1.153(9)

0.4 0.0518(4) 1.760(3) 1.596(1) 1.13(1)

0.5 0.0522(6) 1.759(4) 1.567(2) 1.13(1)

0.6 0.0526(8) 1.760(5) 1.538(1) 1.09(1)

0.7 0.054(1) 1.756(6) 1.506(1) 1.07(1)

0.8 0.054(1) 1.756(6) 1.473(2) 1.05(1)

0.9 0.056(1) 1.750(7) 1.438(2) 1.02(1)

1.0 0.056(1) 1.75(1) 1.400(2) 0.99(1)

1.1 0.057(1) 1.75(1) 1.359(2) 0.95(1)

1.2 0.059(2) 1.74(1) 1.315(2) 0.93(1)

1.3 0.061(3) 1.74(1) 1.266(3) 0.90(1)

1.4 0.064(4) 1.73(2) 1.215(3) 0.84(2)

1.5 0.067(5) 1.73(2) 1.156(4) 0.79(2)

1.6 0.073(7) 1.71(2) 1.088(5) 0.73(3)

1.7 0.08(1) 1.68(3) 1.010(5) 0.62(3)

1.75 0.09(1) 1.67(4) 0.964(6) 0.55(4)

1.8 0.11(2) 1.63(5) 0.914(6) 0.43(4)

1.85 0.14(2) 1.58(5) 0.854(6) 0.28(4)

1.87 0.17(3) 1.54(5) 0.827(5) 0.19(4)

1.9 0.24(2) 1.47(3) 0.778(5) 0.09(3)

1.91 0.29(2) 1.44(2) 0.760(4) 0.03(3)

1.92 0.34(1) 1.41(1) 0.740(3) -0.01(2)

1.93 0.401(5) 1.394(4) 0.717(3) -0.02(2)

1.94 0.46(1) 1.40(1) 0.691(3) -0.02(1)

1.95 0.48(2) 1.48(1) 0.661(3) 0.01(1)

1.951 0.47(2) 1.49(1) 0.658(3) 0.02(1)

1.952 0.46(2) 1.51(1) 0.655(3) 0.01(1)

1.953 0.45(2) 1.53(1) 0.652(3) 0.02(1)

1.954 0.45(2) 1.54(1) 0.649(3) 0.02(1)

1.955 0.44(1) 1.56(1) 0.646(3) 0.02(1)

1.956 0.42(1) 1.59(1) 0.643(3) 0.03(2)

1.957 0.41(1) 1.61(1) 0.640(3) 0.03(1)

1.958 0.40(1) 1.63(1) 0.636(3) 0.04(1)

1.959 0.38(1) 1.66(1) 0.633(3) 0.04(1)

1.96 0.36(1) 1.69(1) 0.629(3) 0.05(2)

1.961 0.34(1) 1.73(1) 0.626(3) 0.05(2)

1.962 0.33(1) 1.76(1) 0.622(3) 0.06(2)

1.963 0.31(2) 1.79(2) 0.619(3) 0.05(2)

1.964 0.29(2) 1.83(2) 0.616(3) 0.06(2)

1.965 0.28(2) 1.86(2) 0.612(3) 0.05(1)

1.966 0.26(2) 1.90(3) 0.608(4) 0.07(2)

1.967 0.24(2) 1.94(3) 0.605(3) 0.06(2)

1.968 0.23(2) 1.97(3) 0.600(3) 0.08(2)

1.969 0.22(3) 2.00(3) 0.597(3) 0.08(1)

1.97 0.21(3) 2.03(4) 0.593(3) 0.08(2)

1.971 0.21(3) 2.05(4) 0.589(3) 0.08(2)

1.972 0.20(3) 2.07(4) 0.585(4) 0.09(2)

1.973 0.20(3) 2.09(4) 0.580(4) 0.10(2)

1.974 0.20(3) 2.10(4) 0.576(4) 0.10(2)

1.975 0.21(3) 2.11(4) 0.571(4) 0.10(2)

1.976 0.21(3) 2.11(4) 0.567(4) 0.11(2)

1.977 0.21(3) 2.13(5) 0.562(5) 0.12(2)

1.978 0.22(4) 2.13(5) 0.557(5) 0.11(2)

1.979 0.22(5) 2.13(6) 0.552(6) 0.13(3)

1.98 0.22(5) 2.14(6) 0.546(7) 0.15(3)

1.981 0.23(7) 2.15(8) 0.540(8) 0.16(3)

1.982 0.24(7) 2.14(8) 0.535(8) 0.16(3)

1.983 0.2(1) 2.2(1) 0.53(1) 0.14(4)

1.984 0.3(1) 2.1(1) 0.52(1) 0.20(4)

1.985 0.2(1) 2.2(1) 0.52(1) 0.18(4)

1.986 0.3(1) 2.1(1) 0.519(8) 0.05(4)

1.987 0.3(1) 2.1(1) 0.51(1) 0.13(4)

1.988 0.3(2) 2.2(1) 0.49(1) 0.31(7)

1.989 0.3(2) 2.1(1) 0.49(1) 0.20(6)

1.99 0.3(3) 2.1(2) 0.48(2) 0.15(9)

1.995 0.3(16) 2.2(4) 0.41(8) 0.5(3)

2.0 15.01(1108) 1.2(4) - -

Таблица 10: Результаты конечномерного анализа магнитной восприимчивости х.

ф

— 7/4 1

117/36 1

"679" 1

★

■

1

2.1 2.0 1.9

7/1/ 1.8 1.7 1.6 1.5 1.4

-- 7/4

"579" | ★

1

I

2 I

1 I 1 ■ , , I

1.90