图

无限深势阱中的高斯波包

预备知识 无限深势阱, 高斯波包

   我们下面来计算无限深势阱中一个高斯波包的运动. 定性来说, 波包会一边移动一边扩散(变宽变矮), 且在两个势阱壁之间来回反弹. 反弹的过程中会发生干涉.

   结果如图 1 图 2 , 动画见百科互动演示

图
图1:束缚态概率分布, $x$ 轴为束缚态的 $n$, $y$ 轴为概率, 求和约等于 1
图
图2:波包遇到势阱壁后发生反弹, 过程中发生干涉

算法

   本词条使用原子单位, 并假设粒子质量为 1. 假设无限深势阱的区间为 $[0, L]$, 能量的本征波函数(本征态)为(式 4

\begin{equation} \psi _n(x) = \sqrt{\frac{2}{L}} \sinRound{k_n x} \end{equation}
\begin{equation} k_n = \frac{n\pi }{L} \end{equation}
能量的本征值为
\begin{equation} E_n = \frac{k_n^2}{2} \end{equation}

   初始时波函数为高斯波包(式 1

\begin{equation} \psi (x,0) = \frac{1}{(2\pi \sigma _x ^2)^{1/4}} \E^{-(x - x_0)^2/(2\sigma _x)^2} \E^{\I k_0x} \end{equation}

   第一步是把初始波函数投影到能量本征态上

\begin{equation} C_n = \int_0^L \psi _n^*(x) \psi (x,0) \dd{x} = \int_0^L \sqrt{\frac{2}{L}} \sinRound{k_n x} \psi (x,0) \dd{x} \end{equation}
那么接下来, 波函数的演化可以表示为
\begin{equation} \psi (x,t) = \sum _i C_i \E^{-\I E_n t} \psi _n(x) \end{equation}

   若我们假设初始波包宽度比势阱小, 使得波函数在势阱外的波函数可以忽略不记, 则式 5 的定积分可以拓展到无穷区间, 即傅里叶变换. 我们已经知道初始高斯波包(指数)傅里叶变换的结果为式 2

\begin{equation} g(k) = \frac{1}{(2\pi \sigma _p^2)^{1/4}} \E^{-(k - k_0)^2/(2\sigma _k)^2} \E^{-\I x_0(k - k_0)} \end{equation}
将正弦函数记为指数的形式为
\begin{equation} \sqrt{\frac{2}{L}} \sinRound{k_n x} = \I \sqrt{\frac{\pi }{L}} \qtyRound{\frac{\E^{-\I k_n x}}{\sqrt{2\pi }} - \frac{\E^{\I k_n x}}{\sqrt{2\pi }}} \end{equation}
所以式 5 变为
\begin{equation} C_n = \I \sqrt{\frac{\pi }{L}} [g(k_n) - g(-k_n)] \end{equation}
带入式 6 即可.

Matlab 代码

  WpkISW.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
% 无限深势阱中的波包

close all; clear; % 关闭所有画图, 清空所有变量

% === 参数 ===
L = 100; % 势阱区间 [0, L]
sig_x = 4; % 波包宽度
x0 = 50; % 波包初始位置
k0 = 2; % 波包初始动量
Nbasis = 100; % 基底数量
Nx = 501; % 画图格点数
tmax = 100; % 演化时间 [0, tmax]
Nt = 201; % 演化步数
% ============

sig_k = 1/(2*sig_x); % 动量谱宽度
Ax = 1/(2*pi*sig_x^2)^0.25; % x 归一化系数
Ak = 1/(2*pi*sig_k^2)^0.25; % k 归一化系数
% x 和 k 的波函数
f = @(x) Ax * exp(-((x - x0)/(2*sig_x)).^2) .* exp(1i*k0*x);
g = @(k) Ak * exp(-((k - k0)/(2*sig_k)).^2) .* exp(-1i*x0*(k-k0));
% 位置,时间动量格点
x = linspace(0, L, Nx);
t = linspace(0, tmax, Nt);
% 初始波包投影系数
k = (1:Nbasis) * pi / L;
coeff = 1i*sqrt(pi/L)*(g(k) - g(-k));
% 画图判断系数是否收敛
coeff2 = abs(coeff).^2;
figure; plot(coeff2);
title(['sum(|coeff|^2) = ' num2str(sum(coeff2))]);

figure;
for it = 1:Nt
    % 计算 t(it) 时刻的波函数
    psi = zeros(1, Nx);
    for i = 1:Nbasis
        Eng = 0.5*k(i)^2;
        psi = psi + coeff(i) * exp(-1i*Eng*t(it)) * sqrt(2/L) * sin(i*pi*x/L);
    end
    % 画出波函数实部和虚部
    clf; plot(x, real(psi));
    hold on; plot(x, imag(psi));
    xlabel x; ylabel \Psi(x);
    title(['t = ' num2str(t(it), '%.1f')]);
    axis([0, L, -0.5, 0.5]);
    drawnow;
    % 取消注释可将每一帧保存为 png 图片(当前目录下)
    saveas(gcf, [num2str(it) '.png']);
end
致读者: 小时物理百科一直以来坚持所有内容免费且不做广告,这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化,例如大量广告,内容付费,会员制,甚至被收购。因此,我们鼓起勇气在此请求广大读者热心捐款,使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元,我们几天内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时(项目创始人)

编辑词条 返回目录 返回主页 捐助项目 © 小时物理百科 保留一切权利