В этой статье речь пойдёт о том, как смоделировать спиновую динамику и увидеть волны намагниченности прямо на экране вашего компьютера!
Математической моделью указанного выше процесса будет выступать уравнение Ландау-Лифшица-Гильберта:
![\begin{equation} \frac{d\mathbf{M}}{dt} = -\frac{\gamma}{1+\alpha^2} \left[ \mathbf{M} \times \mathbf{H}_{\text{eff}} \right] - \frac{\gamma\alpha}{(1+\alpha^2)M_s} \left[ \mathbf{M} \times \frac{d\mathbf{M}}{dt} \right]. \end{equation}](https://habrastorage.org/getpro/habr/formulas/f/fe/fed/fed47e04c7e15d930e9b413af6d797c8.svg)
Где
γ — гиромагнитное отношение,
Ms — намагниченность насыщения,
α— параметр затухания (Gilbert damping),
эффективное поле собирают из вкладов, которые ты можешь явно посчитать на сетке:
Внешнее поле Hext — задаёшь сам.
Обменное поле
— отвечает за то, что соседние спины хотят быть параллельны. В коде это вторые производные по пространству (лапласиан), т. е. соседи на сетке.
Поле размагничивания (demag field) Hd — это решение магнитостатической задачи
В микромагнетизме его считают либо через интеграл по объёму, либо через быстрое преобразование Фурье (FFT) в периодических задачах.
Анизотропия (кристаллографическая, форма образца) — обычно записывается как поле, пропорциональное компонентам M.
Источник: [1]
Для моделирования микроволн в веществе в 1D напишем код на Python:
import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation import os # ================= ПАРАМЕТРЫ (Хабр‑режим) ================= N = 64 # число ячеек dx = 5e-9 # шаг сетки dt = 1e-12 # шаг по времени (1 пс) steps = 800 # сколько шагов считать save_every = 10 # сохранять каждый N‑й кадр для анимации gamma = 2.211e5 # гиромагнитное отношение alpha = 0.0 # затухание выключено (чистая прецессия) A_ex = 0.0 # обмен выключен (для начала) Ms = 8e5 H_ext = 4.0e4 # УВЕЛИЧЕННОЕ поле (4 Т) — чтобы вращение было видно # Инициализация: M = (1, 0, 0) везде M = np.ones((N, 3)) * np.array([1.0, 0.0, 0.0]) center = N // 2 angle = 0.52 # ~30 градусов M[center, 0] = np.cos(angle) M[center, 1] = np.sin(angle) norm = np.sqrt(np.sum(M**2, axis=1, keepdims=True)) M = M / norm # ================= ФУНКЦИИ ================= def compute_heff(M, H_ext): Heff = np.zeros_like(M) Heff[:, 0] = H_ext return Heff def llg_step(M, Heff, gamma, alpha): cross = np.cross(M, Heff) dM_dt = -gamma * cross # без затухания return dM_dt # ================= СБОР КАДРОВ ДЛЯ АНИМАЦИИ ================= frames = [] print("Считаем кадры для анимации...") for step in range(steps): Heff = compute_heff(M, H_ext) dM_dt = llg_step(M, Heff, gamma, alpha) M_new = M + dM_dt * dt # Нормировка norm = np.sqrt(np.sum(M_new**2, axis=1, keepdims=True)) norm[norm == 0] = 1.0 M = M_new / norm if step % save_every == 0: frames.append(M.copy()) print(f"Всего кадров: {len(frames)}") # ================= ВИЗУАЛИЗАЦИЯ: АНИМАЦИЯ ================= fig, ax = plt.subplots(figsize=(8, 4)) x = np.arange(N) # Рисуем стрелки (векторы) quiver = ax.quiver(x, np.zeros(N), M[:, 0], M[:, 1], angles='xy', scale_units='xy', scale=1.0, color='black', width=0.004) ax.set_xlim(-1, N) ax.set_ylim(-1.2, 1.2) ax.axvline(center, color='red', linestyle='--', linewidth=1, label='центр') ax.legend(loc='upper right') ax.set_title("Спиновая прецессия (анимация): вращение вектора M вокруг H") ax.set_xlabel("ячейка (i)") ax.set_ylabel("компоненты Mx, My") ax.grid(True, linestyle=':', alpha=0.3) time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12, bbox=dict(facecolor='white', edgecolor='none', alpha=0.7)) def init(): quiver.set_UVC(M[:, 0], M[:, 1]) time_text.set_text('') return quiver, time_text def update(frame_idx): M_frame = frames[frame_idx] t_ns = frame_idx * save_every * dt * 1e9 quiver.set_UVC(M_frame[:, 0], M_frame[:, 1]) time_text.set_text(f"t = {t_ns:.2f} нс") return quiver, time_text ani = FuncAnimation(fig, update, frames=len(frames), init_func=init, blit=True, interval=40) # Сохранение в GIF (требуется imagemagick или pillow) output_path = "spin_precession.gif" try: ani.save(output_path, writer="pillow", fps=15) print(f"Анимация сохранена в {output_path}") except Exception as e: print(f"Не удалось сохранить GIF: {e}") plt.show()
Результат- красивая анимация:
Отметим, что здесь модель упрощена- отсуствует затухание и обмен, поэтому происходит чистая прецессия по гармоническом закону.
Теперь напишем код для моделирования в 2D (на плоскости):
import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation import os # ================= ПАРАМЕТРЫ (Хабр‑режим) ================= Nx, Ny = 64, 64 # сетка 64x64 dx = 10e-9 # шаг сетки 10 нм (крупнее → устойчивее и быстрее) dt = 5e-13 # шаг по времени 0.5 пс steps = 1200 # сколько шагов считать save_every = 15 # сохранять каждый 15‑й кадр gamma = 2.211e5 # гиромагнитное отношение alpha = 0.005 # небольшое затухание (чтобы волна не «звенела» вечно) A_ex = 1.3e-11 # обменная константа Ms = 8e5 mu0 = 4 * np.pi * 1e-7 H_ext = 4.0e4 # УВЕЛИЧЕННОЕ поле (4 Т) — чтобы динамика была видна # Инициализация: M = (1, 0, 0) везде M = np.ones((Nx, Ny, 3)) * np.array([1.0, 0.0, 0.0]) # Возмущение в центре: отклоняем на 30 градусов в сторону Y cx, cy = Nx // 2, Ny // 2 angle = 0.52 radius_sq = 100 # радиус возмущения 100 ячеек for ix in range(Nx): for iy in range(Ny): dist_sq = (ix - cx)**2 + (iy - cy)**2 if dist_sq < radius_sq: M[ix, iy, 0] = np.cos(angle) M[ix, iy, 1] = np.sin(angle) # Нормировка norm = np.sqrt(np.sum(M**2, axis=2, keepdims=True)) norm[norm == 0] = 1.0 M = M / norm # ================= ВСПОМОГАТЕЛЬНЫЕ ФУНКЦИИ ================= def compute_exchange_field_2d(M, dx, A_ex, Ms, mu0): coeff = 2 * A_ex / (mu0 * Ms**2) H_ex = np.zeros_like(M) # Лапласиан для каждой компоненты (периодические границы через roll) for i in range(3): d2x = (np.roll(M[:, :, i], -1, axis=0) - 2*M[:, :, i] + np.roll(M[:, :, i], 1, axis=0)) / dx**2 d2y = (np.roll(M[:, :, i], -1, axis=1) - 2*M[:, :, i] + np.roll(M[:, :, i], 1, axis=1)) / dx**2 H_ex[:, :, i] = coeff * (d2x + d2y) return H_ex def llg_rhs(M, Heff, gamma, alpha): cross_M_Heff = np.cross(M, Heff) cross_M_cross = np.cross(M, cross_M_Heff) term1 = -(gamma / (1 + alpha**2)) * cross_M_Heff term2 = -(gamma * alpha / (1 + alpha**2)) * cross_M_cross return term1 + term2 # ================= СБОР КАДРОВ ================= frames = [] print("Считаем кадры для анимации...") for step in range(steps): H_ex = compute_exchange_field_2d(M, dx, A_ex, Ms, mu0) Heff = H_ex.copy() Heff[:, :, 0] += H_ext # внешнее поле только по X dM_dt = llg_rhs(M, Heff, gamma, alpha) M_new = M + dM_dt * dt # Нормировка (обязательно!) norm = np.sqrt(np.sum(M_new**2, axis=2, keepdims=True)) norm[norm == 0] = 1.0 M = M_new / norm if step % save_every == 0: frames.append(M.copy()) print(f"Всего кадров: {len(frames)}") # ================= ВИЗУАЛИЗАЦИЯ: АНИМАЦИЯ ================= fig, ax = plt.subplots(figsize=(8, 8)) # Цвет — это Mx (компонента вдоль поля). Так лучше видно волну. im = ax.imshow(M[:, :, 0], cmap='coolwarm', origin='lower', vmin=-1, vmax=1) cbar = fig.colorbar(im, ax=ax) cbar.set_label('Mx') # Стрелки (не все, чтобы не было каши) stride = 4 x_vec, y_vec = np.meshgrid( np.arange(0, Nx, stride), np.arange(0, Ny, stride), indexing='ij' ) quiver = ax.quiver(x_vec, y_vec, np.zeros_like(x_vec), np.zeros_like(y_vec), color='black', width=0.003, scale=30) ax.set_title("2D спиновая волна: распространение возмущения") ax.set_xlabel("x (ячейка)") ax.set_ylabel("y (ячейка)") time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12, bbox=dict(facecolor='white', edgecolor='none', alpha=0.7)) def init(): im.set_data(M[:, :, 0]) quiver.set_UVC(np.zeros_like(x_vec), np.zeros_like(y_vec)) time_text.set_text('') return im, quiver, time_text def update(frame_idx): M_frame = frames[frame_idx] t_ns = frame_idx * save_every * dt * 1e9 # Обновляем цвет (Mx) im.set_data(M_frame[:, :, 0]) # Обновляем стрелки (Mx, My) u = M_frame[::stride, ::stride, 0] v = M_frame[::stride, ::stride, 1] quiver.set_UVC(u, v) time_text.set_text(f"t = {t_ns:.2f} нс") return im, quiver, time_text ani = FuncAnimation(fig, update, frames=len(frames), init_func=init, blit=True, interval=40) output_path = "spin_wave_2d.gif" try: ani.save(output_path, writer="pillow", fps=15) print(f"Анимация сохранена в {output_path}") except Exception as e: print(f"Не удалось сохранить GIF: {e}") plt.show()
Заданные в коде параметры очень близко описывают пермаллой (Permalloy) — сплав на основе железа и никеля (обычно ~80 % Ni, 20 % Fe).
Ms=8⋅105 А/м — попадает в диапазон 7.5–8.5⋅105 А/мЭто классика для Ni80Fe20.
Aex=1.3⋅10−11 Дж/мA— тоже типично для пермаллоя (часто берут 1.3 или 1.4⋅10−11)
γ=2.211⋅105 рад/(с ⋅А/м)— соответствует гиромагнитному отношению для намагниченности насыщения пермаллоя.
α=0.005 — реалистичное значение для качественного пермаллоя в тонких плёнках. В образцах с хорошей кристаллической структурой и малой дефектностью затухание бывает порядка 0.003–0.008
Вывод программы следующий:
Видно как стрелки совершают колебательное движение.
При увеличении сопротивления alpha до 0.5 колебания полностью затухают, не успев даже сделать один оборот:
Увеличив внешнее поле до 1.0е7, получим следующую анимацию:
Но уже при Hext=1.25е7возмущения медленно распространяются по плоскости:
Если же сделать внешнее поле мощным(до 4.0e7), то колебания будут быстро распространяться дальше по всему веществу и становятся хаотичными:
Удивительно следующее: при большем увеличении поля до 10е7 колебания замедляются, но при этом хаос ярко выражен:
При большем увеличении до 15е7 колебания стабилизируются:
Но при ещё большем увеличении опять становятся хаотическими:
Это можно объяснить следующим образом:
Вначале внешнее поле невелико, доминирующим вкладом в H_eff становится само внешнее поле. В этом случае можно приближённо считать, что H_eff ≈ μ₀H_ext (с учётом магнитной постоянной). Тогда прецессия вектора M происходит с частотой, близкой к ларморовской, и частота колебаний линейно от него зависит.
При увеличении поля вклад других составляющих H_eff (анизотропия, обмен) становится сопоставимым с вкладом внешнего поля. Система становится чувствительной к начальным условиям и малым возмущениям. Решение уравнения ЛЛГ в этой области перестаёт быть простым гармоническим колебанием — появляются сложные, нелинейные траектории.
При достаточно большом внешнем поле система может выйти на новый устойчивый режим.
Однако если увеличить сопротивление хотя бы до 50, то даже при очень сильном поле 100е7 возмущения распространяются очень медленно:
Увеличим поле в 10 раз, оставив сопротивление неизменным:
Результат - красивая анимация
Устойчивость динамики намагниченности оценивалась по локальному показателю Ляпунова, рассчитанному численно путём сравнения двух близких траекторий в фазовом пространстве. При малых значениях внешнего поля показатель λ<0, что соответствует устойчивым регулярным колебаниям (режим ускорения). В промежуточной области наблюдается переход к λ>0, указывающий на возникновение хаотической динамики и чувствительность к начальным условиям. При высоких полях показатель снова становится отрицательным (λ<0), что свидетельствует о восстановлении устойчивости и переходе к режиму вынужденной синхронизации.
При этом показатель Ляпунова на всём диапазоне колебался около нуля, но на подавляющем большинстве интервалов оставался отрицательным.
Вывод: Таким образом, в данном материале сделана математическая и компьютерная модель спиновой динамики и волн намагниченности в ферромагнетиках. Показано, как динамика распространения волн зависит от начальных параметров (внешнего поля и сопротивления).Проведён тщательный анализ полученных результатов, получены очень красивые анимации. Код оставлен читателям в статье для экспериментов.
Наблюдаемая картина согласуется с теоретическими представлениями о спиновых волнах в ферромагнетиках. Полученные результаты качественно соответствуют выводам работ по спиновой динамике в низкоразмерных системах (Циберкин К. Б. и др.). Работа демонстрирует эффективность конечно‑разностной схемы с периодическими граничными условиями для моделирования спиновых волн в 2D.
Хотя спиновая динамика традиционно считается сложной темой, требующей глубокого знания теоретической физики, её основные эффекты могут быть наглядно продемонстрированы с помощью простых численных моделей.
Статья носит чисто демонстрационный и научно-популярный характер.
Автор не имеет специального физико‑математического образования и не является профессиональным физиком-теоретиком!
Если такое дадут на ЕГЭ это будет крах всей системы, я реально занервничал!
Литература и научные статьи по теме:
1. Уравнение Ландау — Лифшица (магнетизм) — Википедия 2. Циберкин К. Б. Диссертационное исследование — ПГНИУ 3. Google Scholar: Циберкин К. Б. (публикация 1) 4. Google Scholar: Циберкин К. Б. (публикация 2) 5. Google Scholar: Циберкин К. Б. (публикация 3) 6. Google Scholar: Циберкин К. Б. (публикация 4) 7. А. И. Ахиезер, В. Г. Барьяхтар, М. И. Каганов. «Спиновые волны в ферромагнетиках и антиферромагнетиках. I» (УФН, 1960, т. 71, с. 533–579) 8. С. Г. Гестрин, Е. А. Сальникова. «Математическое моделирование взаимодействия спиновых волн с дислокациями в ферромагнетиках» 9. «Теория спиновых волн в плёночных ферромагнитных многослойных структурах» (Изд. СПбГЭТУ «ЛЭТИ», 2008) 10. Особенности ферромагнитного и спин-волнового резонансов в магнитных структурах — disserCat