Физика Земли, 2023, № 2, стр. 193-211

ВВЕДЕНИЕ

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

Строгий подход к расчету ДС основывается на методах микромагнетизма, основы которого были заложены в работе [Landau, Lifshitz, 1935], однако из-за сильной нелинейности и громоздкости микромагнитных уравнений только относительно простые задачи могли быть решены аналитически. Тем не менее, в рамках этого подхода были получены важные результаты по расчету предела однодоменности (ОД) и обнаружены конфигурации типа flower и vortex, возникающие при увеличении размера частиц несколько выше ОД предела, который для магнетита при комнатной температуре Т_r оценивается в 50–100 нм [Dunlop, Özdemir, 1997], в зависимости от формы частицы и ее структурных особенностей. В дальнейшем, благодаря развитию вычислительной техники, эти результаты были подтверждены численным расчетом микромагнитных уравнений и моделированием трехмерных ДС [Williams, Dunlop, 1989; Shcherbakov et al., 1990b; Fabian et al., 1996; Fukuma, Dunlop, 2006; Щербаков, Сычева, 2001; Dunlop, 1977; Enkin, Dunlop, 1987; Moon, Merrill, 1984; 1985; Shcherbakov et al., 1990a; Белоконь, Харитонский, 1985].

Отметим, что существование структур типа flower обусловлено исключительно неоднородностью размагничивающего поля и по этой причине они возникают с неизбежностью в любых частицах (кроме эллипсоидальных) при условии, что этот размер достаточно мал. При этом их остаточный магнитный момент лишь незначительно отличается от магнитного момента ОД частиц, поэтому на практике такие частицы просто отождествляются с ОД зернами, а за предел однодоменности принимается переход flower → vortex. Классические ОД структуры возникают только в эллипсоидальных частицах и в них за предел однодоменности принимается переход ОД → vortex.

В последние годы возможности микромагнитного моделирования неизмеримо возросли за счет развертывания графических процессоров (graphics processing unit, GPU), которые ускорили вычисления в сотни раз по сравнению с CPU-компьютерами (central processing unit, CPU, центральные процессоры общего назначения), что позволило проводить исследования, которые ранее были невозможны. Достижения последних лет в области компьютерной архитектуры, программирования, численных методов привели к развитию программных кодов для решения задач микромагнетизма – это, в частности, программные пакеты OOMMF [Fu et al., 2016], MuMax3 [Leliaert et al., 2018], Fidimag [Bisotti et al., 2018], MERRILL [Ó Conbhuí et al., 2018] и др. Соответственно, стали появляться теоретические работы [Lascu et al., 2018; Nagy et al., 2017; 2019; и др.] по изучению вихревых структур в частицах магнетита размером до 3 мкм микромагнитным методом. В частности, результатом этих работ оказалось обнаружение двухвихревой структуры (double vortex).

В данной работе мы детально исследуем эволюцию ДС частиц магнетита субмикронных и микронных размеров с температурой при охлаждении образцов от температуры Кюри Т_с до комнатной температуры Т_r и обратном их нагреве до Т_с. Такие расчеты прежде всего важны для понимания особенностей образования полной и парциальной термоостаточной намагниченностей (TRM и pTRM, соответственно), используемых при определении напряженности древнего геомагнитного поля B_anc по методике Телье. С этой целью нами выполнены расчеты методами микромагнетизма ДС зерен магнетита субмикронных и микронных размеров путем численного моделирования ДС зерен при их “охлаждении и нагреве” от температуры Кюри Т_с до T_r и обратно с целью определения степени необратимости изменений ДС при термообработке образцов и их возможного влияния на свойства термоостаточной намагниченности.

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

Пусть (x, y, z) – система координат, а α_x, α_y, α_z – направляющие косинусы вектора спонтанной намагниченности с осями x, y, z. В качестве модели рассмотрим для простоты кубическое зерно магнетита с ребрами вдоль осей (x, y, z).

Полная энергия частицы

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

Для проведения численных расчетов перейдем от континуальной модели к дискретной. Проведем мысленное разбиение кубической частицы с ребром d на N³ меньших кубиков, где N – число разбиений ребра куба, размер ребра одного кубика a = d/N [Щербаков, Сычева, 2001]. Пусть целые числа (i, j, k) являются координатами каждой ячейки в трех ортогональных осях x, y, z. Поскольку в пределах каждой ячейки направление вектора намагниченности предполагается постоянным, и div M_s = 0 (т.е. отсутствуют объемные магнитные заряды), расчет магнитостатической энергии в этом случае определяется поверхностными магнитными зарядами, сосредоточенными на гранях каждой из ячеек. Расчет направляющих косинусов вектора намагниченности M_s производился в сферической системе координат θ, φ:

Детали расчета полной энергии E_n, формулы для всех видов энергий выражения (1), а также выражения для энергий анизотропии, обменной энергии и магнитостатической энергии (i, j, k)-той ячейки приведены в работе [Щербаков, Сычева, 2001], отметим лишь, что в дальнейших расчетах мы пренебрежем вкладом второго члена в формуле для энергии кубической анизотропии в силу малости K₂ в сравнении с K₁.

Температурная зависимость магнитных констант для магнетита была взята в следующем виде: $A\left( T \right) = A\left( {{{T}_{{\text{r}}}}} \right)j_{{\text{s}}}^{{1.7}}$, $K\left( T \right) = K\left( {{{T}_{{\text{r}}}}} \right)j_{{\text{s}}}^{{8.5}}~,$ [Dunlop, Özdemir, 1997]. Здесь: T – температура; T_r – комнатная температура; ${{j}_{{\text{s}}}} = {{{{M}_{{\text{s}}}}\left( T \right)} \mathord{\left/ {\vphantom {{{{M}_{{\text{s}}}}\left( T \right)} {{{M}_{{\text{s}}}}\left( {{{T}_{{\text{r}}}}} \right)}}} \right. \kern-0em} {{{M}_{{\text{s}}}}\left( {{{T}_{{\text{r}}}}} \right)}}$ – приведенная спонтанная намагниченность. Приближенно зависимость величины спонтанной намагниченности магнетита от температуры можно описать формулой $~{{j}_{{\text{s}}}}\left( T \right) = {{\left[ {{{\left( {{{T}_{{\text{c}}}} - T} \right)} \mathord{\left/ {\vphantom {{\left( {{{T}_{{\text{c}}}} - T} \right)} {\left( {{{T}_{{\text{c}}}} - {{T}_{{\text{r}}}}} \right)}}} \right. \kern-0em} {\left( {{{T}_{{\text{c}}}} - {{T}_{{\text{r}}}}} \right)}}} \right]}^{{0.42}}}$ [Moon, Merrill, 1988]. При расчетах были использованы следующие численные значения магнитных параметров: T_c = 585°C, T_r = 20°C, обменная константа A(T_r) = 1.33 × 10^–11 Дж/м [Heider, Williams, 1988], константа кубической анизотропии K₁(T_r) = –1.36 × 10⁴ Дж/м³, намагниченность насыщения M_s(T_r) = 4.8 × 10⁵ А/м [Dunlop, Özdemir, 1997].

Задача определения стабильной конфигурации вектора намагниченности M_s состоит в минимизации полной энергии (1) как функции углов θ(i, j, k) и φ(i, j, k) для всех i, j, k = 1…N (т.е. полное число переменных для всего зерна есть 2N³). На деле для удобства минимизировалась приведенная энергия ${{E}_{n}} = {E \mathord{\left/ {\vphantom {E {\left( {M_{{\text{s}}}^{2}{{d}^{3}}j_{{\text{s}}}^{2}} \right)}}} \right. \kern-0em} {\left( {M_{{\text{s}}}^{2}{{d}^{3}}j_{{\text{s}}}^{2}} \right)}}$. Минимизация проводилась методом сопряженных градиентов. Для проверки стабильности каждого состояния в начальное состояние (магнитные моменты ячеек направлены вдоль легкой оси [111]) случайным образом вносилось небольшое возмущение, при каждом значении температуры после завершения процедуры минимизации полной энергии в величины полярных углов θ(i, j, k) и φ(i, j, k) вновь вносилось случайным образом возмущение, и поиск минимума энергии E_n при данной температуре повторялся, затем при полученной ДС осуществлялся переход к расчету при следующем значении температуры.

Для микромагнитного моделирования важно тонкое пространственное разбиение частицы на ячейки, максимально допустимый размер l_ex которых определяется равенством характерных плотностей магнитостатической ${{{{{{\mu }}}_{0}}M_{{\text{s}}}^{2}} \mathord{\left/ {\vphantom {{{{{{\mu }}}_{0}}M_{{\text{s}}}^{2}} {2~}}} \right. \kern-0em} {2~}}$ и обменной энергий $\frac{A}{{l_{{{\text{ex}}}}^{2}}}$ и, следовательно, описывается в терминах “длины обмена” ${{l}_{{{\text{ex}}}}} = \sqrt {{{2A} \mathord{\left/ {\vphantom {{2A} {\left( {{{{{\mu }}}_{0}}M_{{\text{s}}}^{2}} \right)}}} \right. \kern-0em} {\left( {{{{{\mu }}}_{0}}M_{{\text{s}}}^{2}} \right)}}} $ [Rave et al., 1998]. Здесь μ₀ – магнитная постоянная. Для магнетита при комнатной температуре величина “длины обмена” l_ex составляет 9–10 нм. С учетом цитированной выше температурной зависимости A(T) ∝ $j_{{\text{s}}}^{{1.7}}$, с ростом температуры l_ex медленно растет пропорционально $j_{{\text{s}}}^{{ - 0.15}}$.

РЕАЛИЗАЦИЯ МИКРОМАГНИТНЫХ РАСЧЕТОВ

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

Последовательная часть программ, выполняемых на CPU, разработана на языке высокого уровня с помощью лицензионного программного обеспечения Intel Parallel Studio XE Professional Edition for Fortran Linux. Intel Parallel Studio XE повышает производительность приложений, работающих на платформах Intel, и является лучшим в своем классе компилятором. Для параллельных вычислений на GPU был написан интерфейс к используемым функциям CUDA API с помощью модуля ISO_C_BINDING – части стандарта Fortran 2003 и более поздних версий. Модули программ для параллельных вычислений написаны на языке CUDA С. Технология CUDA (Compute Unified Device Architecture) компании Nvidia – бесплатный эффективный инструмент для программирования на GPU, которая облегчает написание GPU-приложений и предназначена для разработки приложений для параллельных вычислительных устройств. Параллельные расчеты выполнены на GPU Tesla T4, ориентированном на высокопроизводительные вычисления, версия CUDA 11.4, CPU представлен двумя процессорами Intel Xeon Gold 6240.

Быстрое преобразование Фурье (FFT) для расчета магнитостатической энергии выполнено с помощью дополнительной библиотеки, CUDA-адаптированном варианте набора инструментов для осуществления FFT - CUFFT, входящей в набор разработчика CUDA Toolkit. Для оптимизации времени расчетов на GPU с целью минимизации обращений к глобальной памяти GPU в программе также использовалась разделяемая и константная память GPU. Распараллеливание вычислений позволило увеличить скорость счета для массивов 64³ приблизительно в 5 раз, а для массивов 128³ в 30 раз по сравнению с последовательным счетом. Поскольку выигрыш по времени счета на GPU увеличивается с ростом размера массивов, то для массивов 192³ и 256³ получаем выигрыш в скорости счета, предположительно, не менее, чем в 100 раз, что позволило нам провести расчеты для больших трехмерных массивов и получить ДС частиц магнетита размером до 2 мкм. Загрузка GPU при работе программы составляла 60-100% на разных этапах счета. Здесь надо отметить, что, если при расчетах число ячеек вдоль ребра куба N = 255, следовательно, число граней для расчета Е_ms есть (N + 1)³ = 256³, а размер каждого из 6 комплексных массивов с помощью FFT составляет (2(N + 1))³ = 512³.

РЕЗУЛЬТАТЫ

В этом разделе изложены результаты численного моделирования магнитной конфигурации зерен при их “охлаждении и нагреве” от температуры Кюри Т_с до комнатной температуры T_r и обратно. Как известно, для частиц сферической формы существует строгая оценка критического размера однодоменности, а именно переход из ОД состояния в моду curling (vortex) [Eisenstein, Aharoni, 1976]:

Подставляя сюда численные значения параметров, получаем, что при комнатной температуре d₀ = 63.1 нм. При температурах, близких к Т_с, можно пренебречь членом 4K/3 в знаменателе (3) в силу его быстрого убывания с ростом температуры, тогда d₀(T) ∝ l_ex ∝ $j_{{\text{s}}}^{{ - 0.15}}$. Иными словами, в непосредственной окрестности T_c все зерна теоретически в какой-то момент становятся однодоменными и по этой причине во всех расчетах первоначальное доменное состояние полагалось однодоменным. На деле, конечно, ситуация является намного более сложной в силу того, что в игру здесь вступает много разных факторов – размагничивающие поля, термофлуктуации, напряжения, неоднородности состава, примеси и т.д., но здесь мы отвлечемся от этих сложностей, поскольку они мало влияют на дальнейшее изложение.

Расчеты начинались от T = 584.999°С (j_s = 0.0038), уменьшение (увеличение) температуры проводилось небольшими скачками с интервалом Δj_s = = 0.008 (T = 584°С), Δj_s = 0.004 в области температур (584–553.4)°С, Δj_s= 0.02 от 552.4°С до T_r. Из изложенного во ВВЕДЕНИИ ясно, что в этом интервале размеров следует ожидать появления конфигураций типа flower и вихревых структур типа vortex и double vortex. Конкретно в каждом случае возможные конфигурации определяются температурной и магнитной предысторией частицы [Щербаков, Сычева, 2001], так что в одном и том же зерне при одной и той же температуре могут реализоваться с определенной степенью вероятности либо flower, либо vortex как метастабильные состояния. С целью учета такой возможности нами рассчитывались следующие три варианта.

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

2. “Провокация” перехода в моду vortex (V), когда при каждом шаге вносились малые возмущения в виде моды V. Возмущение для всех подобных расчетов выбиралось в виде а × vortex(T_r), где vortex(T_r) есть вихревая конфигурация, полученная в ходе расчетов по варианту 1) при T = T_r, а а есть малый коэффициент, который для определенности в наших расчетах полагался равным 0.1. В этом случае переход в состояние vortex означает преодоление потенциального барьера, отделяющего конфигурацию flower от состояния vortex с помощью внешнего воздействия. Ниже мы кратко обсудим вопрос о том, как величина а влияет на температуру перехода из одного состояния в другое.

3. “Провокация” перехода в мультивихревую моду W (в нашем случае это двойной vortex), величина возмущений та же, что в варианте 2.

Надо отметить, что “провокации” V и W вносились при “охлаждении” частицы от T_c до T_r, а при “нагреве” до T = 584°С – для всех 3 вариантов расчета вносилась провокация перехода в моду flower: а × flower(T_r), a = 0.05, в противном случае, как показали расчеты, конфигурация, полученная для T = T_r, при “нагреве” не изменялась вплоть до T = 584°С.

Визуализация трехмерных структур затруднена в силу большого количества ячеек, на которые разбита частица, поэтому для лучшего понимания полученных структур мы, кроме вывода картины распределения магнитных моментов, применили подход, предложенный в работе [Nagy et al., 2019], изобразив картины срезов анизотропного отклонения моментов ячеек (anisotropic deviation of moment, ADM) в различных плоскостях:

где: –1/3 ≤ Θ ≤ 0, значение Θ = 0 соответствует направлению магнитного момента вдоль оси трудного намагничивания [100], а Θ = –1/3 соответствует направлению магнитных моментов вдоль легкой оси [111] (для магнетита).

Для частиц размером менее 80 нм при расчетах по варианту 1 перехода flower → vortex не происходит, по варианту 2 – переход flower → vortex получен для частиц размером 70 нм и выше. Для зерна размером d = 80 нм изменение конфигурации магнитного момента в зависимости от температуры представлено на рис. 1 и рис. 2. Как можно видеть из рис. 2, при начальном ОД состоянии (T = 584.999°С) магнитный момент частицы направлен строго вдоль легкой оси [111], и такая конфигурация сохраняется до T ≈ 582°С. Далее развивается структура flower, которая при дальнейшем “охлаждении” переходит в состояние с ярко выраженным одним вихрем, то есть возникает конфигурация типа vortex (V). Для варианта 1 (без внесения возмущений) это происходит при Т ≈ 417°С (рис. 1, рис. 2). На самом деле, состояние vortex сосуществует с состоянием flower и при более высоких температурах. Это можно увидеть, если “спровоцировать” переход из метастабильного состояния flower в более стабильную моду vortex путем внесения соответствующего возмущения, как это описывалось выше (вариант 2). В этом случае трансформация ОД→V происходит намного раньше, при Т ≈ 582°С.

Заметим, что при переходе в состояние vortex для частицы размером 80 нм сохраняется достаточно большой остаточный момент M_r ≈ 0.6M_s (рис. 1), постепенно уменьшающийся с температурой, так что при T = T_r отношение M_r/M_s достигает величины ≈ 0.5. Для d = 90 нм при T = T_r отношение M_r/M_s ≈ 0.3, что типично для малых ПОД частиц. Существенно также, что структура vortex для размеров частиц d ≥ 90 нм остается стабильной как при “охлаждении” до T_r, так и при “нагреве” вплоть до T = 584°С и повторном “охлаждении” до T_r.

Распределение магнитного момента в кубической частице идеальной формы в зависимости от температуры, аналогичное описанному выше для частицы размером 80 нм, наблюдается для всех вариантов 1)–2) в частицах размером до 200 нм включительно. При расчетах по провокации W, которые проводились для частиц размером выше 120 нм, в частице размером d = 150 нм образуется конфигурация, состоящая преимущественно из одного крупного вихря вдоль направления [111], но с дополнительными небольшими по объему вихревыми образованиями. Эта конфигурация также остается стабильной при последующих “нагревах–охлаждениях”. На диаграмме рис. 9 эта мода отмечена аббревиатурой SV.

С увеличением размера до 170 нм и выше при расчетах с провокацией W наблюдается уже образование двухвихревой конфигурации W, на рис. 3, рис. 4 представлены результаты расчетов для d = = 200 нм. При расчете без возмущений (вариант 1) переход flower – vortex протекает через серию метастабильных состояний, до T = 559°С развивается flower, при T = 558°С возникает вышеупомянутая структура SV, а эта конфигурация в процессе дальнейшего “охлаждения” к T ≈ 388°С преобразуется в одновихревую структуру vortex, стабильную к последующим “нагревам”–”охлаждениям”. При провокации V при “охлаждении” до T ≈ 583°С мода flower переходит в конфигурацию, близкую к vortex, затем при дальнейшем “охлаждении” плавно трансформируется в vortex. Сравнение энергий E_n(T_r) различных состояний при комнатной температуре показывает, что для вариантов 1 и 2 эти энергии близки, в то время как энергия конфигурации W оказывается заметно выше (рис. 3).

На рис. 5 и рис. 6 представлены результаты расчетов для частицы размером 1 мкм. Конфигурация, образованная при расчетах с провокацией W (рис. 6), состоит, в основном, из двух крупных вихрей вдоль направления [111], закрученных в противоположных относительно друг друга направлениях, ближе к граням куба образуются более мелкие вихревые образования и замыкающие домены. В этом случае M_r/M_s = 0.05. При расчетах без возмущений (вариант 1) получена также конфигурация W, M_r/M_s = 0.07. При расчетах с провокацией V получена одновихревая конфигурация vortex, M_r/M_s = 0.03, что ниже многодоменного (МД) предела 0.05. Отметим, что низкие значения M_r/M_s, возможно, являются артефактом, вызванным неучетом реальной (несимметричной) формы частиц, неоднородностями состава, напряжениями и т.д. Сравнение энергий E_n(T_r) различных состояний при комнатной температуре показывает, что для вариантов 1 и 2 эти энергии близки, как и по результатам расчетов для частицы размером 200 нм, но энергия конфигурации W в этом случае в два раза ниже (рис. 5).

Таким образом, по результатам расчетов получено, что для варианта 1 (без возмущений) при “охлаждении” до комнатной температуры в частицах размерами d = 80–200 нм образуется мода V, d = 300 нм – мода SV, d = 500–2000 нм – мода W. Для варианта 3 (провокация W): d = 120 нм – мода V, d = 150 нм – мода SV, d = 170 – 2000 нм – мода W. Отметим, что конфигурации, объединенные под символом W, несомненно, несколько отличаются для разных размеров частиц и условий расчета, что можно видеть, сравнивая конфигурации на рис. 4 (d = 200 нм) и рис. 6 (d = 1 мкм). Для варианта 2 (провокация V) для частиц размерами d = 70–1500 нм образуется конфигурация vortex, для d = 2 мкм получена конфигурация SV. О возможности существования одновихревых структур в частицах магнетита микронных размеров отмечено в работе [Nagy et al., 2019]. На рис. 7 представлены результаты моделирования с провокацией V для частицы размером d = 1.5 мкм.

Как отмечено выше, возмущение V выбиралось в виде а × vortex(T_r), где a = 0.1. Мы провели расчеты для a < 0.1 с целью выяснить, как величина а влияет на температуру перехода из одного состояния в другое на примере перехода flower → vortex (рис. 8) для частицы размером 80 нм. Как и следовало ожидать, при уменьшении a происходит постепенный сдвиг температуры трансформации T_tr(a) вниз, пока при а → 0 эта температура не совпадет с таковой, полученной без всякого возмущения. Иными словами, величина параметра а в какой-то мере отражает величину потенциального барьера, разделяющего конфигурации flower и vortex при T = T_tr(a).

Как общий результат проведенного моделирования доменной структуры кубических частиц магнетита субмикронных и микронных размеров на рис. 9 приведена диаграмма зависимости температуры переходов “flower – vortex”, “flower - multivortex” от размера частиц и начального возмущения. В области размеров и температур, расположенной левее и ниже серого пунктира мода “flower” является единственно возможной (по результатам данной работы). В области между пунктирными линиями возможно существование как моды flower, так и вихревой структуры – V (до d = 150 нм), для размеров d = 150–300 нм – flower, V, SV, W, при d > 300 нм – моды “flower”, V, W. Правее и выше черной пунктирной линии, до размера d = 150 нм – находится область, где реализуется исключительно состояние единичного вихря V, в интервале размеров d = 150–300 нм – возможно образование как моды V, так и W, SV, а выше d = 300 нм – образуются моды W и V, для d = 2 мкм получена также конфигурация SV (провокация V).

ДИСКУССИЯ

Автор работы [Néel, 1955] предположил, что для ПОД и МД зерен физические механизмы образования термоостаточной намагниченности могут быть разделены на гистерезисные и термофлуктуационные. Гистерезисный механизм пренебрегает тепловыми возмущениями и предполагает, что TRM приобретается за счет переброса вектора магнитного момента по направлению внешнего магнитного поля, либо за счет движения доменных границ (ДГ) во внешнем магнитном поле, приложенном к образцу во время охлаждения породы. Движущей силой движения ДГ в ферримагнитной частице является зависящее от температуры изменение баланса между энергией внешнего поля, энергией размагничивания и энергией доменных стенок. Термофлуктуационный механизм также предполагает переброс вектора магнитного момента или движение ДГ во внешнем магнитном поле как причину образования TRM, но за счет действия термофлуктаций. При этом вводится понятие температуры блокировки T_b, выше которой тепловые возбуждения поддерживают термодинамическое равновесие, а ниже положение ДГ или направление магнитного момента m замораживается [Néel, 1955; Schmidt, 1973]. Однако, на самом деле, ни гистерезисная, ни термофлуктуационная теории не вполне согласуются с наблюдениями [Shcherbakov et al., 1993; Shcherbakov, Shcherbakova, 2001; McClelland, Shcherbakov, 1995; Shcherbakova et al., 2000]. Термофлуктуационная модель, построенная по прямой аналогии с теорией образования TRM в ОД зернах, предсказывает наблюдаемую линейную зависимость TRM для слабых магнитных полей, но не может объяснить экспериментально наблюдаемые в ПОД и МД частицах хвосты парциальных TRM (pTRM) и зависимость интенсивности pTRM от термической предыстории, т.е. разницу между величинами pTRMa, полученной при охлаждении образца от T_с, и pTRMb, полученной при охлаждении образца от T₁ < T_с [Shcherbakov, Shcherbakova, 2001], что очевидным образом противоречит закону Телье аддитивности pTRM, сформулированному в работе [Thellier, Thellier, 1959]. С другой стороны, гистерезисная модель объясняет наличие хвостов, но предсказывает нелинейную зависимость TRM от величины поля и не в состоянии объяснить установленный экспериментально закон аддитивности pTRM, справедливый в определенных условиях (одинаковая температурная история создания всех pTRM) и для ПОД, и для МД зерен.

Поскольку линейность TRM(B) в достаточно слабых магнитных полях является экспериментально установленным фактом для образцов, содержащих не только ОД, но и ПОД и МД частицы, то гистерезисная модель в настоящее время в литературе серьезно не обсуждается, а основные усилия направлены на развитие термофлуктуационной модели применительно к ПОД и МД частицам. Некоторый прогресс в этом направлении был связан с развитием кинематической теории TRM [McClelland, Sugiura, 1987; Shcherbakov et al., 1993] и феноменологической моделью [Fabian, 2001], позволяющей визуализировать некоторые свойства pTRM. Обе эти модели берут за основу положение о том, что блокирующая T_b и деблокирующая T_db температуры в ПОД и МД частицах не равны друг другу. Напомним, что равенство этих температур лежит в основе теоретического обоснования справедливости законов аддитивности и независимости Телье, поэтому нарушение этого равенства тесно связано с упомянутыми выше отклонениями от этих законов для pTRM ПОД и МД зерен.

Физическая природа этих отклонений, очевидно, связана со сложностью ДС ПОД и МД частиц. Как следует из рис. 9, типичная ДС для ПОД частиц, вплоть до d = 2 мкм, представлена модой flower (в непосредственной окрестности T_c) и различными вихревыми конфигурациями на остальном температурном интервале. Тем самым наши расчеты подтверждают результаты работ [Almeida et al., 2014; 2016; Nagy et al., 2017; 2019; Lascu et al., 2018], в которых исследовались вихревые структуры в частицах магнетита размером от 100 нм до 3 мкм.

В работах [Almeida et al., 2014; 2016] экспериментально исследовались частицы магнетита размером 150–250 нм и были построены карты магнитной индукции, полученные по вихревой магнитной структуре в частицах ПОД магнетита размером 200 нм. Эксперимент показал, что вихревая структура устойчива при нагревах до температур, близких к T_c. Этот результат согласуется с нашими расчетами, представленными на рис. 6 и рис. 7. Авторы работы [Khakhalova et al., 2018] также сообщают об экспериментальных наблюдениях вихревых конфигураций, но уже в зернах титаномагнетита значительно бóльшего размера, 1–5 мкм.

Таким образом, наиболее существенным отличием малых ПОД частиц (размером до 2–3 мкм) от однодоменных является наличие в них ряда возможных состояний типа flower и vortex c одним или двумя (возможно, несколькими) вихрями. При этом надо иметь в виду, что количество таких состояний может быть достаточно велико, поскольку различные состояния могут реализовываться вдоль различных осей кубической решетки. Так, в работе [Fabian, Shcherbakov, 2018] авторы показали, что в кубической частице магнетита размером около 50 нм в зависимости от ее предистории могут реализовыватся моды flower и vortex, при этом число метастабильных состояний n_s оказалось равным 60. Используя данные этих расчетов, авторы работы [Fabian, Shcherbakov, 2018], построили модель образования вязкой намагниченности в таких частицах на основе термофлуктуационной модели Нееля с использованием кинетических уравнений вида:

где: t – время; y_i - относительная заселенность уровней с условием нормировки $\sum\nolimits_{i = 1}^{{{n}_{s}}} {{{y}_{i}}} = 1$_;f_ik – частоты переходов из состояния i в состояние k. Аналогичные уравнения для случая n_s = 4 использовались в работе [Афремов, Панов, 2004] для исследования эволюции распределения магнитных состояний двухфазных частиц под действием термофлуктуаций. Подчеркнем, что в обеих работах рассматривались процессы, проходящие при некоторой постоянной температуре, когда коэффициенты f_ik являются постоянными, что позволяет получить полное решение системы (5) в аналитическом виде, например, в форме матричной экспоненты.

Термофлуктуационная модель Нееля TRM ОД частиц также основана на использовании кинетических уравнений типа (5) в простейшем случае n_s = 2, но с коэффициентами f_ik, зависящими от температуры, которая в этом случае изменяется со временем по мере охлаждения или нагрева. Строгий расчет свойств TRM для этой модели был недавно представлен в работе [Shcherbakov et al., 2021] на основе точного решения уравнений (5) для случая n_s = 2 и линейного (или экспоненциального) спада температуры со временем. Очевидно, этот же подход можно пытаться распространить и на ПОД частицы через введение в рассмотрение ряда дискретных магнитных конфигураций, между которыми происходят переходы под совместным влиянием внешнего магнитного поля и термофлуктуаций. Однако изменчивость во времени коэффициентов уравнений (5) сильно затрудняет решение такой системы при n_s > 2, что составляет ее существенное отличие от моделей с постоянными коэффициентами, рассматриваемых в работах [Афремов, Панов, 2004; Fabian, Shcherbakov, 2018], решение которых можно получить в виде матричной экспоненты либо как сумму экспоненциальных членов.

Для понимания некоторых особенностей приобретения TRM в малых ПОД частицах рассмотрим для простоты схему с всего двумя возможными магнитными состояниями, предполагая, например, сосуществование в одной частице только одно- и двух- вихревых конфигураций, направленных по одной и той же оси легкого намагничивания. Тогда можно построить простейшую теоретическую схему образования TRM в малых ПОД частицах, рассматривая систему уравнений (5) для случая четырех возможных магнитных состояний: 1 = V₊; 2 = V_–; 3 = W₊ ; 4 = W_–, где индексы + и – означают направления вдоль или против внешнего поля h. Очевидно, прямые переходы V₊ → W_– и V_– → W₊ и обратно требуют преодоления гораздо бóльшего потенциального барьера, так как эти состояния имеют противоположные направления намагниченности и требуют полной переконфигурации ДС при таком переходе. По этой причине в дальнейшем мы пренебрежем вероятностью этих переходов. Кроме того, по аналогии с работой [Shcherbakov et al., 2021] примем линейную зависимость T(t), так что при охлаждении T = T_c(1 – t/t₀), где t₀ ~ 10² c – характерное время охлаждения. Тогда соответствующие кинетические уравнения можно записать следующим образом:

Здесь и ниже обозначение t следует понимать как безразмерное время, в силу нормировки t → t/t₀. Частота переходов из состояния i в состояние k есть:

где: q = ln(f₀t₀); f₀ ≈ 10⁹ с – частота прецессии (попыток преодоления потенциального барьера); m_i – магнитный момент i-го состояния; B_ik – критическое поле перехода i → k; k_B – постоянная Больцмана. Здесь предполагается, как обычно, что m_ih/k_BT $ \ll $ 1. Это неравенство обеспечивает линейность индуцированной намагниченности по полю, что является и условием линейности TRM(B). Другим фактором, обеспечивающим линейность, может быть магнитостатическое взаимодействие [Щербаков, Щербакова, 1975], но обсуждение этого вопроса лежит далеко за рамками настоящей статьи.

При B $ \ll $ B_ik магнитные моменты частиц в состояниях 1 и 2, а также в состояниях 3 и 4 практически равны друг другу, так что можно ввести m_V = m₁ = m₂ и m_W = m₃ = m₄. Обозначим:

Очевидно, что B₁₂ = B₂₁ и B₃₄ = B₄₃. Тогда коэффициенты (7) можно записать в виде:

Здесь a_ij = exp(q – m_i(T)B_ij(T)/k_BT) (i, j = V,W) есть частоты переходов V → V, V → W, W → V. Теперь, используя условие нормировки $~\sum\nolimits_{i = 1}^4 {{{y}_{i}}} = 1$, можно свести систему (6) к системе трех линейных неоднородных дифференциальных уравнений:

Далее для простоты положим, что зависимость от температуры параметров m_i(T) и B_ij(T) пропорциональна j_s(T) (по аналогии с ОД зернами), то есть:

Здесь m_sd = M_sv – магнитный момент ОД зерна при T = T_r, равного по объему рассматриваемому ПОД зерну. Коэффициенты S_V и S_W (<1) учитывают уменьшение магнитного момента за счет неоднородного распределения намагниченности по его объему в ПОД частицах.¹¹ Тогда:

Пример численного расчета решения системы уравнений (10) с коэффициентами, определенными по формулам (12), приведен на рис. 10 для параметров b_vv = 150, b_vw =20, b_wv = 50, b_ww = 60, b = 0.01, S_V = 0.5, S_W = 0.25. Начальные условия были выбраны в виде y_i = 0.25, i = 1, 3, предполагая для простоты что в окрестности T_c все состояния термодинамически равноправны.

Как видно из рис. 10, приобретение намагниченности в частице с двумя возможными конфигурациями регулируется тремя одновременно идущими процессами: намагничиванием каждой из одно- и двух- вихревой компонент и переходами из одной конфигурации в другую. В самом начале охлаждения мы видим бурное развитие всех этих трех процессов, но затем начинается быстрая перекачка одновихревой конфигурации в двухвихревую, в результате чего происходит видимый спад индуцированной намагниченности y_V(T), ассоциируемой с одновихревой конфигурацией. Приближение к нулю y_V(T) совпадает с практически полной трансформацией моды V в энергетически более выгодную моду W, и дальнейший рост полной намагниченности происходит исключительно за счет роста намагниченности этой моды, находящейся в суперпарамагнитном состоянии. Отметим, что для данного набора параметров блокировка перехода V → W происходит при t ≈ 0.9, а намагниченность y_W(T) блокируется заметно ниже, при t ≈ 0.75.

Рассмотрим теперь процесс терморазмагничивания уже заблокированной намагниченности. Пусть зерно было охлаждено до температуры T₁ (на рис. 10 T₁ = 0.5T_c), чему соответствует нормированная температура τ₁ = 1 – T₁/T_c = 0.5. При нагреве зависимость температуры от времени есть T(t) = T₁ + T_ct, где 0 < t < 1 – τ₁. Соответственно, величины (13) записываются как:

Кривая терморазмагничивания TRM, рассчитанная по уравнениям (10), но при b = 0 и с коэффициентами a_ij, рассчитанными по формуле (14), показана на рис. 10 крупным пунктиром. Параметры b_vv, b_vw, b_wv, b_ww оставались теми же, что и при расчете кривых намагничивания. Начальные условия были выбраны в виде y_i(0) = $y_{i}^{0}$, i = 1,3, где $y_{i}^{0}$ есть значения y_i(t₁), полученные при решении системы (10) для процесса охлаждения (рис. 10). Сравнение кривых терморазмагничивания и намагничивания ясно показывает, что блокирующие и деблокирующие температуры TRM идентичны. Фактически рассмотренный случай мало отличается от картины образования TRM в ОД частицах, поскольку частицы с изначальной V-компонентой быстро трансформируются в более выгодную энергетически W-моду, и при дальнейшем охлаждении практически вся намагниченность оказывается связана именно с этой модой.

На рис. 11 показаны результаты подобного же расчета, но при параметрах b_vv, b_vw, b_wv и b_ww, выбранных таким образом, чтобы энергии этих состояний были примерно одинаковы и не происходило полной трансформации одной моды в другую. Кроме того, за счет большой разницы в величине параметров b_vv и b_ww блокировка намагниченности в подсистемах V и W происходит при заметно разных температурах. Как мы видим, в таком случае и деблокировка происходит в два этапа, но в обоих случаях наблюдается совпадение соответствующих T_b и T_db.

Оговоримся, что представленные здесь расчеты следует рассматривать только как попытку качественного описания механизма приобретения TRM ПОД зернами, никак не претендующее на количественное описание. Действительно, в наших расчетах как величина ключевых параметров a_ij(t) и h_i(t), так и их зависимость от температуры определялись достаточно произвольно, по аналогии с моделью TRM ОД частиц, [Shcherbakov et al., 2021]. Однако ясно, что для реального расчета процесса образования TRM необходимо знать не только возможные ДС и спектр их энергий, но и величину потенциальных барьеров, разделяющих различные ДС, как функцию T. Как показали авторы работы [Fabian, Shcherbakov, 2018] эта задача в принципе может быть решена с помощью метода эластичных лент для нахождения седловых точек (перевалов) на поверхности (многомерной) микромагнитного функционала энергии, но это – отдельная задача, лежащая далеко за рамками настоящей статьи.

Суммируя все сказанное, подчеркнем, что с физической точки зрения принятая здесь модель ПОД частиц предполагает наличие в них конечного набора n_s конфигураций, каждая из которых обладает определенной энергией, зависящей от температуры. Между этими конфигурациями возможны переходы, в чем и заключается фундаментальное отличие ПОД и ОД частиц. Осуществлению перехода мешает некий потенциальный барьер, который при достаточно высоких температурах легко преодолевается за счет термофлуктуаций, приводя ансамбль таких частиц в состояние термодинамического равновесия как между числом частиц в различных ДС, так и между состояниями вдоль и против направления внешнего поля в одной и той же конфигурации – расширенный аналог суперпарамагнитного состояния. По мере охлаждения ансамбля зерен переходы между конфигурациями (как и перевороты магнитных моментов внутри конфигураций) постепенно замораживаются. Иными словами, блокировка в такой модели осуществляется поэтапно, как это проиллюстрировано на рис. 10 и рис. 11. Существенно, что при нагреве частиц деблокировка заблокированных намагниченностей и конфигураций происходит в обратном порядке при той же температуре, что и их блокировка, что означает равенство блокирующих и деблокирующих температур. Действительно, в парадигме Нееля термофлуктуационного механизма приобретения TRM в ПОД зернах мы по-прежнему остаемся в рамках модели, описываемой конечной системой n_s кинетических уравнений (10), и хотя коэффициенты этой системы зависят от температуры, но эта зависимость не изменяется от того, идет ли охлаждение или нагрев ансамбля зерен, что и обеспечивает равенство T_b и T_db.

Вывод о справедливости законов Телье аддитивности и независимости pTRM в малых ПОД частицах субмикронного размера не противоречит экспериментальным данным [Большаков, Щербакова, 1979; Shcherbakova et al., 2000; Shcherbakov, Shcherbakova, 2001]. Однако остается вопрос, что же является причиной нарушения этих законов при увеличении размеров частиц? Можно предположить, что решающую роль здесь играет взрывной рост числа возможных магнитных конфигураций при росте размеров зерен – как уже говорилось выше, число возможных ДС в зернах, несколько превышающих ОД размер, уже может достигать нескольких десятков, а при дальнейшем увеличении размера их число может расти неограниченно. Можно предположить, что при очень большом числе различных ДС такая система потеряет термодинамическое равновесие и обратимость процессов перестройки ДС при охлаждении и нагреве, что и приведет к нарушению равенства T_b и T_db.

ВЫВОДЫ

1. Рассчитана магнитная конфигурация магнетитовых частиц кубической формы субмикронных и микронных размеров и ее эволюция с изменением температуры при охлаждении образцов от температуры Кюри Т_с до комнатной температуры и обратном их нагреве до Т_с.

2. Показано, что типичные магнитные конфигурации в частицах размера до 2 мкм представлены формами flower и vortex с одним или двумя вихрями.

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

4. Согласно этой модели, в малых ПОД частицах наблюдается равенство блокирующих и деблокирующих температур, и соответственно в них должны выполняться законы Телье аддитивности и независимости pTRM.