Studying the causes of similarity between blocking temperature spectra of chemical and thermoremanent magnetization by computer simulation

Full Text

ВВЕДЕНИЕ

Задачи палеомагнетизма решаются путем расшифровки информации о характеристиках геомагнитного поля, содержащейся в векторе естественной остаточной намагниченности образцов (NRM), такая расшифровка возможна лишь при условии сохранности первичной намагниченности образцов исследуемых горных пород. По этой причине при интерпретации палеомагнитных данных всегда возникает проблема их достоверности. При этом для целей надежного определения на вулканических породах палеонапряженности и палеонаправлений, кроме хорошей сохранности остаточной намагниченности и первичности выделяемого сигнала, требуется доказать еще и термоостаточную природу (TRM) естественной остаточной намагниченности пород (NRM). В то же время хорошо известно, что в природных условиях стабильная и значительная по величине остаточная намагниченность может возникать в горных породах и в результате физико-химических изменений магнитных минералов $–$ носителей остаточной намагниченности.

Химической остаточной намагниченностью (CRM) принято называть остаточную намагниченность, возникающую при длительной выдержке образца во внешнем поле B при некоторой постоянной температуре T_CRM вследствие минералогических изменений [Draeger et al., 2006]. Другое ее название $–$ кристаллизационная намагниченность [Doell, 1957]. Возможны два физических механизма возникновения CRM: a) в результате роста размера ферри (ферро)-магнитных суперпарамагнитных (СПМ) зерен, когда при выдержке при некоторой температуре Т_CRM и в присутствии внешнего поля В их объем вырастает до блокирующего объема v_b (модель Хайга [Хайг, 1962]); б) в результате роста температуры Кюри T_c вследствие окисления титаномагнетитовых зерен.

В работе [Smirnov, Tarduno, 2005] авторы предположили, что термоостаточная намагниченность TRM и CRM могут иметь схожую температурную стабильность, если обе намагниченности имеют узкий интервал блокирующих температур T_b. Если это так, то они будут обладать похожими спектрами T_b, и тем самым будут плохо различимы с точки зрения экспериментов по методикам Телье$–$Коэ или Вилсона$–$Буракова, давая в итоге ложные определения палеонапряженности. Однако с теоретической точки зрения эта гипотеза встречает возражения. Действительно, еще в книге [Stacey, Banerjee, 1974], в рамках модели [Хайг, 1962] была показана справедливость соотношения:

$\frac{\text{CRM}}{\text{TRM}}=\frac{K\left(T_b\right)}{K\left(T_{\text{CRM}}\right)}\frac{M_s\left(T_{\text{CRM}}\right)}{M_s\left(T_b\right)},$

где: M_s $–$ спонтанная намагниченность; T_b $–$ блокирующая температура частиц ансамбля при образовании TRM; K $–$ константа анизотропии; T_CRM $–$ температура кристаллизации. Принимая для простоты, что константа анизотропии определяется, главным образом, анизотропией формы K = NdMs²/2, где N_d $–$ параметр анизотропии, получим:

 

            CRM/TRM = M_s (T_b)/M_s(T_CRM) < 1,   (1)

поскольку, очевидно, что T_CRM < T_b и, стало быть, M_s(T_b) < M_s(T_CRM). (На самом деле, неравенство (1) справедливо и в общем случае [Stacey, Banerjee, 1974]). Из соотношения (1) следует, что отношение CRM/TRM падает с ростом T_b, то есть с ростом коэрцитивности и/или объема, а это означает, что теоретически спектры T_b для TRM и CRM, приобретенной по первому механизму, неидентичны [Shcherbakova et al., 2014].

Однако необходимо отметить, что соотношение (1) получено для невзаимодействующих однодоменных (ОД) частиц, а как показали теоретические расчеты [Афремов, Харитонский, 1988; Shcherbakov et al., 1996], интенсивность CRM существенно определяется магнитостатическим взаимодействием между зернами, средняя величина которого пропорциональна M_s и относительной объемной концентрации ферримагнетика c (≈ 5 Msс). В результате взаимодействие приводит к заметному снижению интенсивности как CRM, так и TRM. Эти выводы были подтверждены и результатами численного моделирования процессов образования CRM (по первому механизму) и TRM (на уже выросших до конечного объема зернах), выполненного уже с учетом магнитостатического взаимодействия [Shcherbakov et al., 1996; Щербаков и др., 2017].

В работе [Щербаков и др., 2017] было выполнено численное моделирование экспериментов определения палеонапряженности по методу Телье на ансамбле взаимодействующих ОД частиц, приобретших CRM за счет роста их размера от СПМ до ОД. Результаты моделирования показали, что практически на всех диаграммах Араи$–$Нагата, построенных по результатам такого моделирования, можно выделить протяженный, очень близкий к линейному, участок, по которому можно определить «палеонапряженность» Banc. Немаловажно отметить, что такой участок можно выделить и при моделировании экспериментов определения палеонапряженности по методу Телье на невзаимодействующих ОД зернах, что видимым образом находится в противоречии с соотношением (1). Подчеркнем, что линейность диаграммы Араи$–$Нагата в том или ином температурном интервале как раз и означает подобие в нем спектров Tb. Иными словами, данные компьютерного моделирования говорят в пользу справедливости гипотезы [Smirnov, Tarduno, 2005].

Такой же вывод о наличии квазилинейного участка на диаграммах Араи$–$Нагата был сделан и по ряду экспериментальных работ [Draeger et al., 2006; Грибов, 2015; 2016; 2017; Щербаков и др., 2017; Gribov et al., 2018]. Все эти эксперименты проводились на образцах, содержащих титаномагнетитовые (ТМ) зерна с Tc ≈ 200 °C при длительной выдержке при T_CRM ≈ (350$–$530) °C. В настоящей работе мы предприняли новый цикл компьютерного моделирования образования CRM- и TRM-экспериментов определения палеонапряженности по методам Телье и Вилсона$–$Буракова в системе взаимодействующих СПМ и ОД частиц с целью сравнения температурных интервалов, на которых наблюдается подобие спектров T_b при работе тем или иным методом, и выяснения возможных причин их возникновения.

МОДЕЛИРОВАНИЕ ОБРАЗОВАНИЯ ХИМИЧЕСКОЙ НАМАГНИЧЕННОСТИ

Подробно методика моделирования приобретения TRM и CRM описана в статье [Щербаков и др., 2017]. Ниже для удобства читателя кратко излагаются необходимые для дальнейшего изложения сведения.

Пусть x_i $–$ вероятность того, что магнитный момент i-й частицы m_i = M_sv_iµ_i ориентирован параллельно ее легкой оси l, направление которой выбрано так, чтобы произведение l · B было ≥ 0. Здесь: v_i $–$ объем частицы; µ_i $–$ единичный вектор, параллельный m_i, т. е. m_i =m_i µ_i, m_i = M_sv_i. Тогда функция x (t) подчиняется известному кинетическому уравнению:

$\frac{dx}{dt}=-\frac{x}{{\tau }_1}+\frac{1-x}{{\tau }_2},$ ${\tau }_{\mathrm{1,2}}\approx {\tau }_0\exp \left(E_{b\mathrm{1,}b2}/kT\right),$  (2)

где: τ₀ ≈ 10^-9 с $–$ характерный период попыток преодоления потенциального барьера; τ₁ и τ₂ $–$ времена релаксации для достижения равновесия в состояниях m↑↑l и m↑↓l, соответственно; t $–$ время; k $–$ постоянная Больцмана. В предположении доминирующей анизотропии формы и дипольного характера поля взаимодействия B_int, действующего на частицу с магнитным моментом m_i, имеем:

 

	$\left.\begin{array}{c}{E}_{b1}\\ E_{b2}\end{array}\right\}={\mu }_0\frac{N_d{M}_s^{2}v}{2}+{\mu }_0\mathbf{m}\cdot \left(\mathbf{B}+{\mathbf{B}}_{\mathrm{int}}\right),$
	$$B_{\mathrm{int}}(i)={\displaystyle \sum _{\begin{array}{l}j=1\\ i\ne j\end{array}}^{N_p}\left(\frac{m_j}{R_{ij}^3}-\frac{3R_{ij}\cdot (m_j{R}_{ij})}{R_{ij}^5}\right)},$$	(3)

где: M_s $–$ спонтанная намагниченность; R_ij $–$ расстояние между i-й и j-й частицами. Всюду в дальнейших расчетах принималось M_s(T) = M_s(0)j_s(T). Здесь j_s(T) = (1$–$T/T_c)^d $–$ относительная спонтанная намагниченность. Для магнетита d = 0.42 [Dunlop, Özdemir, 1997]. Из (2) следует, что равновесное значение x = τ₁/(τ₁ + τ₂), откуда равновесное значение намагниченности частицы есть 2x $–$1 = (τ₁$–$τ₂)/(τ₁ + τ₂) = tanh[m_i(B + B_int)/kT, в соответствии с общеизвестным выражением, получаемым из больцмановской статистики.

Если к системе дифференциальных и алгебраических уравнений (2) и (3) для совокупности N_p взаимодействующих частиц добавить еще и законы роста объема каждой индивидуальной частицы v_i = v_i(t), то эта система решает, в принципе, задачу нахождения намагниченности ансамбля СПМ зерен:

$$\begin{array}{c}M\left(t\right)=M_s\left(T\right){\displaystyle \sum }_{i=1}^{N_p}{v}_i\left(2x_i-1\right)\frac{\left(l\cdot \frac{B}{B}\right)}{L^3}=\\ =M_s\left(T\right)c{\displaystyle \sum }_{i=1}^{N_p}{v}_i\left(2x_i-1\right)\frac{\left(l\cdot \frac{B}{B}\right)}{{{\displaystyle \sum }}_{i=1}^{N_p}{v}_i},\end{array}$$

(4)

где с $–$ относительная объемная концентрация ферримагнетика. Отметим, что величина M(t) является полной намагниченностью, то есть она включает в себя вклад как СПМ, так и уже заблокированных зерен с τ₁ и τ₂ > t.

На практике, однако, для N_p >> 1 магнитостатически взаимодействующих частиц, такой подход представляется невозможной задачей даже для мощных вычислительных машин. По этой причине, как это часто происходит в таких сложных случаях, для моделирования процессов приобретения TRM и CRM в системе взаимодействующих одноосных ОД зерен был использован метод численного статистического моделирования Монте$–$Карло, в основе которого лежит использование случайных чисел для имитации различных вероятностных распределений.

При расчете рассматривался ансамбль N_p взаимодействующих зерен. Число N_p принималось максимально возможным при доступных вычислительных мощностях, конкретно, N_p = 5000 для взаимодействующих и N_p = 50 000 для невзаимодействующих частиц. Пространственные координаты частиц для каждой конфигурации задавались посредством случайного их выбора (с процедурой исключенного объема) в кубе размера L, выбранного таким образом, чтобы выполнялось равенство ${\displaystyle \sum }_{i=1}^{N_p}{v}_i=c{L}^3.$ Для каждой частицы случайным образом задавались пространственная ориентация легкой оси и направление магнитного момента вдоль этой оси. При c → 0 магнитостатическим взаимодействием, очевидно, можно пренебречь, и в этом варианте численные расчеты намного упрощаются просто принятием условия c = 0.

Моделирование методом Монте$–$Карло стохастических процессов, описываемых вероятностными формулами типа (2), заключается в генерации последовательности случайных событий, рассматриваемых начиная с t = 0 с некоторым временным интервалом Δt. Для этого сначала заметим, что на коротком интервале времени можно пренебречь изменением объема частицы. Тогда, полагая τ₁ и τ₂ постоянными величинами и решая (2) при условии x(0) = 1, получаем:

$x(t)=\exp \left(-\frac{1}{{\tau }_1}-\frac{1}{{\tau }_2}\right)+\frac{{\tau }_1}{{\tau }_1+{\tau }_2}\left(1-\exp \left(-\frac{1}{{\tau }_1}-\frac{1}{{\tau }_2}\right)\right).$       (5)

Затем программным генератором псевдослучайных чисел вырабатывается число r в интервале (0,1). Если r < x_i (Δt), момент не изменяется, но если r > x_i (Δt), m переворачивается. Продолжительность временнóго интервала Δt согласно оригинальной версии метода Монте$–$Карло (схема Метрополиса) должна быть меньше времен релаксации τ₁ и τ₂. Следуя этому, при каждом следующем испытании Монте$–$Карло полагалось Δt = τ_min/2, где τ_min $–$ минимальное время релаксации, рассчитанное по всем частицам. Для сокращения времени счета процедура Монте$–$Карло выполнялась для частиц ансамбля, времена релаксации которых превышали минимальное τ_min не более, чем в 100 раз, в противном случае частица считалась заблокированной, и ее магнитный момент не изменялся. Расчеты проводились при температуре кристаллизации Т_CRM = 400 °C с начальным размером частиц 34$–$37 нм и Т_CRM = 500 °C с начальным размером частиц 44$–$47 нм. Такой спектр размеров частиц был выбран для того, чтобы на лабораторных масштабах времени частицы были суперпарамагнитными с временем релаксации (0.1$–$10) секунд, при этом время релаксации рассчитывалось согласно (2), (3), где для простоты расчета принималось B_int = 0 (поскольку в нашем приближении mB/kT < 1, и в этом случае слагаемыми с B в формулах (2), (3) можно пренебречь). По ходу численного эксперимента размер каждой частицы случайно изменялся от цикла к циклу согласно соотношению d = (Δd₀/t₀) · Δt · random, где random $–$ случайное число от 0 до 1, t₀ = 86 400 с, Δd₀ = 15 нм. Таким образом, максимально возможный прирост размера частицы за одни сутки составлял 15 нм при «длительности» всего «эксперимента» t₁ от 3 до 10 суток, с тем, чтобы на конечном этапе время релаксации было не менее одного года. Заметим, что при T_CRM = 400 °C и N_d = 1 время релаксации 1 год имеет частица размером 43 нм, т. е. при начальных размерах 34$–$37 нм частице достаточно вырасти всего на 6$–$9 нм, чтобы стать надежно заблокированной на интервале времени в несколько суток. На рис. 1 представлены спектры объемов частиц после образования CRM и спектры деблокирующих температур.

 

Рис. 1. Спектры объемов частиц после образования CRM (а), (в) и спектры деблокирующих температур (б), (г): (а), (б) – для температуры кристаллизации T_CRM = 400 °C; (в), (г) – для T_CRM = 500 °C. За единицу объема принят объем частицы с Tb = T_CRM. Статистика приведена для случая невзаимодействующих частиц для 50 Монте–Карло-экспериментов (т. е. общее число частиц в данном случае N_p = 2 500 000).

 

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

M(t) = Ms(T)сj, $j=\left({\displaystyle \sum _{i=1}^{i=N_p}\left({\mu }_i\cdot B\right)}\cdot v_i/B\right)/{\displaystyle \sum _{i=1}^{i=N_p}{v}_i},$ (6)

где j $–$ относительная намагниченность, вычисляемая как сумма по i произведений косинуса угла между вектором магнитного момента i-й частицы m_i = M_sv_iμ_i и вектором внешнего поля В на объем i-й частицы vi, деленная на сумму объемов всех частиц ансамбля. Нетрудно видеть, что запись (4) есть дискретный аналог формулы (3), поскольку для ансамбля, содержащего бесконечное число частиц с легкой осью l_i (направленной параллельно либо антипараллельно μ_i) и объемом v_i, средняя величина случайной переменной (μ_i · B)/B как раз и есть функция (2x_i $–$ 1) (l · B/B).

Результаты численного моделирования процесса образования CRM представлены на рис. 2 и рис. 3 для T_CRM = 400 и 500 °C, соответственно, и различных значений с. Усреднение проведено по числу вариантов от 50 (для невзаимодействующих частиц) до 70$–$100 (при учете взаимодействия). Для определения среднего по всем экспериментам применялась интерполяция рядов данных B-сплайнами с использованием процедур математической библиотеки IMSL для Фортрана. Как уже подчеркивалось выше, M(t) включает в себя вклад СПМ зерен, что проиллюстрировано на рис. 2 и рис. 3 тонкими линиями, представляющими кривые относительного числа переворотов частиц в единицу времени. Отметим, что частота переворотов ведет себя зеркально по отношению к поведению j(t) и на финальном этапе расчетов, когда частицы блокируются, число переворотов, естественно, падает до нуля.

 

Рис. 2. Результаты численного моделирования методом Монте–Карло. Кинетика приобретения j[lg(t)] в ансамбле ОД зерен в поле B = 80 А/м при различных значениях относительной объемной концентрации ферримагнетика c для температуры кристаллизации T_CRM = 400 °C. Тонкими линиями нарисованы графики отношения числа переворотов частиц в секунду к общему числу частиц в зависимости от времени (правая ось ординат). Цифры при кривых дают объемную концентрация частиц после завершения их роста, в скобках приводится начальная концентрация.

 

Рис. 3. То же самое, что и на рис. 2, но для температуры кристаллизации T_CRM = 500 °C.

 

При взгляде на кривые намагничивания j(t) бросается в глаза их неожиданно немонотонное поведение со временем, хорошо наблюдаемое для невзаимодействующих частиц (c = 0%) и в какой-то мере проявляющее себя при начальной концентрации 0.05% (рис. 2 и рис. 3). Действительно, вначале мы видим их быстрый рост на интервале от 0 до ≈ 100 сек., затем стагнацию (медленный рост), продолжающуюся до (3$–$4) тыс с. Этот этап сменяется новым ростом j, совпадающим с падением числа переворотов частиц вплоть до их полной блокировки и заканчивающимся, вследствие этого, окончательной стабилизацией j на конечном этапе моделирования, что соответствует стадии формирования CRM.

Для понимания физических причин такого поведения нормированной намагниченности мы проанализировали численное решение уравнения (2) для x(t) для невзаимодействующих частиц (B_int = 0) при начальном условии x(0) = 1/2, предполагая, что рост размера частиц d со временем происходит по закону d = a₀ + bt. Пример такого расчета приведен на рис. 4, на котором ясно видны те же 4 стадии эволюции величины намагниченности со временем.

 

Рис. 4. Серая линия – зависимость намагниченности j(t), рассчитанной при следующих параметрах: b = 10^–4 нм/сек, a₀ = 34 нм. Черная линия маркирует величину CRM после блокировки магнитного момента частиц. На врезке показана зависимость размера образца (серая линия) и блокирующего размера (черная линия) от времени t. Пересечение этих кривых соответствует моменту блокировки магнитного момента частиц.

 

Очевидно, быстрый первоначальный рост j(t) обусловлен тем, что в начале расчета частицы полагались суперпарамагнитными (v << v_b, рис. 3) с нулевой суммарной исходной намагниченностью (x = 0.5), поэтому в присутствии ненулевого внешнего поля B намагниченность ансамбля 2 x $–$ 1 быстро достигла равновесной величины (τ₁$–$τ₂)/(τ₁ + τ₂). Заметим, однако, что длительность этого процесса (время релаксации) зависит от первоначального размера частицы и, как показывает расчет, если при a₀ = 38 нм τ₁ (и τ₂) ≈ 100 с, то уже при a₀ = 34 нм этот процесс занимает всего 0.1 с. Дальнейшая стагнация j (t) вплоть до t ≈ 10⁴ с обусловлена медленным ростом размеров зерен на этом интервале времени (врезка к рис. 4), так что равновесное значение намагниченности ансамбля (серая линия) практически не меняется со временем. Но при t ≈ 10⁴ с начинается значительный рост объема частиц, что отражается в росте равновесного значения намагниченности j (серая линия на рис. 4). Однако рост размера частиц, естественно, приводит к тому, что их объем сначала сравнивается (при t ≈ 10^4.75 с), а затем и превышает блокирующий объем v_b (врезка к рис. 4), в результате чего происходит фиксация CRM (черная линия на рис. 4). Выражение для блокирующего объема здесь было получено согласно известному выражению [Dunlop, Özdemir, 1997]:

$$v_b=\frac{2\ln \left(t/{\tau }_0\right)kT}{{\mu }_0{N}_d{M}_s^2}.$$

(7)

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЭКСПЕРИМЕНТОВ ОПРЕДЕЛЕНИЯ ПАЛЕОНАПРЯЖЕННОСТИ ПО МЕТОДАМ ТЕЛЬЕ$–$КОЭ И ВИЛСОНА$–$БУРАКОВА

Для оценки степени сходства температурной стабильности CRM и TRM в ансамблях взаимодействующих частиц нами проведено сравнение их спектров блокирующих температур путем численного моделирования экспериментов определения палеонапряженности методами Вилсона$–$Буракова и Телье$–$Коэ. Для этого вначале моделировалось приобретение CRM на ансамбле ОД частиц с последующим его «охлаждением» в поле 100 мТ от T_CRM до комнатной температуры T_r (образование pTRM). Затем полученная остаточная намагниченность CRM + pTRM «терморазмагничивалась» при «нагреве» ансамбля до температуры Кюри магнетита T_c = 585 °C. Образование и «терморазмагничивание» TRM имитировалось путем «охлаждения» ансамбля в поле 100 мТ от T_c до T_r с последующим «нагревом» до T_c. Сравнение этих термокривых и составляет суть метода Вилсона$–$Буракова.

Моделирование экспериментов определения палеонапряженности по методу Телье$–$Коэ на TRM и CRM + pTRM проводилось по схеме, описанной в работе [Щербаков и др., 2017]. Каждый численный эксперимент включал в себя 9 парных «нагревов$–$охлаждений», по результатам которых строились диаграммы Араи$–$Нагата. Результаты численных экспериментов показаны на рис. 5$–$рис. 9 при значениях относительной объемной концентрации с = 0, 0.2, 0.32, 0.61, 1.27% (T_CRM = 400 °C) и 0%, 0.11%, 0.25%, 0.48%, 1.13% (T_CRM = 500 °C). Диаграммы (а) и (г) на всех этих рисунках дают стандартное представление метода Вилсона$–$Буракова путем сопоставления кривых терморазмагничивания CRM (T) (нижняя кривая) и TRM (T) (верхняя кривая), в то время как диаграмма {CRM, TRM} на рис. 6а, 6б и 6д построена по аналогии с методикой Телье путем построения непрерывной зависимости CRM (TRM), где величины CRM и TRM отвечают текущей температуре Т. Такая методика представления результатов метода Вилсона$–$Буракова была недавно предложена в работе [Muxworthy et al., 2011], и она имеет то преимущество, что на ней наглядно видны участки с подобными спектрами T_b, отвечающие прямолинейной зависимости CRM (TRM). Поскольку одной из задач данного исследования является как раз выяснение степени подобия спектров CRM и TRM, такие диаграммы были построены для всех просчитанных случаев с целью обнаружения и фиксации таких прямолинейных участков.

 

Рис. 5. Результаты численного моделирования для невзаимодействующих частиц (с = 0 %): (а), (б), (в) – T_CRM = 400 °C; (г), (д), (е) – T_CRM = 500 °C; (а) и (д) – кривые терморазмагничивания CRM (нижняя кривая) и TRM (верхняя кривая); (б) и (д) – метод Вилсона; (в) и (е) – метод Телье. Для построения графиков использованы усредненные данные по 50 Монте–Карло-экспериментам. Все графики приведены только на интервале от температур выше T_CRM до T_c, поскольку ниже T_CRM как CRM, так и TRM не изменяются по величине. Величина k, приведенная на графиках, есть k = B_anc/B. Ошибки расчета величин CRM и TRM на диаграмме Араи–Нагата укладываются в размер квадрата, представляющего соответствующую точку.

 

Рис. 6. То же, что и на рис. 5, но для взаимодействующих частиц (с = 0.2 (0.05)%, T_CRM = 400 °C и с = 0.11 (0.05) %, T_CRM = 500 °C). Отрезки прямых вокруг точек на диаграмме Араи–Нагата представляют собой ошибки расчета соответствующих точек.

 

Как видно из рис. 5$–$рис. 9, на всех диаграммах {CRM, TRM}, действительно, присутствуют прямолинейные участки (обозначены пунктирными линиями) в области температур, нижняя граница которой лежит заметно выше температуры образования CRM, так что область подобия спектров сдвинута в сторону T_c и занимает достаточно узкий температурный интервал (T₁, T₂) шириной от 12 до 60 °C. На диаграмме {CRM, TRM} это составляет примерно половину ее полной длины, что позволяет выполнить формальное определение величины древнего поля (соответствующие величины приведены на диаграммах). Вместе с тем, на всем остальном температурном интервале кривая CRM (TRM) оказывается существенно выгнутой, что говорит об отсутствии здесь подобия спектров.

 

Рис. 7. То же, что и на рис. 6, но при с = 0.32 (0.1)%, T_CRM = 400 °C и с = 0.25 (0.1)%, T_CRM = 500 °C.

 

На рисунках (в) и (е) представлены расчетные диаграммы Араи$–$Нагата {CRM, pTRM} и выполнено определение Banc на том же интервале (T₁, T₂), где диаграммы Вилсона$–$Буракова демонстрируют подобие спектров. Как видно, во всех случаях наблюдается хорошее соответствие между этими методами как по величине коэффициента k, так и по интервалу прямолинейности диаграмм (T₁, T₂).

 

Рис. 8. То же, что и на рис. 6, но при с = 0.61 (0.2)%, T_CRM = 400 °C и с = 0.48 (0.2)%, T_CRM = 500 °C.

 

Рис. 9. То же, что и на рис. 6, но при с = 1.27 (0.5) %, T_CRM = 400 °C и с = 1.13 (0.5) %, T_CRM = 500 °C.

 

Рис. 10. Зависимость коэффициента линейной аппроксимации k, рассчитанный по диаграммам Араи–Нагата (квадраты) и Вилсона–Буракова, (треугольники). Пóлые фигуры – T_CRM = 400 °C, заштрихованные – T_CRM = 500 °C.

 

Как видно из рис. 10, коэффициент определения k изменяется в пределах от 0.19 до 0.51 для различных методов и температур реакции, при этом не обнаруживается какой-либо однозначной зависимости k от концентрации c. В то время как при c < 0.4 % видна тенденция к уменьшению k при росте c, при больших c мы видим даже увеличение k, которое, скорее всего, следует связать с ошибками расчета, поскольку величина полей взаимодействия при таких концентрациях становится слишком большой и метод Монте$–$Карло начинает работать неустойчиво, результаты от цикла к циклу могут различаться очень сильно, так что в этом случае необходимо использование очень мощных вычислительных устройств, чего не было в нашем распоряжении.

ДИСКУССИЯ

Как видно из рис. 5$–$рис. 9, на всех диаграммах Араи$–$Нагата и Вилсона$–$Буракова есть область подобия спектров, по которой можно выполнить формальное определение Banc, при этом его величина оказывается в 2$–$5 раз меньше реальной величины поля, используемого при моделировании (рис. 10). Этот вывод соответствует как имеющимся экспериментальным данным [Draeger et al., 2006; Грибов и др., 2015; 2016; 2017; Щербаков и др., 2017; Gribov et al., 2018], так и нашим предыдущим расчетам. Новизна представленных расчетов заключается в определении области подобия спектров, которая, как оказалось, занимает далеко не весь интервал (T_CRM, T_c), а лишь достаточно узкий температурный интервал (T₁, T₂) шириной от 12 до 60 °C в окрестности T_c. На остальном же интервале CRM разрушается быстрее, чем TRM, что приводит к заметной вогнутости диаграмм Араи$–$Нагата и Вилсона$–$Буракова.

Из диаграмм, приведенных на рис. 5$–$рис. 9, можно также сделать вывод о том, что результаты, полученные для различных концентраций частиц (отражающих степень их магнитостатического взаимодействия), показывают их большое сходство, в смысле общего поведения, диаграмм Араи$–$Нагата и Вилсона$–$Буракова: на всех этих диаграммах отмечаются заметная их кривизна в начале этих диаграмм на интервале (T_CRM, T₁) и близкое к линейному поведение при увеличении температур нагрева до T₁ и выше. При этом следует отметить, что на диаграммах для невзаимодействующих частиц некоторая кривизна заметна даже и на условно линейном интервале (T₁, T₂).

Таким образом, наши расчеты согласуются с гипотезой [Smirnov, Tarduno, 2005] о том, что подобие спектров может быть обусловлено узостью интервала блокирующих температур TRM и CRM, но с той поправкой, что оно наблюдается лишь на части полного интервала спектра T_b CRM. Вопрос состоит в том, что же именно обеспечивает это подобие и какова его реальная степень, учитывая, что на интервале (T₁, T₂) сосредоточена все же значительная часть полной интенсивности CRM.

Для расчета кривизны диаграмм Араи$–$Нагата и Вилсона$–$Буракова хотя бы в случае невзаимодействующих частиц, когда вычисления можно провести в аналитическом виде, нами была рассчитана термокривая CRM (T), созданная по механизму Хайга роста зерен с последующим расчетом термокривой TRM (T), созданной на этом же ансамбле частиц, по следующей схеме. Восприимчивость ансамбля неориентированных невзаимодействующих ОД СПМ частиц с объемом v есть:

$$\chi (T)=\frac{{[m(T)]}^{2}n}{3kT}=\frac{{[M_s(T)]}^{2}vc}{3kT}.$$

(8)

Принимая в (7) для определенности характерное время опыта t ≈ 100 с, получим стандартную оценку ln (2 t/t₀) = 50. Заменяя теперь в (8) v на v_b из (7) и учитывая рост M_s при охлаждении до T₀, найдем $\text{TRM}=\frac{50cB{j}_s(T_0)}{N_d{j}_s(T_b)}$ [Shcherbakov et al., 1996], где j_s(T) = M_s(T)/M_s(T₀) $–$ приведенная спонтанная намагниченность. Распространяя это выражение на ансамбль частиц с функцией распределения по объемам f(v), найдем для кривой терморазмагничивания:

$$\text{TRM}(T)=\frac{50cB{j}_s(T_0)}{3N_d}{\displaystyle \underset{v_b(T)}{\overset{\infty }{\int }}\frac{f(v)}{j_s(T_b)}dv}.$$

(9)

Заменяя здесь T_b на T_CRM, получим аналогичное выражение для кривой терморазмагничивания CRM (T):

$$\text{CRM}(T)=\frac{50cB{j}_s(T_0)}{3N_d{j}_s(T_{\text{CRM}})}{\displaystyle \underset{v_b(T)}{\overset{\infty }{\int }}f(v)dv}.$$

(10)

Из (9) и (10) следует, что производная диаграммы Араи$–$Нагата:

$$\begin{array}{c}\frac{d\text{CRM}}{dp\text{TRM}}=\frac{d\text{CRM}(T)/dT}{-d\text{TRM}(T)/dT}=\\ =\frac{-j_s(T)}{j_s(T_{\text{CRM}})}=\frac{-{(1-T/T_c)}^d}{{(1-T_0/T_c)}^d},\end{array}$$

(11)

где pTRM (T) = TRM (T_r) $–$ TRM (T). Этот результат, разумеется, находится в полном соответствии с формулой (1), записанной для отдельного зерна с блокирующей температурой Tb. Любопытно, что производная (11), представленная как функция T/T_CRM имеет универсальную форму, не зависящую от T_CRM и спектра размеров ансамбля. Если же перейти к интересующей нас производной диаграммы Араи$–$Нагата, то она в параметрическом виде записывается как пара {pTRM(T), dCRM(T)/dTRM(T)} и в этом случае форма этой кривой уже зависит от формы спектра T_b (рис. 11), но по-прежнему не зависит от T_CRM.

 

Рис. 11. (а) – Функции распределения объемов частиц f(v), использованные при расчете графиков (б), (в), (г), (д) параметры логнормального распределения: 1 – σ = 0.5, μ = 0.5; 2 – σ = 0.5, μ = 1.6. За единицу объема принят объем частицы с блокирующей температурой T_CRM – v_b(T_CRM); (б), (г) – производные dCRM/dpTRM как функции pTRM (x)/TRM(xr); (в), (д) – диаграммы Араи–Нагата. Жирной черной линией отмечены плато на кривых dCRM/dpTRM и участки квазилинейного поведения этих диаграмм. Нумерация кривых на рисунках (б)–(д) соответствует их нумерации на рис. (а).

 

При полном подобии CRM и TRM производная dCRM/dpTRM должна быть константой, что, как видно из формулы (11) и рис. 11, не выполняется нигде. Тем не менее, на всех графиках производной на рис. 11 можно выделить температурный интервал, близкий к Tc, на котором эта производная изменяется не более чем в полтора раза, и соответствующие участки на диаграммах Араи$–$Нагата и Вилсона$–$Буракова можно принять за линейные и попытаться определить на них «древнее поле», как это проиллюстрировано на рис. 11. Появление такого квазилинейного участка объясняется достаточно просто: как видно из формулы (11) и рис. 11б, производная dCRM/dpTRM монотонно изменяется от 1 при T = T_CRM до 0 при T = T_с, при этом необходимо отметить очевидную тенденцию: чем шире спектр размеров частиц в сторону частиц бóльшего размера (соответственно, бóльших T_b), тем резче спад dCRM/dpTRM при T, близких к T_CRM, и тем более плоскими и низкими по величине эти кривые становятся при приближении к T_CRM. Именно наличие такого плато на кривых dCRM/dpTRM и порождает квазилинейные участки на диаграммах Араи$–$Нагата в случае невзаимодействующих частиц. При этом ясно, что чем ниже значения производной, тем больше занижение величины древнего поля, определяемого по квазилинейному участку. Иными словами, чем больше спектр T_b удален от T_CRM и сосредоточен вблизи T_c, тем больше должно быть занижение Banc.

Переходя теперь к обсуждению результатов численных расчетов, прежде всего заметим, что диаграммы Араи$–$Нагата и Вилсона$–$Буракова, приведенные на рис. 5$–$рис. 9 для различных концентраций частиц (отражающих степень их магнитостатического взаимодействия), указывают на их большое сходство в смысле общего поведения этих диаграмм $–$ на всех диаграммах отмечается заметная кривизна на интервале (T_CRM, T₁) и, близкое к линейному, поведение на интервале (T₁, T₂). Но при этом, справедливости ради, следует отметить, что на диаграммах для невзаимодействующих частиц некоторая кривизна заметна даже и на условно линейном интервале (T₁, T₂). Здесь уместно вспомнить о гипотезе о том, что именно взаимодействие приводит к сходству спектров T_b CRM и TRM [Draeger et al., 2006]. К сожалению, из наших расчетов нельзя сделать определенное заключение о справедливости этой гипотезы, поскольку, с одной стороны, разница в форме диаграмм Араи$–$Нагата и Вилсона$–$Буракова для взаимодействующих и невзаимодействующих частиц на рис. 5$–$рис. 9 невелика, а с другой, точность наших расчетов (ошибки определения) для взаимодействующих частиц не позволяет с уверенностью оценить степень возможной кривизны этих диаграмм на всем температурном интервале.

Из вышесказанного следует, что предположение [Draeger et al., 2006] о том, что подобие спектров T_b для TRM и CRM обеспечивается магнитостатическим взаимодействием между ОД-зернами, не вполне верно: достаточная степень подобия может наблюдаться и для ансамбля невзаимодействующих частиц в силу указанных выше причин. Вместе с тем, остается открытым вопрос о том, может ли такое взаимодействие улучшить степень подобия. Исходя из наших расчетов, можно предположить, что дело именно так и обстоит. Действительно, из рис. 5 видно, что диаграммы Араи$–$Нагата и Вилсона$–$Буракова для невзаимодействующих частиц демонстрируют хоть и небольшую, но ощутимую кривизну (в соответствии с теоретическими расчетами этого раздела) даже и на тех участках, которые мы принимаем за квазилинейные. В то же время на соответствующих диаграммах для взаимодействующих частиц кривизна на таких участках практически отсутствует. Однако ошибки расчетов, показанные вертикальными и горизонтальными отрезками на этих диаграммах, слишком велики, чтобы с уверенностью говорить об отсутствии этой кривизны, так что утверждать с определенностью, что именно магнитостатическое взаимодействие приводит к линейности диаграммы Араи$–$Нагата и Вилсона$–$Буракова, пока нельзя $–$ здесь необходимы длительные и более тщательные расчеты с использованием мощных компьютеров.

Перейдем теперь к сравнению моделированных и экспериментальных данных, насколько это возможно, исходя из имеющегося материала. Наиболее подробно экспериментальный материал по лабораторному созданию и исследованию свойств CRM, включая методику Телье, представлен в работах [Draeger et al., 2006; Грибов, 2015; 2016; 2017; Щербаков и др., 2017; Gribov et al., 2018]. Из этих экспериментов следует, что, как правило, (квази)линейный участок на диаграммах Араи$–$Нагата для CRM, присутствует всегда и занимает почти весь интервал (T_CRM, T_c) (при T, близком к T_c становятся заметными минералогические изменения, искажающие результат). Хотя справедливости ради надо отметить, что в окрестности T_CRM обычно присутствует чуть ощутимая степень кривизны, и строго линейной диаграмма Араи$–$Нагата становится лишь при T ≈ (475$–$500) °C (или позже, если T_CRM = 500 или 530 °C). Таким образом, эти данные не согласуются с гипотезой [Smirnov, Tarduno, 2005], поскольку интервал подобия спектров оказывается отнюдь не узким. Но они не согласуются в полной мере и с нашими численными расчетами, поскольку заметной кривизны на экспериментальных диаграммах не наблюдается даже и в самом их начале.

По-видимому, дело заключается в том, что лабораторные эксперименты имеют дело далеко не с идеализированным случаем образования CRM по механизму Хайга, как это предполагалось в теоретических расчетах. Действительно, на самом деле в этих экспериментах приобретение CRM в образцах происходило еще на стадии однофазного окисления [Gribov et al., 2018], то есть при образовании CRM были задействованы оба механизма, указанные во Введении, и поэтому прямое сравнение теории и эксперимента здесь неправомерно.

ВЫВОДЫ

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

Результаты моделирования экспериментов определения палеонапряженности по методике Телье показали, что на всех диаграммах Араи$–$Нагата {CRM, TRM} присутствует область приближенного подобия спектров, нижняя граница которой лежит заметно выше T_CRM, а верхняя граница приближается к T_c. Эта область занимает довольно узкий температурный интервал (T₁, T₂) шириной от 12 до 60 °C, составляя при этом примерно половину ее полной длины диаграммы {CRM, TRM}. Аналитический разбор случая невзаимодействующих частиц (очень малая концентрация зерен) показал, что (квази)линейность в этом случае объясняется наличием плато на кривых производной dCRM/dpTRM, возникающей при узком спектре блокирующих температур.

По результатам численных экспериментов можно предположить, что при концентрации частиц, превышающей 0.2%, магнитостатическое взаимодействие приводит к практически полной линейности диаграммы Араи$–$Нагата на значительном интервале полного спектра блокирующих температур CRM. Вместе с тем, на остальном температурном интервале как диаграммы Араи$–$Нагата, так и кривые зависимости CRM (TRM) оказываются существенно вогнутыми, что говорит об отсутствии подобия спектров.

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

Для дальнейшего прогресса в понимании механизма возникновения CRM и корректного описания ее свойств, необходимо: а) выполнение экспериментов строго по модели Хайга, то есть путем выращивания ферримагнитных зародышей до ОД состояния в немагнитной матрице; б) проведение численного моделирования с использованием мощных компьютеров; в) развитие теории образования CRM по второму механизму, то есть в результате роста температуры Кюри титаномагнетитовых зерен вследствие их окисления.

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

Работа выполнена в рамках госзадания ИФЗ РАН $№$ 0144-2014-00117 и при поддержке гранта РФФИ $№$ 17-05-00259.

About the authors

V. P. Shcherbakov

Borok Geophysical Observatory, Schmidt Institute of Physics of the Earth, Russian Academy of Sciences

Author for correspondence.
Email: shcherbakovv@list.ru
Russian Federation, Borok

N. K. Sycheva

Borok Geophysical Observatory, Schmidt Institute of Physics of the Earth, Russian Academy of Sciences

Email: shcherbakovv@list.ru
Russian Federation, Borok

References

Supplementary files