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

贡献者： addis

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

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