Визуализация спиновой динамики: моделируем волны намагниченности в ферромагнетиках на Python

4 часа назад Прокомментировать Просмотры: 6

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

Теоретическим фундаментом для описания поведения магнитных моментов служит уравнение Ландау — Лифшица — Гильберта (LLG):

$\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}$

Основные физические величины:

γ — гиромагнитное отношение,
M_s — намагниченность насыщения материала,
α — безразмерный параметр затухания Гильберта (Gilbert damping).

$\mathbf{H}_{\text{eff}} = \mathbf{H}_{\text{ext}} + \frac{2A}{\mu_0 M_s^2} \nabla^2\mathbf{M} + \mathbf{H}_d + \mathbf{H}_{\text{anis}}.$

Эффективное поле формируется суперпозицией отдельных физических вкладов, которые рассчитываются в узлах расчетной сетки:

Внешнее магнитное поле H_ext — задается в качестве граничных или начальных условий эксперимента.
Обменное поле

$\frac{2A}{\mu_0 M_s^2} \nabla^2\mathbf{M}$

— обусловлено квантовомеханическим стремлением соседних спинов выстраиваться сонаправленно. В дискретной схеме оно выражается через оператор Лапласа (вторые пространственные производные по соседним узлам сетки).
Размагничивающее (дипольное) поле H_d — находится из уравнений магнитостатики:

$\nabla\cdot\mathbf{B}=0 , \mathbf{B}=\mu_0(\mathbf{H}+\mathbf{M})$

В вычислительном микромагнетизме его определяют либо прямым интегрированием по объему, либо с помощью алгоритмов быстрого преобразования Фурье (FFT) для периодических структур.

Магнитная анизотропия (кристаллографическая или анизотропия формы) — как правило, моделируется эффективным полем, зависящим от ориентации компонент вектора M.

Первоисточник: [1]

Для начала реализуем базовую одномерную (1D) модель спиновой динамики на Python:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import os
================= ПАРАМЕТРЫ (SE7EN‑режим) =================
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()


В результате получаем наглядную динамическую визуализацию:
Волна намагниченности в 1D
Стоит отметить, что в данном приближении модель упрощена: диссипация и обменное поле временно исключены, за счет чего наблюдается идеальная гармоническая прецессия вокруг вектора внешнего поля.
Теперь перейдем к более сложной постановке — моделированию процессов на двумерной плоскости (2D):
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import os
================= ПАРАМЕТРЫ (SE7EN‑режим) =================
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 = 9         # радиус возмущения 3 ячейки
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 + alpha2))  cross_M_Heff

term2 = -(gamma  alpha / (1 + alpha2)) * 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()


Графический вывод работы скрипта демонстрирует двумерную картину распространения:
Волна на плоскости
На визуализации отчетливо прослеживается периодическое колебательное движение локальных векторов намагниченности.
При значительном увеличении коэффициента затухания (α = 0.5) система переходит в режим сильной диссипации — колебания практически мгновенно гаснут, не успевая завершить даже первый период:
Затухающие колебания
Если поднять величину внешнего магнитного поля до значения 1.0·10⁷, динамика системы трансформируется следующим образом:
Магнитное поле 10е6
При значении H_ext = 1.25·10⁷ фронт возмущения начинает плавно и равномерно распространяться по всей расчетной области:
умеренное внешнее поле
Дальнейшее форсирование внешнего поля (до 4.0·10⁷) приводит к резкому росту групповой скорости волн и возникновению выраженного нелинейного хаотического отклика:
распространение по всей поверхности вещества
Любопытно, что при переходе к значениям порядка 10·10⁷ наблюдается кажущееся замедление фронта при сохранении сложной турбулентной структуры фазового пространства:
Ещё увеличиваем поле
При росте напряженности до 15·10⁷ волновой режим вновь демонстрирует тенденцию к квазистационарной стабилизации:
стабилизирующиеся колебания
Однако экстремальные значения поля снова провоцируют развитие нелинейной неустойчивости:
Огромное поле
При введении высокого коэффициента диссипации (α ~ 50) даже сверхинтенсивные внешние поля (100·10⁷) приводят лишь к крайне медленной диффузии магнитного возмущения:
увеличим сопротивление
Увеличение поля еще на порядок при сохранении высокого затухания формирует регулярную пространственную картину:
Красивое распространение волны
Заключение: В рамках данной работы построена математическая модель и реализован вычислительный алгоритм для анализа спиновой динамики и волновых процессов намагниченности в ферромагнитных средах. Проведен детальный параметрический анализ, продемонстрировавший чувствительность волнового режима к амплитуде внешнего поля и параметрам диссипации. Весь исходный код доступен для воспроизведения и дальнейших экспериментов.
Полученные результаты находятся в качественном согласии с теоретическими концепциями физики магнетизма и выводами профильных исследований по спиновой динамике в низкоразмерных магнитных структурах (Циберкин К. Б. и соавт.). Численный эксперимент наглядно подтверждает применимость конечно-разностной схемы с периодическими граничными условиями для 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




 
Источник