Радиотехника и электроника, 2023, T. 68, № 1, стр. 75-82

ВВЕДЕНИЕ

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

В нелинейных решеточных системах определяющим фактором возникновения локализованных волн является дискретность. Было показано, что в таких системах существуют пространственно локализованные, периодические по времени колебательные режимы, которые получили название дискретные бризеры (ДБ), или внутренние локализованные нелинейные моды [1, 2]. При изучении ДБ основное внимание сфокусировано на кристаллических решетках, однако ведется поиск ДБ и в спиновых системах [3–9]. В частности, с помощью методов математического моделирования была продемонстрирована возможность их возникновения в гейзенберговских цепочках с одноионной магнитной анизотропией [4]. В настоящее время рассматриваются модели спиновых цепочек с учетом более сложных магнитных взаимодействий, например антисимметричного обмена [10].

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

1. МОДЕЛЬ И ОСНОВНЫЕ УРАВНЕНИЯ

Рассмотрим одномерную цепочку спинов конечной длины L. Пусть $~{{\vec {S}}_{n}}$ – спиновый вектор n‑го узла цепочки. Уравнение движения Гайзенберга для циклической координаты

нормированной на S − длину спинового вектора

имеет вид

где J, ${{H}_{0}}$, D, B – некоторые константы. Уравнение (1) можно рассматривать в качестве модели, являющейся прототипом моноаксиального хирального гелимагнетика [11]. Тогда J – константа обменного взаимодействия, H₀ – величина внешнего магнитного поля, В – приведенная константа одноионной анизотропии, D – константа антисимметричного обмена ДМ. В этом случае должно выполняться условие

параметр $~\left( {{D \mathord{\left/ {\vphantom {D {2J}}} \right. \kern-0em} {2J}}} \right)$ во всех расчетах принимает значение 0.16, что типично для этого класса соединений. Исследуем решения уравнения (1) в виде ДБ или внутренних локализованных мод спиновой цепочки. Случай D = 0, B > 0 как раз соответствует гайзенберговской спиновой цепочке, рассматривавшейся в работе [4]. Ищем решение в виде, предложенном в указанной работе:

причем амплитуда отклонения нормированного спина от оси цепочки ${{s}_{n}}\left( t \right)~$ должна оставаться почти постоянной по всей длине цепочки, за исключением небольших участков. Подставляя в уравнение (1) $s_{n}^{ + }$ в виде (2) и отделяя действительную и мнимую части, получаем систему дискретных уравнений:

Уравнения (1) задают форму волны, уравнения (2) описывают ее движение. Введены следующие обозначения: a – расстояние между соседними узлами, $\Omega = {{\left( {\hbar \omega - {{H}_{0}}} \right)} \mathord{\left/ {\vphantom {{\left( {\hbar \omega - {{H}_{0}}} \right)} {\left( {2JS} \right)}}} \right. \kern-0em} {\left( {2JS} \right)}}$ – эффективная частота, ${\text{tg}}\delta ~\,\, = \,\,\left( {~{D \mathord{\left/ {\vphantom {D {2J}}} \right. \kern-0em} {2J}}~} \right)$, безразмерное время $\tau ~ = {{~t} \mathord{\left/ {\vphantom {{~t} {~{{t}_{0}}}}} \right. \kern-0em} {~{{t}_{0}}}}~$, где ${{t}_{0}} = {\hbar \mathord{\left/ {\vphantom {\hbar {\left( {2JS} \right)}}} \right. \kern-0em} {\left( {2JS} \right)}}$, $\hbar $ – постоянная Планка.

2. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ МАГНИТНЫХ ДИСКРЕТНЫХ БРИЗЕРОВ

Найдем статическое решение системы (3). Из уравнения (3б) в случае статического решения следует соотношение $k~a\,\,~ = ~--{\text{arctg}}\left( {{D \mathord{\left/ {\vphantom {D {2J}}} \right. \kern-0em} {2J}}~} \right)$, означающее, что ДБ будут обладать топологическим зарядом, который придаст им дополнительную устойчивость. Решение первой группы дискретных уравнений (3а) определим с помощью численных методов. Для численного моделирования были использованы цепочки длиной $L = 100$ и 101 узел с открытыми граничными условиями. При открытых граничных условиях уравнения для спинов на концах цепочки будут иметь вид, отличающийся от (3а):

Для заданных значений параметров и приведенной частоты Ω задавалось пробное значение спинового отклонения на крайнем левом узле цепочки s₀, по формуле для левого края цепочки (4а) рассчитывали амплитуду s₁. Далее по общей формуле для внутренних узлов (3а) вычисляли значения s_n до правого края цепочки. Для ускорения расчетов s_n, которые приходилось проводить миллионы раз, была максимально упрощена процедура. Вместо численного решения уравнения (3а) было получено рекуррентное соотношение для расчета s_{n + 1} через s_{n – 1} и s_n. Действительно, относительно s_{n + 1} уравнение (3а) представляет из себя квадратное уравнение:

где

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

Критерием правильности выбранного значения начального отклонения ${{s}_{0}}~$ служила проверка выполнения уравнения для граничных спинов ${{s}_{{L - 1}}}$ и ${{s}_{L}}$ (4б). Если уравнение (4б) выполнялось с заданной точностью, то решение считалось построенным, если нет – бралось новое пробное значение ${{s}_{0}}$.

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

Как известно, состояния солитонного типа возникают при строго определенном наборе параметров В, Ω и стартового значения. Для получения корректных результатов требовалась, как правило, погрешность менее 0.000001, что “де факто” соответствует разнице в 0.0001% между левой и правой частью критерия. Дальнейшее увеличение точности приводило к необходимости рассматривать диапазон с более мелким шагом от стартового значения S, что существенно увеличивало время расчета.

В процессе численного моделирования было установлено, что при фиксированной точности проверки критерия для большей части бризерных решений будет достаточно подобрать стартовое значение s₀ с точностью до 0.000001…0.00000001, но, однако, существуют решения, для которых критерий точности выполняется только после подбора 18-го знака после запятой. Ситуация также усложняется тем фактом, что стартовое значение S₀ может лежать в достаточно большом диапазоне, от 0 до 0.7, который, как было отмечено выше, иногда приходится проверять на корректность с очень маленьким шагом, при этом не следует забывать, что программа проводит моделирование всей длины цепочки (100…101 элемент). Для ускорения вычислений был использован пакет для параллельных вычислений на видеокарте Nvidia CUDA.

В алгоритме [12], разработанном для решения системы уравнений (3а), программа запускалась для фиксированных В, Ω и глобального значения краевого ${{s}_{0}}$, от которого начиналась проверка корректности в сторону увеличения значения. Как правило, значение s₀ выбиралось близким к нулю. Величину глобального шага алгоритма h подбирали под конкретный вид решения, диапазон ее изменений был от ${{10}^{{ - 2}}}$ до ${{10}^{{ - 9}}}$, в зависимости от требований к итоговой величине шага каждой ветви. Шаг алгоритма выглядел следующим образом: программа запускала $2 \times {{10}^{9}}$ параллельных веток, каждая из которых проверяла корректность своего стартового значения, которое вычислялось по формуле ${{s}_{p}} = \,\,~{{s}_{0}}~\,\, + {{hp} \mathord{\left/ {\vphantom {{hp} {(2 \times {{{10}}^{9}})}}} \right. \kern-0em} {(2 \times {{{10}}^{9}})}}$, где p – номер ветви. Для каждой ветви проверялась корректность выполнения граничного уравнения с заданной точностью. Если уравнение выполнялось, программа записывала значение s_p и выводила его на экран. После окончания расчета и проверки всех ветвей к s₀ прибавлялось значение h и алгоритм вновь запускал параллельные ветви уже с новым

С целью апробация правильности работы программы сначала были воспроизведены результаты моделирования бризерных решений для случая D = 0, приведенные в работе [4]. В результате проведенного численного моделирования для случая D, неравного нулю, были выделены группы значений параметров, при которых возможно построение искомых дискретных бризерных решений нескольких типов. При фиксированных значениях В и Ω возможны несколько определенных значений для начального отклонения, которые задают бризеры с различным числом периодов – от одного до многопериодных “бризерных решеток”. Примеры всех полученных нами типов дискретных бризеров приводятся ниже. На графиках по оси абсцисс отложен номер узла, отсчитанный от середины цепочки для более наглядного представления решения, а по оси ординат – значение амплитуды ${{s}_{n}}$.

На рис. 1а–1г представлен пример группы решений, для которых амплитуды спиновых отклонений соседних узлов имеют одинаковые знаки. Решения, представленные на рис. 1а и 1б, получены при одинаковом наборе параметров, но для разного начального отклонения: одиночный кинк и бризерная решетка. Также были найдены решения, для которых не происходит смены знака амплитуды спинового отклонения по всей длине цепочки, как в виде одиночной пары кинк-антикинк (рис. 1в), так и бризерной решетки (рис. 1г). Во всех случаях расчеты проводились на цепочке, содержащей 101 узел.

Рис. 1.

Значение амплитуды s_n в зависимости от номера узла при Ω = –0.28, B = 0.15 (а, б) и Ω = 1.95, B = –1 (в, г): s₀ = = 0.495373257505 (а), 0.4858544687 (б), 0.000308187350000 (в), 0.016876627550000 (г).

Решения, для которых амплитуды спиновых отклонений соседних узлов отличаются знаком, представлены на рис. 2. Среди них также есть и однопериодные (рис. 2а), и бризерные решетки (рис. 2б). Для наглядности расчетные точки на этих графиках соединены линией. Такое “антиферромагнитное” упорядочение спиновых отклонений соседних узлов характерно для больших значений константы анизотропии В. Также используется цепочка, содержащая 101 узел.

Рис. 2.

Значение амплитуды s_n в зависимости от номера узла при Ω = –3.95, B = 4: s₀ = 0.000018542840220 (а), 0.003323453831000 (б).

Численные расчеты проводились на цепочках из 100 и 101 узла, что позволило проследить зависимость вида решения от четности числа узлов в цепочке. Наблюдались как случаи, когда смена четности приводила к качественному изменению вида решения, так и противоположные.

Так, например, для набора параметров Ω = 1.95, B = –1 на цепочке с нечетным числом узлов получаются решения без смены знака амплитуды, приведенные выше (см. рис. 1г), а на цепочке из 100 узлов генерируется решение с периодической сменой знака амплитуды (рис. 3). Оба решения представляют из себя бризерную решетку из двух периодов, но за счет разной четности числа узлов в цепочке первое решение имеет максимумы на узлах (мода Такено–Сиверса [13]), а второе – в междоузлиях (мода Пейджа [14]). В следующем разделе будет показано, что эти решения соответствуют разным траекториям фазового портрета уравнения Дюффинга.

Рис. 3.

Значение амплитуды s_n в зависимости от номера узла при Ω = 1.95, B = –1, s₀ = 0.24708307570; цепочка длиной 100 узлов.

Противоположный пример: при Ω = –3.95, B = 4 на цепочках той и другой четности реализуются однотипные решения (рис. 4а, 4б), также относящиеся к различным модам.

Рис. 4.

Значение амплитуды s_n в зависимости от номера узла при Ω = –3.95, В = 4: а) s₀ = 0.00646293212, длина цепочки 100 узлов, мода Пейджа, б) s₀ = 0.003323453831000, длина цепочки 101 узел, мода Такено–Сиверса.

Рассмотрим кратко, какие изменения по сравнению с результатами работы [4] внесло в дискретные бризерные решения дополнительное слагаемое в уравнении (1). Кроме уже отмеченного в разд. 1 качественного отличия – появления у рассматриваемых решений топологического заряда – имеются и количественные. Приведем сравнение численных решений, полученных при D = 0 и ненулевом D для одного того же набора параметров (рис. 5). Видно, что при сохранении формы решения амплитуда спиновых отклонений уменьшается.

Рис. 5.

Огибающая функция в зависимости от номера узла. Ω = –3.94, В = 4, длина цепочки 100: кривая 1 – D = 0, s₀ = 0.002177; кривая 2 – D = 0.16, s₀ = = 0.00048823.

3. КОНТИНУАЛЬНЫЙ ПРЕДЕЛ

Результаты численного счета показывают, что решения в виде дискретных бризерных мод получаются не для любого набора параметров уравнения (1). Объяснить это можно с помощью приближенного решения системы дискретных уравнений (2а) в континуальном пределе. Для получения уравнения движения спинов в континуальном пределе перейдем от дискретной переменной s_n к непрерывной переменной ψ(z), которую будем считать медленно меняющейся от узла к узлу. Как видим из графиков решений, для ДБ, у которых амплитуда спиновых отклонений соседних узлов не меняет знака (см. рис. 1), роль такой переменной может исполнять сама амплитуда s_n, а для остальных (см. рис. 2) следует взять функцию огибающей $~{{s}_{n}} = {{\left( { - 1} \right)}^{n}}{{\psi }}\left( z \right)$. Выполняя преобразования

в уравнении (1) и пренебрегая нелинейными членами, включающими ${{d{{\psi }}\left( z \right)} \mathord{\left/ {\vphantom {{d{{\psi }}\left( z \right)} {dz}}} \right. \kern-0em} {dz}}$ и $~{{{{d}^{2}}\psi } \mathord{\left/ {\vphantom {{{{d}^{2}}\psi } {d{{z}^{2}}}}} \right. \kern-0em} {d{{z}^{2}}}}$, получаем для непрерывной амплитуды известное уравнение Дюффинга:

Коэффициенты α и β выражаются через параметры цепочки и приведенную частоту [15]. Бризерные моды соответствуют ограниченным периодическим решениям (6). Фазовая диаграмма для уравнения Дюффинга [16] приведена на рис. 6.

Рис. 6.

Фазовый портрет уравнения Дюффинга (6). Траектории а) соответствуют “ярким” решениям (α > 0, β > 0); траектории б) – “темным” (α < 0, β < 0).

Легко видеть, что каждой из найденных путем численного моделирования бризерных мод можно сопоставить одно из возможных ограниченных решений уравнения Дюффинга. Например, решения, приведенные на рис. 1а, 1б, 1в, 1г, отвечают траекториям соответственно на рис. 6б кривая 2, кривая 1, и на рис. 6а кривая 2, кривая 3, решения, приведенные на рис. 2а и 2б – траекториям рис. 6а кривая 2 и кривая 1 соответственно.

Решениями уравнения Дюффинга для различных значений его коэффициентов являются эллиптические функции Якоби (cn, sn, nd). Следовательно, найденные бризерные моды являются кноидальными волнами.

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

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

Рис. 7.

Значение амплитуды s_n в зависимости от номера узла при Ω = –0.28, B = 0.15, s₀ = 0.4858544687, аналитическое решение имеет вид s(z) = sn(0.1609z, 0.89); цепочка длиной 101 узел.

Для континуального описания дискретных бризерных мод мы использовали огибающую функцию ψ. Известно, что поведение огибающей функции для движущихся солитонов в сплошных средах определяется конкуренцией коэффициентов дисперсии и нелинейности. Например, такая ситуация возникает для нелинейного уравнения Клейна–Гордона в модели ${{\varphi }^{4}}$:

которое сводится к нелинейному уравнению Шредингера (НУШ)

с помощью подстановки:

где $~q = {{3\beta } \mathord{\left/ {\vphantom {{3\beta } {2{{\omega }_{0}}}}} \right. \kern-0em} {2{{\omega }_{0}}}}$ и $~\omega _{0}^{2} = \omega _{c}^{2} + {{c}^{2}}{{k}^{2}}$ [16]. В НУШ коэффициенты дисперсии и нелинейности определяются вторым и последним слагаемыми соответственно. Стабильность движущегося волнового пакета в НУШ регулируется критерием Лайтхилла: если $\omega _{0}^{{''}}q > 0$ ($\omega _{0}^{{''}} = {{{{d}^{2}}{{\omega }_{0}}} \mathord{\left/ {\vphantom {{{{d}^{2}}{{\omega }_{0}}} {d{{k}^{2}}}}} \right. \kern-0em} {d{{k}^{2}}}}$), то возникает модуляционная нестабильность, т.е. возможный распад формы волны на последовательность импульсов. Этот случай соответствует ярким солитонам по терминологии, введенной в [17]. В обратном пределе $\omega _{0}^{{''}}q < 0~$модуляционной неустойчивости нет, и этот случай соответствует темным солитонам.

В нашем случае уравнение Дюффинга не содержит “классических” коэффициентов дисперсии и нелинейности. Причина в том, что мы рассматриваем режим не движущихся, а стационарных волн, и рассматриваемые нами ДБ, по сути, представляют собой нелинейные стоячие волны. Как следствие, мы получаем набор дискретных частот вместо закона дисперсии ω = ω(k). Это исключает возможность прямого применения критерия Лайтхилла. Вместо этого для исследования устойчивости дискретных бризерных мод нужно использовать линейный анализ Флоке. Тем не менее, существует аналог критерия Лайтхилла в том смысле, что можно предсказать тип решения, опираясь на знаки коэффициентов α, β в уравнении Дюффинга. Так, если α > 0, β > 0, мы имеем дело с “яркими” модами, при α < 0, β < 0 появляются “темные” моды.

ЗАКЛЮЧЕНИЕ

В работе изложен численный метод решения задачи о нахождении магнитных дискретных бризеров, или внутренних локализованных мод для гайзенберговской спиновой цепочки конечной длины с учетом взаимодействия Дзялошинского–Мория. Для поиска решений использовалось уравнение движения для спинового оператора и открытые граничные условия. Длина цепочки при расчетах имела четное и нечетное число узлов (100 и 101 узел). При проведении численных расчетов использовалась оригинальная программа, максимально упрощающая расчет спиновых отклонений по цепочке и допускающая возможность использования технологий параллельных вычислений. Найдено, что в рассмотренной спиновой цепочке конечной длины существуют несколько видов возможных пространственно-периодических дискретных бризерных решений. Наличие взаимодействия ДМ не препятствует их возникновению, но уменьшает их амплитуду. Результаты численного моделирование возбуждений магнитных дискретных бризеров хорошо согласуются с континуальным приближением, связывающим моды решетки дискретных бризеров и периодические решения нелинейного уравнения Дюффинга. Цепочка магнитных дискретных бризеров является системой, в которой области запасенной энергии связаны с областями аккумуляции поперечных спиновых компонент, т.е. представляет собой одномерный периодический массив резонаторов, запасенная энергия которых может принимать только дискретные значения, и определяется двумя степенями свободы, числом кинков и частотой. Из-за нелинейности количество накопленной энергии в этих “бризерных конденсаторах” будет значительно превосходить емкость линейных стоячих спиновых волн. Эта функциональность может быть использована для создания спинтронных резонаторов на основе киральных гелимагнетиков.

Авторы заявляют об отсутствии конфликта интересов.