Письма в ЖЭТФ, том 110, вып. 10, с. 700 - 705

К численному расчету фрустраций в модели Изинга

А. Г. Макаров^+∗1), К. В. Макарова^+∗, Ю. А. Шевченко^+∗, П. Д. Андрющенко^+∗, В. Ю. Капитан^+∗,

К.С.Солдатов^+∗, А.В.Пержу^+∗, А.Е.Рыбин^+∗, Д.Ю.Капитан^+∗, Е.В.Васильев^+∗, Р.А.Волотовский^+∗

Ю. В. Чубов⁺, К. В. Нефедев^+∗

⁺Школа естественных наук, Дальневосточный федеральный университет, 690091 Владивосток, Россия

^∗Институт прикладной математики, Дальневосточное отделение РАН, 690041 Владивосток, Россия

Поступила в редакцию 22 сентября 2019 г.

После переработки 16 октября 2019 г.

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

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

носительное число возбужденных парных взаимодействий в гамильтониане системы. Нами был модер-

низирован алгоритм Метрополиса до гибридного мультиспинового Монте-Карло метода, использующего

квазимарковскую цепь случайных событий. Объединение канонического и мультиканонического семпли-

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

висимости приращения числа возбуждений и приращения энтропии гексагональной решетки искусствен-

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

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

гексагональной решетки точечных диполей. Предложенная нами методика расчета количественной ме-

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

DOI: 10.1134/S0370274X19220120

В работе [1] введена простая эмпирическая ме-

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

ра фрустрации, которая впоследствии стала широко

ми фрустраций магнитный порядок и макроскопиче-

применяться [2-6]. Необходимо отметить два суще-

ское вырождение нижайшей энергии отсутствовали.

ственных момента, которые нужно учитывать при

Ванье показал [7], что антиферромагнитная модель

определении температуры Кюри-Вейса θ_CW . Пер-

Изинга на треугольной решетке имеет очень большое

вый состоит в том, что экспериментатор может

вырождение основных состояний. Иногда такое вы-

работать только в ограниченном интервале конеч-

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

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

деление, характеристику фрустраций [10].

клонной линейной асимптоты y

= aT + b функ-

Фрустрации порой определяются как конкурен-

ции 1/χ(T ) необходимо одновременное существова-

ция между взаимодействиями, которая приводит к

ние двух конечных пределов lim (χ(T )T )^-1 = a

T →+∞

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

(

)

и lim

χ^-1(T) - aT

= b. Если хотя бы один из

[11, 10]. Впервые термин “фрустрация” был применен

T →+∞

к магнетикам Джерардом Тулузом в 1977 г. в рабо-

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

тах [12, 13]. Энергия взаимодействия пары магнит-

ствует. Второй момент состоит в том, что экспери-

ных моментов в общем случае зависит от расстояния

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

между ними. Суммарная энергия всех парных взаи-

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

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

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

или конфигурацию системы. Для

Ĥ (2, 2) векторных

ра напряженности магнитного поля. Магнитное поле

моделей, см. классификацию в работе [14], значение

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

энергии каждой пары взаимодействующих моментов

го влияния на фрустрированные пары, суммарный

фиксировано, может измениться лишь знак в зависи-

магнитный момент которых равен нулю.

мости от их взаимной ориентации. В модели Изинга

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

пара магнитных моментов является фрустрирован-

чены точные решения модели антиферромагнетика

ной, имеет неудовлетворенную связь, если энергия

Изинга на треугольной [7, 8], кагоме [9] решетках

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

ким образом, фрустрация в этой векторной модели

¹⁾e-mail: makarov.ag@dvfu.ru

есть элементарное возбуждение.

700

Письма в ЖЭТФ том 110 вып. 9 - 10

2019

К численному расчету фрустраций в модели Изинга

701

Особыми микросостояниями или конфигурация-

ный мультиспиновый метод, который является логи-

ми являются основные состояния. Они соответству-

ческим продолжением алгоритма Метрополиса. Рас-

ют значению энергии E_gs относительного миниму-

смотрим случай ферромагнитной J > 0 системы, со-

ма (минимум по отношению ко всем существующим

стоящей из N спинов Изинга S_i = ±1, расположен-

разрешенным значениям энергии). В конфигурации

ных на квадратной решетке, и гамильтонианом

основного состояния число возбуждений или фруст-

∑

H = -J S_iS_j,

(1)

рации является минимальным. Значение суммарной

〈i,j〉

энергии взаимодействия E_min системы магнитных

моментов будет абсолютным минимумом, если все

где суммирование 〈i, j〉 проводится по уникальным

вклады будут отрицательными. Абсолютный мини-

парам ближайших четырех соседей z = 4.

мум энергии для фрустрированных систем лежит

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

в запрещенной зоне, т.е. не существует конфигура-

фигурации на шаге t в односпиновом алгоритме

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

Метрополиса-Гастингса [16] определяется как

ствия, а основное состояние в этом случае является

(

)

возбужденным. Для вычисления количественной ме-

p(E_{t})

p(E_{t} → E_{t+1}) = min

ры фрустраций необходимо знать значения энергии

p(E_{t+1})

(

абсолютного минимума и основного состояния.

[E_{t} -E_{t+1}])

= min

1, exp

(2)

Решение задачи поиска основных состояний и ис-

k_BT

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

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

Необходимо отметить, что формула (2) верна для

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

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

численных расчетов. Эта задача может оказаться

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

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

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

взгляд [15]. Нахождение алгоритма поиска основно-

ведем Монте-Карло семплирование пространства со-

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

го состояния является определяющей проблемой в

теории фрустрированного магнетизма. Основное со-

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

численное с помощью формулы (2) идентичны при

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

ствием ближайших соседей на решетках специфиче-

бесконечном числе МК-шагов для любой температу-

ских геометрий бросает вызов всем теоретическим

ры.

Выделим область (ядро) размером n = 2 × 2 от-

и численным методам, используемым до настоящего

времени.

носительно большой решетки с периодическими гра-

ничными условиями, красные стрелочки на рис. 1.

Существующие сегодня МК- (Монте-Карло) ме-

Взаимодействие kernel-спинов (или k-спинов) с си-

тоды [16-20] имеют свои тонкие места. Односпи-

стемой происходит через граничные border-спины

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

(или b-спины), синие стрелочки на рис. 1. В границу

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

включены b = 8 спинов, с которыми могут напря-

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

мую взаимодействовать спины ядра. Для каждой из

чете термодинамики систем относительно большого

2^b конфигураций граничных спинов рассчитываются

размера. Использование односпиновых МК-методов

все возможные конфигурации ядра (2ⁿ). Для указан-

(например, алгоритм Метрополиса) для расчета ос-

ной на рис. 1 подсистемы внутреннюю энергию кон-

новного состояния систем с грубым ландшафтом

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

энергии сопряжено с трудноразрешимыми проблема-

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

ми. Для преодоления больших энергетических барье-

ров, разделяющих квазизрожденные конфигурации

E_tot = E⁰ + E^∗,

(3)

фрустрированного магнетика Изинга, которые меша-





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

(S₄S₅ + S₄S₈ + S₅S₉ + S₈S₉ ) +

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



E^∗ = -J ( S₁S₄ + S₂S₅ + S₃S₄ + S₇S₈ +

, (4)

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

S₆S₅ + S₁₀S₉ + S₁₁S₈ + S₁₂S₉ )

тиспиновых кластеров.

Для решения задачи о термодинамике фрустри-

где в E^∗ включены только те парные взаимодей-

рованных векторных моделей сложных систем мно-

ствия, в которых есть хотя бы один k-спин. Осталь-

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

ные пары выделены в константную часть E⁰ и пока

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

не рассматриваются. Первые 4 парных произведения

Письма в ЖЭТФ том 110 вып. 9 - 10

2019

702

А. Г. Макаров, К. В. Макарова, Ю. А. Шевченко, П. Д. Андрющенко и др.

общем получаем 2^n+b всевозможных конфигураций

ядра и граничных спинов.

Перед началом численного расчета для каждой

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

храняется локальная плотность состояний ядра фик-

сированного размера. Для каждого t цикла алгорит-

ма Метрополиса:

1. Выбирается случайная область системы - ядро.

2. Извлекается информация о плотности состо-

яний для заданной конфигурации граничных

спинов с ядром, выбранным на шаге 1. Вычис-

ляется Z_k и p(E^∗i) для всех i.

3. Вычисляются средние термодинамические ха-

рактеристики для текущего МК-шага, опреде-

ляя их по формуле полной вероятности. Напри-

Рис. 1. (Цветной онлайн) Блок для прохода по системе

мер, средняя энергия системы:

двумерной модели Изинга. Красным выделены спины,

входящие в ядро (k-спины), синим - спины, гранича-

∑

щие с ядром (b-спины)

〈E〉_t = E⁰ + p(E^∗i)E^∗i.

(8)

в E^∗ представляют парные энергии взаимодействия

4. Итоговая конфигурация ядра выбирается с ве-

внутри выделенной области, и оставшиеся 8 - взаи-

роятностью p(E^∗i).

модействие подсистемы с границей.

Вероятность появления энергии E^∗i в ядре задана

По завершении q циклов термодинамически усред-

распределением Гиббса и определяется как:

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

[

]

ожидание:

E^∗i

p(E^∗i) =

exp -

(5)

∑

Z_k

k_BT

〈E〉 =

E_t.

(9)

где

[

]

∑

E^∗j

Температурные зависимости теплоемкости C(T )

Z_k =

exp -

(6)

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

k_BT

j=1

ной восприимчивости χ(T ) при H → 0, см. рис. 2,

представляет статистическую сумму ядра.

рассчитывались на один спин Изинга

Подставив уравнения (3) и (6) в (5), получим:

1 〈E²〉 - 〈E〉

[

]

[

]

C(T ) =

(10)

E^∗i

∑

E^∗j

N k_BT²

p(E_i) = exp -

exp -

(7)

k_BT

j=1

〈M²〉 - 〈M〉²

χ(T )|_H→0 =

(11)

N k_BT

Из уравнения (7) видно, что вероятность появле-

ния энергии в ядре зависит только от парных взаи-

простой квадратной решетки размером N = 10⁴ спи-

модействий, в которых участвуют спины ядра E^∗ и

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

не зависит от остальных E⁰. При известном распре-

ми. Размер ядра был 3 × 3. Для сравнения резуль-

делении p(E_i) мы можем определить конфигурацию,

татов данные были получены с помощью аналити-

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

ческого решения Фердинанда [21], алгоритма парал-

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

лельного отжига и алгоритма Ванга-Ландау. Кри-

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

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

ческие суммы для каждой границы и используем их

щины линии.

в ходе МК-семплирования. Так как E^∗ зависит и от

Мы применили наш метод для решения ак-

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

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

p(E_i) для каждой из 2^b граничных конфигураций. В

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

Письма в ЖЭТФ том 110 вып. 9 - 10

2019

К численному расчету фрустраций в модели Изинга

703

Рис. 2. (Цветной онлайн) Теплоемкость простой квад-

Рис. 3. (Цветной онлайн) Блок для прохода по систе-

ратной решетки 10⁴ спинов Изинга: гибридный - CHM ,

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

энтропийный - CW L, реплично-обменный - CP T и

взаимодействием. Красным выделены спины входящие

аналитический - CF методы. Магнитная восприимчи-

в ядро, синим - граничные спины ядра. Пунктирными

вость: гибридный χHM и реплично-обменный χP T ме-

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

тоды

действия.

(ИСЛ) в модели диполей. Рассчитывалась гекса-

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

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

поли расположены на ребрах гексагона, магнитный

вектор i-диполя задается как m_i

= S_im^∗i, таким

образом момент имеет только две ориентации ±m_i.

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

между магнитными моментами рассчитывалась по

следующей формуле

]

∑

[ (m_im_j)

(m_ir_ij )(m_j r_ij )

E=Da³

-3

(12)

r^3/2

r^5/2

i<j

где D = µ²/a³ есть константа дипольного взаимодей-

ствия, a - постоянная решетки и r_ij - радиус-вектор,

соединяющий спины i и j. Радиус диполь-дипольного

Рис. 4. (Цветной онлайн) Одно из возможных основных

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

состояний гексагональной решетки

вого льда ограничен координационной сферой R₃ (14

ближайших соседей) ввиду несущественного влия-

Поскольку основное состояние известно, можно

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

ввести количественную меру фрустрации в модели

ческое поведение системы [22]. Ядро представлено на

Изинга, как относительное число возбуждений в ос-

рис. 3 и включает в себя 6 диполей, в то время как

новном состоянии

граничных спинов 24.

E_max + 〈E〉(T)

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

P_f (T) =

(13)

2E_max

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

фигураций основного состояния для гексагональной

где P_f ∈ [0, 1], 〈E〉(T → 0) = E_gs, модуль значения

∑

решетки магнитных моментов. На рисунке 4 схема-

абсолютного минимума E_max = -E_min =

|E_ij |,

〈i,j〉

тично изображена конфигурация основного состоя-

суммирование 〈i, j〉 проводится по всем учитывае-

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

мым парам диполей. Для 1D, 2D и 3D ферромагне-

ки. Всего учитывалось N = 11250 диполей с перио-

тиков Изинга P_f (T = 0) = 0, фрустрации отсутству-

дическими граничными условиями.

ют. Расчеты показывают, что для системы диполей

Письма в ЖЭТФ том 110 вып. 9 - 10

2019

704

А. Г. Макаров, К. В. Макарова, Ю. А. Шевченко, П. Д. Андрющенко и др.

на гексагональной решетке с R₃ радиусом взаимодей-

т.е. изменение параметра фрустраций P_f (T ) и вы-

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

числить долю возбуждений при известном значении

солютный минимум энергии E_min = -6.802 N ·D, то-

абсолютного минимума.

гда как энергия основного состояния E_gs = -2.441 N ·

В заключение необходимо отметить, что разра-

D. Таким образом, в соответствии с (13), основное

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

состояние, реализуемое за счет конфигураций рис.4,

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

фрустрировано на 32 %.

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

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

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

тельного числа возбуждений

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

стемы спинов Изинга к основному состоянию при

ΔP_f (T ) = P_f (T₂) - P_f (T₁)

(14)

T → 0, и его возможности необходимо исследовать

более глубоко. Также метод легко обобщается на спи-

и энтропии

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

∫ T2

C(T )

вольным Гамильтонианом, произвольным разбавле-

ΔS(T ) = S(T₂) - S(T₁) =

(15)

нием и внешним полем.

Работы выполнены при финансовой поддерж-

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

ке Министерства науки и высшего образования,

на рис. 5.

государственные задания

#075-00400-19-01 и

3.7383.2017/8.9. Реализация мультиспинового ги-

бридного метода выполнена А. Макаровым при

финансовой поддержке Российского фонда фун-

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

проекта

#18-32-00713. Часть исследования, вы-

полняемого П. Андрющенко, была поддержана

проектом Российского фонда фундаментальных ис-

следований # 18-32-00557. Для выполнения расчетов

были использованы вычислительные ресурсы ЦКП

“Центр данных ДВО РАН” [25].

1. A. P. Ramirez, Annual Review of Materials Science 24,

453 (1994).

Рис. 5. (Цветной онлайн) Температурное поведение

2. A. Mizuno, Y. Shuku, M. M. Matsushita, M. Tsuchiizu,

приращения энтропии и параметра фрустраций гекса-

Y. Hara, N. Wada, Y. Shimizu, and K. Awaga, Phys.

гональной решетки 11250 диполей

Rev. Lett. 119, 057201 (2017).

3. J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, Rev.

Как видно из рис.5, существует интервал темпе-

Mod. Phys. 82, 53 (2010).

ратур, где скорость роста числа фрустраций увели-

4. J. Reuther, R. Thomale, and S. Trebst, Phys. Rev. B

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

84, 100406 (2011).

медляется. Введенная в нашей работе количествен-

5. I. S. Hagemann, Q. Huang, X. P. A. Gao, A. P. Ramirez,

ная мера фрустраций может использоваться для экс-

and R. J. Cava, Phys. Rev. Lett. 86, 894 (2001).

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

6. L. C. Chapon, P. G. Radaelli, H. Zheng, and

нового льда экзотических решеток наноостровковых

J. F. Mitchell, Phys. Rev. B 74, 172401 (2006).

магнетиков, таких как, например, см. здесь [23, 24].

7. G. H. Wannier, Phys. Rev. 79, 357 (1950).

При выполнении необходимых условий к экспери-

8. G. H. Wannier, Phys. Rev. B 7, 5017 (1973).

ментальным образцам, исследователь может вычис-

9. K. Kanô and S. Naya, Prog. Theor. Phys. 10, 158

лить E_min в модели взаимодействующих точечных

(1953).

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

10. L. Balents, Nature 464, 199 (2010).

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

11. M. J. Harris, S. T. Bramwell, D. F. McMorrow,

менение внутренней энергии

T. Zeiske, and K. W. Godfrey, Phys. Rev. Lett. 79,

∫ T2

2554 (1997).

E(T₂) - E(T₁) =

C(T )dT,

(16)

12. G. Toulouse, Commun. Phys. 2, 115 (1977).

Письма в ЖЭТФ том 110 вып. 9 - 10

2019