一维有限深方势阱中的光电离模拟（Matlab）

贡献者： addis

本文处于草稿阶段。

预备知识　Crank-Nicolson 算法解一维含时薛定谔方程（Matlab)

　　本文使用原子单位制。一维含时薛定谔方程为（式 9 ）

\begin{matrix} (1) & - \frac{1}{2} \frac{\partial^{2} ψ}{\partial x^{2}} + V (x, t) ψ = i \frac{\partial ψ}{\partial t} . \end{matrix}

若在一维薛定谔方程中加入电场，则势能可以分为与时间无关的势阱与含时的电场项两部分

\begin{matrix} (2) & V (x, t) = V_{0} (x) - q x E (t) . \end{matrix}

在这里我们把

V_{0} (x)

设为有限深势阱，电场

E (t)

设为高斯波包。

图 1：运行结果。黑色线条为势阱形状，蓝色红色分别为波函数的实部和虚部。

未完成：动画，注意开始时会有噪音波包

未完成：做傅里叶变换验证

ℏ ω - V_{0} = E_{k}

未完成：和一阶微扰的结果做比较

未完成：是否可以做多光子电离？HHG？

　　本代码从 “Crank-Nicolson 算法解一维含时薛定谔方程（Matlab)” 的 TDSE_cn1d.m 修改而来，依赖 “有限深方势阱束缚态程序（Matlab)” 中的 FSW_bound_eng.m 和 FSW_bound_psi.m。

代码 1：TDSE_FSW_PI.m

% 一维有限深方势阱中的光电离模拟
% 等间距网格，稀疏矩阵

function TDSE_FSW_PI
% ==== 参数设置 ======
xmin = -300; xmax = 300; Nx = 6001; % x 网格
tmin = 0; tmax = 225; Nt = 2251; % 时间网格
Nplot = 10; % 画图步数
ax = [xmin, xmax, -0.006, 0.006];
% 方势阱 [-L/2, L/2], 深 -V0
L = 2; V0 = 0.87; m = 1; q = -1;
E0 = 0.01; a = 0.02; w0 = 1; tc = 20; % 电场波包参数
Ef = @(t) E0*exp(-a*(t-tc).^2).*sin(w0*t); % 电场函数
V = @(x,t) V_fun(x,t,Ef,L,V0,q); % 势能函数
n = 1; % 束缚态编号
% ===================
close all;
% 初始波函数
[k, ka] = FSW_bound_eng(L, V0, m);

x = linspace(xmin, xmax, Nx); dx = (xmax-xmin)/(Nx-1);
t = linspace(tmin, tmax, Nt); dt = (tmax-tmin)/(Nt-1);
% 画电场波包
figure; plot(t, Ef(t)); xlabel 'time [au]', ylabel 'E [au]';
psi = FSW_bound_psi(L, k(n), ka(n), mod(n+1,2), x).';
% 准备稀疏矩阵 [对角线, 上对角线, 下对角线]
ind1 = [1:Nx, 1:Nx-1, 2:Nx]; % 行标
ind2 = [1:Nx, 2:Nx, 1:Nx-1]; % 列标
% 动能矩阵非零元
T = -1/(2*m)*[-2*ones(1,Nx), ones(1,2*Nx-2)]/dx^2;
A = (1i*dt/4)*T; % 对角元稍后更新
figure;
set(gcf, 'Unit', 'Normalized', ...
    'Position', [0.1, 0.1, 0.6, 0.3]);
plot(x, real(psi)); hold on;
plot(x, imag(psi));
for it = 1:Nt
    % 更新对角元
    A(1:Nx) = 0.5 + (1i*dt/4)*(T(1:Nx) + V(x, t(it)+dt/2));
    Asp = sparse(ind1, ind2, A, Nx, Nx);
    tmp = Asp \ psi;
    % err = norm(Asp*tmp - psi);
    psi = tmp - psi;
    if mod(it, Nplot) == 0
        clf;
        plot(x, real(psi)); hold on;
        plot(x, imag(psi));
        plot([xmin, -L/2, -L/2, L/2, L/2, xmax],...
            [0, 0, ax(3)/2, ax(3)/2, 0, 0], 'k');
        axis(ax);
        xlabel('x [au]'); ylabel('\psi');
        title(['t = ' num2str(t(it), '%.2f')]);
        pause(0.2);
    end
end
end

% 势能函数：方势垒 + 电场
function ret = V_fun(x,t,Ef,L,V0,q)
ret = zeros(size(x));
ret(x > -L/2 & x < L/2) = -V0;
% 电场势能
ret = ret - (q*Ef(t))*x;
end

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