自旋波频率与波数关系绘图及MATLAB代码优化问题咨询

一、原代码核心逻辑修正 #

原LLG时间循环中，for i=2:N-1的循环内部错误重复更新了第1个和最后1个自旋，导致边界条件处理混乱。正确的做法是将边界自旋（i=1和i=N）的更新与中间自旋分开，放在外层循环中：

% ---- 修正后的时间演化循环 ----
for t = 1:steps
    S_new = S;
    
    % 更新中间自旋（i=2到N-1）
    for i = 2:N-1
        % 有效场：交换相互作用（左右邻居）
        H_eff = J * (S(i+1,:) + S(i-1,:));
        % LLG方程：进动+阻尼
        precession = -cross(S(i,:), H_eff);
        damping_term = alpha * cross(S(i,:), precession);
        dS = precession + damping_term;
        % 欧拉更新+归一化
        S_new(i,:) = S(i,:) + dt * dS;
        S_new(i,:) = S_new(i,:) / norm(S_new(i,:));
    end
    
    % 更新第1个自旋（周期性边界：邻居是N和2）
    H_eff = J * (S(N,:) + S(2,:));
    precession = -cross(S(1,:), H_eff);
    damping_term = alpha * cross(S(1,:), precession);
    S_new(1,:) = S(1,:) + dt * (precession + damping_term);
    S_new(1,:) = S_new(1,:) / norm(S_new(1,:));
    
    % 更新最后1个自旋（周期性边界：邻居是N-1和1）
    H_eff = J * (S(N-1,:) + S(1,:));
    precession = -cross(S(N,:), H_eff);
    damping_term = alpha * cross(S(N,:), precession);
    S_new(N,:) = S(N,:) + dt * (precession + damping_term);
    S_new(N,:) = S_new(N,:) / norm(S_new(N,:));
    
    S = S_new;
    % ... 其余绘图代码不变
end

二、自旋波频率的获取 #

自旋波频率是横向自旋分量的时间振荡频率，需通过时间序列的FFT提取：

实现步骤 #

1. 存储时间序列：在模拟前初始化数组，记录每个时间步的横向自旋分量（如Sx）：

   % 初始化时间序列存储
   Sx_time = zeros(steps, N);
   Sy_time = zeros(steps, N);
   

   在时间循环中添加记录：

   Sx_time(t,:) = S(:,1);
   Sy_time(t,:) = S(:,2);
   

2. 时间FFT计算频率：模拟结束后，选取单个自旋（或所有自旋的平均）的时间序列，做FFT找到峰值对应的频率：

   function omega = spinwave_frequency_fft(Sx_t, dt)
       Nt = numel(Sx_t);
       s = Sx_t - mean(Sx_t); % 去直流分量
       Y = fftshift(fft(s));
       P = abs(Y).^2;
       
       % 频率网格（Hz，或无量纲频率）
       fgrid = fftshift((0:Nt-1)-floor(Nt/2)) / (Nt*dt);
       mask = fgrid > 0; % 只取正频率
       fpos = fgrid(mask); Ppos = P(mask);
       
       [~, idx] = max(Ppos);
       f_peak = fpos(idx);
       omega = 2*pi*f_peak; % 转换为角频率（rad/时间单位）
   end
   

   调用示例：取中间自旋的Sx时间序列计算频率

   mid_spin = floor(N/2);
   omega = spinwave_frequency_fft(Sx_time(:,mid_spin), dt);
   

三、批量收集ω-k数据并绘图 #

要绘制ω-k关系图，需构造具有确定波数k的初始自旋构型（线性小振幅近似，避免非线性效应），然后遍历不同k值计算对应ω：

1. 封装模拟为函数 #

将整个模拟逻辑封装为函数，输入波数k，输出对应的角频率ω：

function omega = simulate_spinwave_k(N, a, J, alpha, dt, T, k)
    steps = round(T/dt);
    x = (0:N-1)*a;
    % 初始构型：小振幅正弦波（线性区）
    A = 0.05; % 小振幅，保证线性响应
    S = zeros(N,3);
    S(:,1) = A*sin(k*x);
    S(:,3) = sqrt(1 - S(:,1).^2); % 归一化
    
    Sx_time = zeros(steps, N);
    for t = 1:steps
        S_new = S;
        % 中间自旋更新
        for i=2:N-1
            H_eff = J*(S(i+1,:)+S(i-1,:));
            precession = -cross(S(i,:), H_eff);
            damping_term = alpha*cross(S(i,:), precession);
            S_new(i,:) = S(i,:) + dt*(precession + damping_term);
            S_new(i,:) = S_new(i,:)/norm(S_new(i,:));
        end
        % 边界自旋（周期性）
        H_eff = J*(S(N,:)+S(2,:));
        precession = -cross(S(1,:), H_eff);
        damping_term = alpha*cross(S(1,:), precession);
        S_new(1,:) = S(1,:) + dt*(precession + damping_term);
        S_new(1,:) = S_new(1,:)/norm(S_new(1,:));
        
        H_eff = J*(S(N-1,:)+S(1,:));
        precession = -cross(S(N,:), H_eff);
        damping_term = alpha*cross(S(N,:), precession);
        S_new(N,:) = S(N,:) + dt*(precession + damping_term);
        S_new(N,:) = S_new(N,:)/norm(S_new(N,:));
        
        S = S_new;
        Sx_time(t,:) = S(:,1);
    end
    % 计算频率：取所有自旋的平均Sx时间序列
    avg_Sx_t = mean(Sx_time, 2);
    omega = spinwave_frequency_fft(avg_Sx_t, dt);
end

2. 遍历k值并绘图 #

% 参数设置
N = 100; a = 0.5; J = 2.0; alpha = 0.03; dt = 0.01; T = 50;
% 生成k值范围（第一布里渊区：k ∈ [0, π/a]）
k_list = linspace(0.1, pi/a, 20); % 避开k=0（直流分量）
omega_list = zeros(size(k_list));

% 批量计算
for idx = 1:length(k_list)
    omega_list(idx) = simulate_spinwave_k(N, a, J, alpha, dt, T, k_list(idx));
    fprintf('k=%.3f, ω=%.3f\n', k_list(idx), omega_list(idx));
end

% 绘制ω-k关系图
figure('Color','w');
plot(k_list, omega_list, 'bo-', 'LineWidth',1.5, 'MarkerSize',6);
xlabel('波数 k (rad/单位长度)');
ylabel('角频率 ω (rad/时间单位)');
title('一维铁磁自旋波 ω-k 色散关系');
grid on;

四、关键说明 #

• 线性区要求：初始振幅A必须足够小（如0.05），确保自旋波处于线性响应区，此时ω-k关系满足理论色散关系：$\omega = 4J\sin^2(ka/2)$（无量纲单位下）
• 周期性边界：k的取值需满足周期性条件：$k = 2\pi n/(Na)$，n为整数，这样初始构型才符合周期性边界
• 阻尼的影响：小阻尼α（如0.03）会衰减自旋波，但不改变线性区的色散关系，仅降低信号振幅

内容的提问来源于stack exchange，提问作者Tommy Wigenton