简谐振子中的高斯波包模拟（Matlab）

贡献者： addis

预备知识　量子简谐振子（级数法），薛定谔方程 2（单粒子多维）

图 1：运行结果，动画见这里。

　　本文使用原子单位制。令无限深势阱的束缚态为 $ψ_{n} (x)$ ，含时波函数为 $ψ (x, t)$ ，已知 $ψ (x, 0)$ ，那么由式 9 得

\begin{matrix} (1) & ψ (x, t) = \sum_{n = 1}^{\infty} C_{n} ψ_{n} (x) e^{- i E_{n} t} . \end{matrix}

其中

\begin{matrix} (2) & C_{n} = \int_{- \infty}^{+ \infty} ψ_{n} (x) ψ (x, 0) d x (n = 0, 1, \dots) . \end{matrix}

　　注意经过一个周期后，波包模方的形状并不会改变，只是相位改变了。这和经典力学是相符的，根据经典的简谐振子模型，无论初始时刻质点的位置和速度是多少，一个周期后都会回到原位。

　　一般来说，波包扩散是因为初始时存在一定的速度分布。例如中心动量为零的波包经过较长时间 $t$ 后，位置分布就会趋近于初始动量分布乘以 $t$ 。

代码 1：SHOqm.m

% 高斯波包在简谐振子势阱中运动
function SHOqm
% === 参数 ===
m = 1; w = 1; % 质量，角频率
xmin = -13; xmax = 13; Nx = 500; % x 网格
tmin = 0; tmax = 2*pi/w; Nt = 200; % 时间网格
Nn = 70; % 束缚态的数量
x0 = 0; sig_x = 2; % 初始位置， 波包宽度
p0 = 5; % 初始动量
t0 = 0; % 高斯波包的初始时间
ax = [xmin, xmax, -1.2, 1.2]; % 坐标范围
% ===========
x = linspace(xmin, xmax, Nx);
psi_bound = zeros(Nn, Nx);

%% 计算并画前几个束缚态
figure;
for n = 0:Nn-1
    psi_bound(n+1,:) = psi_SHO(m, w, n, x);
    if n <= 4
        plot(x, psi_bound(n+1,:)); hold on; grid on;
    end
end
axis([-5, 5, -0.7, 0.8]);
xlabel x; title '\psi_n';
legend({'n=0','n=1','n=2','n=3','n=4'});

%% 计算高斯波包的初始展开系数
C = zeros(1, Nn);
psi_gs = @(x) 1/(2*pi*sig_x^2)^0.25/sqrt(1 + 1i*t0/(2*m*sig_x^2))...
      *exp(-(x-x0-p0*t0/m).^2/(2*sig_x)^2/(1 + 1i*t0/(2*m*sig_x^2)))...
      .*exp(1i*p0*(x-p0*t0/(2*m)));
plot(x, real(psi_gs(x)));

for n = 0:Nn-1
    psi_sho = @(x) psi_SHO(m, w, n, x);
    C(n+1) = integral(@(x) psi_sho(x).*psi_gs(x), xmin, xmax);
end
figure; plot(abs(C).^2);

%% 含时动画
figure; set(gcf, 'Unit', 'Normalized', 'Position', [0.1, 0.1, 0.4, 0.3]);
t = linspace(tmin, tmax, Nt);
for it = 1:Nt
    psi = zeros(1, Nx);
    for n = 0:Nn-1
        En = (0.5 + n)*w;
        psi = psi + C(n+1)*exp(-1i*En*t(it))*psi_bound(n+1,:);
    end
    clf;
    plot(x, real(psi)); hold on;
    plot(x, imag(psi)); axis(ax);
    xlabel x; title(['t = ' num2str(t(it), '%.2f')]);
    set(gca, 'FontSize', 12);
    saveas(gcf, [num2str(it), '.png']);
end
end

% 简谐振子束缚态波函数
function psi = psi_SHO(m, w, n, x)
alpha = sqrt(m*w); u = alpha * x;
psi = 1/sqrt(2^n*factorial(n)) *...
    (alpha^2/pi)^0.25 * hermiteH(n, u) .* exp(-u.^2/2);
end

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。