Микромагнетизм одноосных ферромагнетиков

Размещено на http://www.allbest.ru

Размещено на http://www.allbest.ru

Введение

Актуальность темы. Началом теории микромагнетизма считается опубликованная в 1935 году работа Ландау и Лифшица [1]. С тех пор микромагнитный подход получил широкое распространение как для получения фундаментальных представлений о закономерностях формирования и эволюции магнитных структур в ферромагнетиках, так и для оценки потенциальных возможностей их практического использования в различных областях техники. Рассматриваемые в теории микромагнетизма магнитные материалы обычно называют ферромагнетиками, хотя микромагнитный подход применим и для других магнетиков, магнитные моменты в которых ориентированы упорядоченным образом, а усредненный магнитный момент не равен нулю и не зависит (слабо зависит) от внешнего поля [2].

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

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

Рис. 1. Данные из работы [3] о росте плотности записи с 1980 года в применяемых (традиционных) методах магнитной записи

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

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

- анализ влияния обменного и магнитостатического взаимодействий на характер термического намагничивания многослойной стохастической системы;

- выяснение механизма термического намагничивания в бесконечно длинных монокристаллах;

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

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

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

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

Научная новизна работы заключается в том, что впервые:

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

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

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

На основе разработанных методов получены следующие новые результаты:

- установлена возможная причина появления «седловых точек» при минимизации функционала свободной энергии и предложен метод их устранения;

- установлено существование при некоторых условиях доменных структур Ландау в бесконечно длинных монокристаллических призмах;

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

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

- показан механизм появления эффекта термического намагничивания в бесконечно длинном монокристаллическом стержне высокоанизотропного магнетика в некотором диапазоне поперечных размеров;

- установлено, что при симметричных внутри полубесконечного монокристалла доменных структурах, на торцевой поверхности в широком диапазоне изменения анизотропии возможно устойчивое существование как симметричных (или близких к симметричным), так и асимметричных доменных структур;

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

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

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

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

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

На защиту выносятся:

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

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

- критерий ориентации однородного поля намагниченности в бесконечно длинном монокристалле;

- механизмы термического намагничивания в многослойной стохастической системе и в бесконечно длинном монокристалле;

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

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

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

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

Апробация работы. Материалы диссертации докладывались и обсуждались на следующих конференциях: на 1-ой Всероссийской научно - практической конференции «Измерения, автоматизация и моделирование в промышленности и научных исследованиях» (Бийск, 2000); на XIII Международной конференции по постоянным магнитам (Суздаль,2000); на 2-ой Всероссийской научно - практической конференции «Информационные технологии в экономике, науке и образовании» (Бийск, 2001); на 4-й Всероссийской научной конференции «Краевые задачи и математическое моделирование» (Новокузнецк, 2001); на 2-ой Всероссийской научно - практической конференции «Измерения, автоматизация и моделирование в промышленности и научных исследованиях» (Бийск, 2001); на 3-ей Всероссийской научно - практической конференции «Измерения, автоматизация и моделирование в промышленности и научных исследованиях» (Бийск, 2002); на Всероссийской школе-семинаре «Магнитная анизотропия и гистерезисные свойства редкоземельных сплавов» (Тверь, 2002); на 7-й международной школе-семинаре «Эволюция дефектных структур в конденсированных средах» (Барнаул, 2003); на 2-ой Всероссийской научно-технической конференции «Физические свойства металлов и сплавов» (Екатеринбург, 2003); на 5-ой Всероссийской научно-технической конференции «Измерения, автоматизация и моделирование в промышленности и научных исследованиях» (Бийск, 2004); на 4-ой Всероссийской научно-технической конференции «Информационные технологии в экономике, науке и образовании» (Бийск, 2004); на VIII Международной школе-семинаре «Эволюция дефектных структур в конденсированных средах» (Барнаул, 2005); на Всероссийской научно-практической конференции «Фундаментальные науки и образование» (Бийск, 2006); на Российско-Японском семинаре «Магнитные явления в физикохимии молекулярных систем» (Оренбург, 2006); на 2-ом Российско-Японском семинаре «Магнитные явления в физикохимии молекулярных систем» (Оренбург, 2007); на II Всероссийской научно-практической конференции «Фундаментальные науки и образование» (Бийск, 2008); на Всероссийской научно-практической конференции «Многопрофильный университет как региональный центр образования и науки» (Оренбург, 2009); на Международной конференции «Фотоника молекулярных наноструктур» (Оренбург, 2009).

Публикации. По материалам диссертации опубликовано 49 печатных работ, в том числе 15 в изданиях, рекомендованных перечнем ВАК РФ. Получено 3 свидетельства о государственной регистрации программ для ЭВМ.

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

1. Расчеты полей намагниченности методами численного моделирования

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

2. Численное моделирование полей намагниченности в монокристаллах бесконечной длины

Представлены следующие результаты:

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

- расчеты полей намагниченности в бесконечно длинных монокристаллических призмах Co, Nd2Fe14B и Ni80Fe20 с квадратными поперечными сечениями разных размеров;

- критерий ориентации однородного поля намагниченности в бесконечно длинном монокристалле квадратного сечения;

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

Широко распространенные методы вычисления распределения намагниченности в микромагнитной системе основаны на минимизации дискретного аналога функционала свободной энергии. Для этой цели обычно используются градиентные методы, которые требуют вычисления частных производных по независимым переменным. В качестве переменных выбираются или компоненты вектора намагниченности М, или углы, определяющие ориентацию М. Если потенциал размагничивающего поля вычисляется с помощью уравнения Пуассона, то сложность вычисления частных производных вызвана неявной зависимостью функционала от поля М через потенциал. При минимизации функционала без учёта неявной зависимости встречаются точки (не являющиеся точками локальных минимумов), в которых любое смещение в направлении, противоположном «градиенту», приводит к росту функционала. Можно предположить, что такие особенности в ряде работ называют "седловыми точками" или "точками бифуркации", т.е. рассматривают их как точки лабильного равновесия. Однако они могут быть исключены, если при вычислении частных производных принять во внимание неявную зависимость функционала от поля M. Конструкция расчетного метода рассматривается на примере задачи о распределении намагниченности в монокристаллической призме бесконечной высоты с квадратным поперечным сечением . При выборе в качестве характерного линейного размера стороны квадрата и ориентации оси легкого намагничивания (ОЛН) по координатной оси Ox (рис. 1) функционал свободной энергии [2] можно представить в следующей безразмерной форме:

(1)

где в подынтегральном выражении записана сумма плотностей энергии обмена, энергии анизотропии, энергии системы во внешнем и собственном (размагничивающем) магнитных полях; - единичный вектор; - константа обмена; D = - безразмерное поперечное сечение; - константы анизотропии; - напряженности внешнего и собственного полей. Безразмерное поле можно выразить через потенциал U : .



Рис. 2

Потенциал находится из решения задачи:

, (2)

на границе кристалла, (3)

n - единичная внешняя нормаль к поверхности кристалла. Здесь (в формуле (3)) и далее использовано определение производной по направлению для функции , действующей из линейного нормированного пространства в линейное нормированное пространство : , где - точка, в которой вычисляется производная; - единичный вектор ; ; - число. В случае формулы (3) является пространством геометрических векторов, а является числовым полем. Предполагается, что поле является двумерным, т.е. выполняется условие .

Задача, состоящая в поиске локального минимума функционала (1) при условиях (2) - (3), математически некорректна, поэтому необходимо дополнить постановку правилом выбора локального минимума. Можно предположить, что эволюция системы происходит по траектории наискорейшего понижения функционала. При численных методах решения задачи данное предположение приводит к применению градиентного спуска для дискретного аналога функционала (1) - . - функция конечного числа переменных, в качестве которых в настоящей работе выбраны углы и (рис. 3). Для того, чтобы найти , необходимо вычисление частных производных по выбранным переменным, которое можно провести по обычным правилам дифференцирования для дискретных аналогов всех членов выражения (1) кроме последнего, сложность которого в неявной зависимости потенциала U от векторного поля m, задаваемой уравнениями (2) - (3). Рассмотрим возможный способ решения этой проблемы.

Рис. 3

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

,

где - дискретный аналог

. (4)

Здесь - индексы сеточных узлов в области D; - компоненты вектора m в узле ; - линейные функции, равные произведению дискретного аналога производной от потенциала U в узле (по координате x или y соответственно) на долю ячейки d в области D, приходящуюся на узел. Рассматриваются три типа узлов, показанных на рис. 1. Если узел внутренний (A), то обе производные аппроксимируются центральными разностными отношениями и d=1. Если узел расположен на регулярном (не угловом) участке границы (B), то производная по линии нормали односторонняя, ориентированная во внутрь D, а производная вдоль границы - центральная, d=0.5. В угловых точках (C) обе производные односторонние, ориентированные по границе D, а d=0.25. Так, например, для ячеек A, B, C на рис. 1 соответственно:

Сеточные компоненты вектора m можно выразить через углы (рис. 2):

. (5)

После такой замены (и, в частности, ) становятся функциями 2 (n+1)2 независимых переменных i,j и i,j.

Для реализации градиентного спуска необходимо найти частные производные . Пусть (i0, j0) - любой фиксированный узел сетки в области D. Из (4) и (5) следует, что:

(6)

(7)

Здесь учтено, что изменение углов и в одной точке влияет на поле потенциала U во всей области D. Получим уравнение для вычисления сеточной функции , которая входит в правые части выражений (4), (6) и (7).

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

, (8)

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

. (9)

Рассмотрим объёмную ячейку (рис. 4) с центром в точке , где - любой узел сетки внутри, на границе или вне области D, - шаг сетки в размерных единицах (в см). Интегрирование (9) по объему ячейки и последующая замена по теореме Гаусса объемных интегралов поверхностными приводит к формуле:

, (10)

где - поверхность выделенной ячейки, - единичная внешняя нормаль на поверхности . Из предположений и (что эквивалентно и ) следует, что интеграл в левой части (10) и по верхней, и по нижней поверхностям пластинки равен нулю, а в правой части равна нулю сумма интегралов по верхней и нижней поверхностям. Таким образом, при вычислении (10) достаточно ограничиться боковой поверхностью ячейки. Пусть центр ячейки - любой сеточный узел внутри области D. Замена интегралов в (10) дискретными аналогами дает сеточное уравнение:

,

Где:

. (11)

Рис. 4

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

,

где .

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

, (12)

где - внутри области D и - вне области D.

В регулярных точках границы D:

В угловых точках:

Здесь использованы обозначения:

где - любая сеточная функция.

Формула (12) аппроксимирует уравнение (2) с порядком внутри и снаружи области D, а условие (3) с порядком в регулярных точках границы. В соответствие с теорией потенциала [4] сеточную функцию можно рассматривать в качестве дискретного аналога объемной плотности источника потенциала. Причем, если внутри области D эта функция является обычной аппроксимацией объемного источника , то на ее границах - эффективный источник, учитывающий поверхностный источник с плотностью .

Дифференцируя (12) по переменным , получаем уравнения для вычисления :

. (13)

В выражениях (13) и всюду, кроме точек, упомянутых ниже.

Так, если , то ненулевые и следующие:

Если i0=n, 0 <j0<n, то:

Если i0=0, 0 <j0<n, то:

Если 0 <i0<n, j0=n, то:

Если 0 <i0<n, j0=0, то:

Если i0=0, j0=n, то:

Если i0=n, j0=0, то:

Если i0=0, j0=0, то:

Если i0=n, j0=n, то:

Здесь:

Решение уравнения (12) в области D сводится к задаче Дирихле. Для этого снаружи области D, на некотором расстоянии от границы D, выделяется контур R (R - граница квадрата, в центре которого расположена область D), в узловых точках которого (iR,jR) потенциал U вычисляется по формуле:

, (14)

где xi, yj и , - соответственно координаты узлов в области D и на выделенном контуре. Формула (14) является дискретным аналогом логарифмического потенциала [5], создаваемого источником , записанным в уравнении (12). При известных значениях потенциала на границе контура , внутри контура, содержащего область D, потенциал находится решением уравнения (12) методом последовательной верхней релаксации. Таким же методом можно решать уравнения (13), но вычисление производных значительно упрощается, если уравнения (13) рассматривать на некотором квадрате R* с центром в точке (рис. 2). Предполагается, что область R* смещается относительно области D при перемещении точки в области D. Источники полей имеют локальный характер. Поэтому при достаточно большой области R* на её границе можно задать условия:

. (15)

Точное решение уравнений (13) с условиями (15) известно [6] и в векторной форме имеет вид:

.

(16)

Здесь U' и f - векторы левых и правых частей уравнений (13), P2 - количество ячеек в квадрате R* (рис. 1), ukl - полная ортогональная система собственных функций дискретного аналога оператора Лапласа, обращающихся в 0 на границе R*. Верхний двойной индекс обозначает имя собственной функции. kl - собственное число функции ukl. - норма конечномерного гильбертова пространства с обычным определением скалярного произведения. В компонентах решение (16) выглядит так:

для 0 <i0<n, 0 <j0<n: ;

для i0=n, 0 <j0 <n: ;

для 0 <i0 <n, j0=n: ;

для i0=n, j0=n: ;

для i0=0, j0=n: ;

для i0=0, 0 <j0 <n: ; (17)

для 0 <i0 <n, j0=0: ;

для i0=0, j0=0: ;

для i0=n, j0=0: .

(18)

Здесь - индексы узлов в квадрате R*, , . В случае первого уравнения (13) - компоненты на сетке квадрата R*, а . Для второго уравнения (13) соответственно - компоненты , а .

Формулы (18) зависят только от числа p, задающего размер области R* в единицах h, поэтому коэффициенты (18) можно вычислить один раз при достаточно большом p. Затем они могут многократно использоваться в расчетах компонент U' по простым формулам (17).

Сравнительные расчёты проводились минимизацией функционала с учётом и без учёта неявной зависимости от поля m. В последнем случае суммы в правых частях (6)-(7) равны нулю, что позволяет избежать вычисления полей . Оказалось, что в ряде случаев предлагаемый метод не дает существенных отличий от моделирования без учёта неявной зависимости. Однако были обнаружены примеры, когда результаты моделирования качественно различались. На рис. 5 показаны расчёты распределения намагниченности, в которых минимизация функционала проводилась с учётом и без учёта неявной зависимости от векторного поля m. В вычислениях использовались значения магнитных параметров, соответствующих : эрг/см, Гс, , . Размер монокристалла нм. Поле m в начальном состоянии ориентировано по координатной оси x.



Рис. 5. Распределение намагниченности в монокристалле Со (сверху вниз: общий вид, фрагмент, схема) рассчитанное: а) минимизацией функционала без учёта неявной зависимости от поля намагниченности; б) с учётом неявной зависимости. Показана проекция поля m на поперечное сечение монокристалла

При минимизации с без учёта неявной зависимости модуль градиента уменьшается от начального значения более чем в 105 раз до величины, близкой к 0, и дальнейшее смещение в направлении приводит к росту . Векторы намагниченности отклоняются от начального положения только в окрестностях угловых точек (рис. 5а). Компонента намагниченности монокристалла уменьшается незначительно (от 1 до 0.98), что создает иллюзию однодоменного состояния. «Включение» неявной зависимости позволяет продолжать спуск, приводя к многодоменному состоянию c (рис. 5б). С математической точки зрения состояние на рис. 5а соответствует точке локального минимума функционала по некоторому направлению в пространстве независимых переменных . В общем случае эта точка не является локальным минимумом функционала , что и имело место в рассмотренном примере.

Расчёты намагниченности во внешнем поле 2130 Э, ориентированном против оси х, приводят к одинаковым результатам для обоих методов минимизации. Распределение намагниченности в этом случае близко к показанному на рис. 5б и отличается тем, что ориентированные по полю домены вырастают за счёт ориентированных против поля.

Работу Ландау и Лифшица [1], в которой доменная структура бесконечного монокристалла прямоугольного сечения представлена в виде, показанном на рис. 6, часто называют пионерской, основополагающей и т.д. в теории микромагнетизма. Однако данный вид доменной структуры получен не в результате расчетов, а построен из общих физических соображений и являлся исходным предположением для вычисления количественных характеристик структуры. Физические соображения, по которым построена структура Ландау, состоят в следующем. Предполагается, что данная структура должна обеспечивать нулевую магнитостатическую энергию монокристалла ввиду отсутствия поверхностных и внутренних источников собственного поля. Отсутствие поверхностных источников обеспечивается прилегающими к поверхности доменами треугольного сечения, в которых векторы намагниченности параллельны поверхности. Внутри монокристалла блоховские доменные границы, параллельные ОЛН, источников поля не содер- жат. Предполагалось, что внутренние границы доменов треугольного сечения ориентированы так, что нормальные составляющие векторов по обе стороны от границы равны. Из этого предположения следует, что угол на рис. 5 равен 900, а интеграл от источника поля , вычисленный по отрезку нормали от одной стороны границы до другой, равен нулю для любого типа границы, при условии, что векторное поле однородно в любом сечении граничного слоя, параллельном поверхности границы. Равенство нулю интеграла означает, что внутри границы источники магнитостатического поля взаимно компенсируются в пределах ширины границы, и можно считать, что достаточно узкая граница не создаёт источников поля. Исходя из изложенных выше предположений в работе [1] установлено, что минимум свободной энергии достигается при ширине доменов , где A и K - соответственно константы обмена и анизотропии, L - размер монокристалла в направлении ОЛН (рис. 6).

Рис. 6

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

, (19)

где ; - время; - гиромагнитное отношение электрона; - безразмерный параметр, определяющий вклад диссипативного члена. Приводимые ниже результаты получены при . - безразмерный вектор эффективного поля с компонентами:

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

, .

В расчетах использовались следующие значения магнитных параметров: а) для Ni80Fe20 - А=1.310-6 Эрг/см, Ms=800 Гс, K1=5103 Эрг/см3, K2=0 Эрг/см3; б) для Nd2Fe14B - А=1.710-6 Эрг/см, Ms=1275 Гс, K1=4.5107 Эрг/см3, K2=6.6106 Эрг/см3. Параметры Co приведены на стр. 14. Значения для Ni80Fe20, Co и Nd2Fe14B составляют 7.05 нм, 63.3 нм и 301 нм соответственно.

Здесь и далее полученные решения проверялись заданием случайных возмущений в равновесные векторные поля с последующим контролем их возвращения в исходные состояния и проецированием рассчитанных равновесных полей на более мелкие сетки. Расчеты проводились на сетках, содержащих от 8181 до 10411041 точек. Во всех случаях шаг сетки выбирался меньше безразмерной характерной ширины блоховской доменной границы /L. Векторное поле, рассчитанное на мелких сетках, изображалось точками следующим способом: если вектор m отклонялся от направления оси Ox не более 450, то на рисунках он изображался серой точкой. При противоположных направлениях - светлой. В остальных случаях (ориентации преимущественно поперёк ОЛН) - тёмной. На некоторых рисунках показана проекция поля m на крупную сетку. Фрагменты поля m показаны на расчетных сетках.

Рис. 7

Для монокристалла Co из начального однородно намагниченного состояния, когда векторы m направлены по оси Ох (вдоль ОЛН), получены следующие результаты. При размере монокристалла оба метода расчета (минимизация дискретного аналога функционала (1) и решение дискретного аналога уравнения (19)) показали одинаковые результаты: в монокристалле образуется вихревая структура, показанная на рис. 7. В этой конфигурации компоненты намагниченности монокристалла и близки к нулю, а .

Уменьшение размера L до 0.2 приводит к увеличению до 0.96, что означает квазиоднодоменное состояние с намагниченностью вдоль оси Oz.

При решение уравнения Ландау-Лифшица привело к доменной структуре, а минимизация функционала дала качественно другой результат: доменную структуру, сходную с показанной на рис. 5б. Первую будем обозначать II2, вторую II3 (по числу полосовых доменов). Обе структуры соответствуют локальным минимумам свободной энергии , причем , а (компоненты энергии записаны в таблице). Равновесное состояние (среди обнаруженных), в котором энергия системы минимальна, будем условно называть стабильным. В данном случае это II2.

При и оба метода дают одинаковый результат - трехполосную структуру II3 (рис. 5б), в которой относительный размер треугольных доменов уменьшается с ростом L. Для расчеты проводились только с помощью уравнения (19). При получена структура, показанная на рис. 9а, которую обозначим III3. Структуры типа III представляют собой структуры типа II с внедренными в полосы клиновидными доменами. Нижний индекс обозначает количество полосовых доменов. Увеличение размера до и дает аналогичную структуру с клиновидными доменами большего размера при относительном уменьшении треугольных доменов.

При данный подход позволил обнаружить состояния, энергетически более выгодные по сравнению с полученными из однородно намагниченных по ОЛН. Из таблицы следует, что для появление в трехполосной структуре II3 клиновидных доменов значительно снижает магнитостатическую энергию при некотором повышении энергии обмена

Таблица 1. Результаты расчетов для монокристаллов Co. - энергия обмена, энергия в собственном поле и энергия анизотропии

Начальное условие	Результаты расчетов
Размер монокристалла (L)	Поле m в начальном состоянии	Тип доменной структуры	Компоненты энергии	Значение	Тип состояния
1 dc	по оси Оx	II2	0,35	0,39	0,66	1,40	стабильное
1 dc	по оси Оx	II3	0,67	0,17	1,00	1,84	метастабильное
1 dc		II2	0,35	0,39	0,66	1,40	стабильное
3 dc	по оси Оx	II3	0,29	0,20	0,43	0,92	стабильное
3 dc	II2	II2	0,11	0,39	0,48	0,98	метастабильное
3 dc		II4	0,45	0,15	0,48	1,08	метастабильное
3 dc		II5	0,63	0,12	0,57	1,32	метастабильное
3 dc		II5	0,63	0,12	0,57	1,32	метастабильное
10 dc	по оси Оx	II3	0,08	0,19	0,29	0,56	метастабильное
10 dc		II4	0,13	0,13	0,29	0,52	стабильное
10 dc		II5	0,19	0,10	0,26	0,55	метастабильное
10 dc		II6	0,22	0,08	0,28	0,58	метастабильное
15 dc	по оси Оx	III3	0,09	0,13	0,26	0,48	метастабильное
15 dc	II3 (рис. 4б)	II3	0,05	0,19	0,27	0,51	метастабильное
15 dc		II4	0,09	0,12	0,23	0,44	метастабильное
15 dc		II5	0,13	0,09	0,21	0,43	стабильное
15 dc		II6	0,17	0,07	0,22	0,46	метастабильное
20 dc	по оси Оx	III3	0,08	0,12	0,24	0,44	метастабильное
20 dc	по оси Оy	II11	0,27	0,04	0,23	0,54	метастабильное
20 dc	по оси Оz	III3	0,08	0,12	0,24	0,44	метастабильное
20 dc	II4	II4	0,07	0,13	0,21	0,41	метастабильное
20 dc	II5	II5	0,10	0,09	0,19	0,38	метастабильное
20 dc		III5	0,11	0,08	0,19	0,38	стабильное
20 dc		III6	0,13	0,07	0,19	0,39	метастабильное
40 dc	по оси Оx	III5	0,07	0,06	0,15	0,28	-

и понижении энергии анизотропии . Повышение объясняется увеличением количества доменов и, следовательно, площади доменных границ, а понижение - уменьшением размера треугольных замыкающих доменов, в которых векторы m ориентированы поперек ОЛН. клиновидные домены дают такой же эффект (структуры II5 и III5 в таблице), но более слабый из-за малых размеров клиньев в пятиполосной структуре. В монокристалле размером из начальной намагниченности по оси Ох получается пятиполосная структура, похожая на структуру в монокристалле с , но с увеличенными размерами клиньев. По данным работы [7] клиновидные домены в монокристалле Со в стабильной структуре должны появляться при. В приведенных расчётах они появились при в метастабильной структуре (рис. 9а) и только при в стабильной. Однако подстановка показанной на рис. 8а структуры в качестве начального состояния для показывает, что клиновидные домены в метастабильном состоянии существуют и при меньших размерах монокристалла.

Рис. 8

Рис. 9

В высокоанизотропном магнетике Nd2Fe14B из-за узкой доменной границы расчеты проводились только для . Оказалось, что однородно намагниченный вдоль ОЛН монокристалл этих размеров, в отличие от Со, под действием собственного поля не размагничивается. Но при включении достаточно большого внешнего поля Hext в направлении, противоположном m, в углах появляются зародыши перемагничивания, которые, прорастая внутрь монокристалла, приводят к перемагничиванию монокристалла в направлении Hext.

Возникающие в постоянном внешнем поле доменные структуры неравновесны, т.к. величины поля, в котором появляются зародыши, оказывается достаточно для полного перемагничивания монокристалла. Однако при импульсном включении поля на короткое время зародыши перемагничивания могут либо исчезать (при малом размере), либо прорастать, образуя трехполосную равновесную структуру типа II3. Другой способ получения многодоменнной структуры состоит в намагничивании монокристалла в направлениях, отличных от ОЛН, с последующим отключением внешнего поля. При малых размерах (), когда в равновесии возможно только однодоменное состояние, поле намагниченности в Nd2Fe14B, в отличие от Со, ориентируется поперек монокристалла (по ОЛН). Для рассматриваемых монокристаллов можно получить критерий, показывающий, когда в однородно намагниченном состоянии поле намагниченности будет ориентировано вдоль монокристалла, а когда поперек. Если поле m ориентировано по оси Ox, то магнитостатическое поле создают две бесконечные полосы (грани призмы) и , на которых поверхностная плотность источника равняется, соответственно, -1 и 1. Потенциал поля в произвольной точке сечения с координатами ,, можно вычислить как сумму вклада полос:

.

Используя выражение для потенциала, получаем магнитостатическую энергию системы:

.

Безразмерная величина и не зависит ни от размера L, ни от магнитных параметров монокристалла. При продольной ориентации однородного поля m энергия системы , а при поперечной (по ОЛН) . Отсюда получаем критерий: если , то в однодоменном состоянии энергетически выгодна продольная ориентация однородного поля m, если , то поперечная. Для Со и Nd2Fe14B величина составляет 2.58 и 31.7 соответственно. Таким образом, результаты расчета поля m в мелких кристаллах соответствуют полученному критерию.

Даже в сравнительно крупном монокристалле Ni80Fe20 () поле m из состояния однородной намагниченности по оси x из-за низкой анизотропии ориентируется в однородное, направленное вдоль монокристалла. Увеличение размера до (5.16 мкм) приводит из того же начального состояния к вихревому состоянию, близкому к показанному на рис. 7. При размере (11.61 мкм) из начального состояния получается структура Ландау , сходная с показанной на рис. 5б. В работе [8] приводятся экспериментальные и расчетные данные о существовании стабильного состояния доменной структуры Ландау в пластинке нм с ОЛН, ориентированной вдоль длинной стороны. По данным настоящей работы стабильным состоянием бесконечно длинного монокристалла с нм является однородно намагниченное по оси z. Расхождение с результатами работы [7] объясняется влиянием поверхностей пластинки, ортогональных оси z, которых нет у бесконечно длинного монокристалла.

Таким образом, несмотря на некоторые расхождения с доменной структурой на рис. 6, результаты проведенного численного моделирования в монокристалле Со показали возможность существования доменных структур Ландау по крайней мере в диапазоне размеров . При меньших размерах структура Ландау трансформируется в вихревую (рис. 7), а затем в однодоменную. Как видно на рисунках, появление клиновидных доменов приводит к уменьшению относительных размеров треугольных доменов, в которых намагниченность ориентирована поперек ОЛН. Так что и в больших монокристаллах структура Ландау должна вырождаться. Расхождения рассчитанных доменных структур от показанной на рис. 6 состоят также в расположении и форме треугольных доменов, которые не примыкают друг к другу, а углы для Со могут значительно отличаться от 900. Из таблицы видно, что во многих случаях компонента дает существенный вклад в свободную энергию , что противоречит исходным соображениям, на основе которых построена структура на рис. 6.

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

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

3. Микромагнитное моделирование эффекта термического намагничивания

Моделируются механизмы термического намагничивания (ТН) в многослойной системе (многослойной пленке) и в бесконечно длинном монокристалле, методы расчета полей намагниченности в котором представлены в главе 2.

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

Показано, что эффект ТН отсутствует в системе, состоящей из слоев Nd2Fe14B со случайной ориентацией ОЛН в каждом слое и имеющей немагнитные прослойки, исключающие обменное взаимодействие между слоями. Если прослойки исключить, то появляется эффект ТН, зависящий от ориентации легких осей. Наибольшая величина эффекта наблюдается при ориентации ОЛН в плоскости слоев, когда собственное поле отсутствует. При увеличении случайно выбираемых углов отклонения от плоскости слоев собственное поле возрастает, а эффект ТН уменьшается. Таким образом, в многослойной системе эффект ТН появляется при включении обменного взаимодействия между слоями, а собственное поле препятствует термическому намагничиванию.

Численное моделирование показало существование эффекта ТН и бесконечно длинном монокристалле Nd2Fe14B с поперечной ориентацией ОЛН. В некотором диапазоне поперечных размеров между и стабильным является однодоменное состояние, но возможно и метастабильное размагниченное. При нагревании монокристалла в метастабильном состоянии коэффициенты анизотропии уменьшаются, и монокристалл переходит в стабильное (намагниченное) состояние с продольной ориентацией поля m. Намагничивание сопровождается уменьшением энергии в собственном поле и энергии обмена при росте энергии анизотропии. Последующее охлаждение приводит только к повороту поля m к направлению ОЛН с сохранением намагниченного состояния. В монокристалле с возможно только однодоменное состояние, а при поле m после охлаждения возвращается в исходное состояние, и монокристалл размагничивается.

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

Предлагается метод расчета распределения намагниченности в монокристалле. Метод основан на разделении монокристалла на две области. Конечную, прилегающую к торцу монокристалла, с трехмерным полем намагниченности и остальную (полубесконечную), в которой поле намагниченности предполагается двумерным. При расчете намагниченности в конечной области учитывается влияние магнитостатического поля полубесконечной части монокристалла.

Рассматривается бесконечно длинная монокристаллическая призма одноосного ферромагнетика с квадратным поперечным сечением . Расположение координатной системы показано на рис. 14а. В области D, прилегающей к торцу монокристалла, распределение намагниченности предполагается трехмерным и находится путем вычисления стационарного решения уравнения Ландау-Лифшица (19), при условии, что в полубесконечной части монокристалла ниже области D поле m в любом поперечном сечении совпадает с полем m на нижней границе области D. В трехмерном случае эффективное поле H в (19) можно записать в виде:

, (20)

где w - единичный вектор направления ОЛН. На гранях призмы D для уравнения (19) ставилось граничное условие: , где n - единичная внешняя нормаль к поверхности кристалла. На верхней и боковых гранях это условие означает отсутствие поверхностной анизотропии, на нижней - неизменность распределения намагниченности m вдоль монокристалла (в этом случае n - нормаль к нижней грани призмы D).

Рис. 10

При численном решении задачи координатное пространство заполняется равномерной сеткой с шагом h, ячейка которой с центром в узле показана на рис. 10 б. Предполагается, что индексы i, j и k возрастают в направлениях координат x, y и z соответственно. В дальнейшем центр ячейки будем обозначать одной буквой - N. Вычисление потенциала собственного поля строилось на основе формулы (10), записанной для кубической ячейки (рис. 10 б). В качестве расчетной области для потенциала выбиралась призма , верхняя и боковые грани которой располагались в пустом пространстве параллельно соответствующим граням призмы D на расстоянии нескольких шагов h, а нижняя грань содержала нижнюю грань D. Сеточное уравнение для потенциала получено методом, который использовался в главе 2. В данном случае замена интегралов в (10) дискретными аналогами дает сеточное уравнение в безразмерном виде:

. (21)

Для внутренних сеточных узлов расчетной области функцию в правой части уравнения (21) можно записать, используя обозначения , , следующим образом. Для точек N, лежащих внутри призмы D: . Для точек N в пустом пространстве: . Для не лежащих на ребрах (регулярных) точек верхней и боковых граней призмы D: , где p и q индексы координатных осей, ортогональных к вектору нормали . Так, например, если узел N лежит на правой грани, то , . В результате получаем: . Если узел лежит на боковых или верхних ребрах (но не является вершиной трехгранного угла), то , где - внешние нормали граней, образующих ребро, а p - индекс координатной оси, ортогональной векторам. Например, для правого верхнего ребра: . В результате: . Для четырех трехгранных углов, примыкающих к верхней грани: , гдеи - внешние нормали граней, образующих трехгранный угол с вершиной в точке N. Например, для переднего правого угла: и .

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

Формула (21) аппроксимирует уравнение (2) с порядком внутри и снаружи призмы D, а условие (3) с порядком в регулярных точках граней. Сеточную функцию можно рассматривать в качестве дискретного аналога объемной плотности источника потенциала . Причем, если внутри призмы D функция является обычной аппроксимацией объемного источника , то на ее гранях - эффективный источник, учитывающий поверхностный источник с плотностью .

Для уравнения (21) на границе расчетной области ставилось условие Дирихле (задавался потенциал), который вычислялся как сумма:

, (22)

где - вклад области D, а - вклад полубесконечной части стержня, расположенной ниже области D. Поскольку при построении сеточного уравнения (21) поверхностный источник потенциала учитывался через эффективный объемный источник, то вычисление функций и проводилось по формулам ньютоновского потенциала [4] как результат вклада только объемного источника . Для вклада области D формула ньютоновского потенциала имеет вид:

, (23)

где - точка (радиус-вектор) в которой вычисляется потенциал, а переменная интегрирования из области D.

Для вычисления рассмотрим потенциал , создаваемый отрезком стержня с поперечным сечением G, в котором функция q зависит только от координат x и y. В соответствие с формулой ньютоновского потенциала можно записать:

, (24)

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

. (25)

Для рассматриваемой задачи сечением G является безразмерный квадрат , а в качестве координаты выбиралась z координата нижней грани области D. При построении дискретных аналогов интегралов (23) и (25) объемная плотность q заменялась сеточной функцией , причем для интеграла (25) вместо функции использовались значения , лежащие на шаг h выше нижней грани области D. Учитывалось, что в точках , лежащих на нижней грани области D, эти интегралы имеют особенность (несобственные второго рода).

Предлагаемый метод можно применить и для конечных стержней. В этом случае вклад в потенциал отрезка стержня , в котором намагниченность зависит только от координат x и y, должен вычисляться по формуле (24). В случае коротких стержней (пластинок), в которых отсутствует участок с двумерным полем m, рассматривается только конечная область D с трехмерным полем, в которой нижняя грань с прилегающими к ней двугранными и трехгранными углами не отличается от остальных граней призмы.

В качестве примера применения расчетного метода рассматривается монокристалл Со размером . Коэффициент диссипации . Предполагается, что ось легкого намагничивания ориентирована по координатной оси x. В глубине монокристалла (на нижней грани D) распределение намагниченности соответствует распределению, полученному методом минимизации свободной энергии в бесконечно длинном монокристалле такого размера (рис. 5б). На торцевой поверхности (верхней грани D) наблюдается вполне ожидаемое незначительное изменение поля намагниченности из-за влияния поверхностных источников магнитостатического поля. Однако если в состоянии коэффициенты и уменьшить в 10 раз, а затем увеличить до исходного значения, то распределение намагниченности на поверхности станет качественно отличным от распределения в глубине. По виду торцевой поверхности можно сделать вывод, что образец находится в однодоменном состоянии, хотя в глубине монокристалла распределение намагниченности неоднородно. Поскольку коэффициенты анизотропии уменьшаются при нагревании образца, полученную выше трансформацию доменной структуры можно рассматривать как результат цикла нагревание - охлаждение. Приведенные в главе 4 данные других расчетов показывают, что полученный на частном примере результат: в стабильном состоянии поле намагниченности на поверхности образца соответствует полю в глубине, а в метастабильном - не соответствует, не является общим.

5. Микромагнитное моделирование распределения намагниченности в тонких пластинках

Предлагается конечно-разностный метод расчета распределения намагниченности в тонкой монокристаллической пластинке. На основе разработанного метода рассчитываются доменные структуры в пластике Nd2Fe14B, которые появляются при изменении констант анизотропии. Проведено сравнение с экспериментально наблюдаемыми доменными структурами в пленках Nd2Fe14B и . Методами, изложенными в главах 4 и 5, исследуется изменение структуры доменных границ в пластинке Со при увеличении ее толщины и моделируется магнитная запись на доменных границах монокристаллической пленки Nd2Fe14B.

Рис. 11

В ряде случаев в расчетах распределения намагниченности m в достаточно тонких пластинках можно предполагать, что поле m изменяется только в плоскости пластинки, т.е. является двумерным. Поскольку поле намагниченности предполагается одинаковым в любом сечении, параллельном плоскости пластинки, в качестве расчетной области рассматривается среднее сечение D, в центре которого расположено начало координатной системы (рис. 16а). Потенциал U и его градиент , в отличие от намагниченности m, являются функциями трех координат. Рассмотренные в главе 1 методы предполагали, что для вычисления собственного поля требуется либо вычисление U на трехмерной сетке, либо вычисление трех компонент в каждом сеточном узле как трех интегральных сумм - вклада отдельных ячеек. Описанный ниже метод вычисления магнитостатического поля позволяет ограничиться двумя суммами и двумерной сеткой.

Поле намагниченности в области D рассчитывается путем нахождения стационарных решений уравнения Ландау-Лифшица (19), в котором эффективное поле имеет вид (20). Потенциал магнитостатического поля U вычисляется как сумма вкладов объемных источников с плотностью и поверхностных источников с плотностью , где n - внешняя нормаль к поверхности пластинки. Используя теорему Гаусса, получаем следующее выражение для потенциала:

. (26)

Здесь - точка (радиус-вектор), в которой вычисляется потенциал, - переменная интегрирования, которая в тройных интегралах изменяется по объему пластинки V , а в поверхностном - по ее поверхности S. На боковых гранях пластинки для уравнения (19) задается граничное условие , которое означает отсутствие поверхностной энергии. На верхней и нижней гранях пластинки выполнение этого условия следует из постановки задачи.

При численном решении задачи среднее сечение пластинки D покрывалось равномерной сеткой с шагом h, на которой индексы сеточных узлов i и j возрастают в направлениях координат x и y соответственно (рис. 16б). Аппроксимация трехмерного вектора на двумерной сетке выполнялась следующим способом. Во внутренних и регулярных граничных (не угловых) узлах сетки методом, описанным ниже, вычислялся сеточный аналог потенциала . Затем по обычным конечно-разностным формулам во внутренних узлах находились первые две компоненты сеточного аналога вектора : , . Третья компонента , как и сеточная функция , вычислялись следующим способом. Из формулы (26) при следует:

...