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

贡献者： addis

本文处于草稿阶段。

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

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

\begin{equation} -\frac12 \frac{\partial^{2}{\psi}}{\partial{x}^{2}} + V(x,t)\psi = \mathrm{i} \frac{\partial \psi}{\partial t} ~. \end{equation}

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

\begin{equation} V(x,t) = V_0(x) - qx\mathcal E(t)~. \end{equation}

在这里我们把 $V_0(x)$ 设为有限深势阱，电场 $\mathcal E(t)$ 设为高斯波包。

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

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

未完成：做傅里叶变换验证 $\hbar\omega - 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