QWZ模型的波包演化模拟
可以用来模拟外加势场的 QWZ 模型。
% 定义晶格和时间
L = 1; % 晶格常数
N = 750; % 设置 N*N 格点
x = linspace(-N*L/2,(N-2)*L/2,N)'; % 1维位矢数组
dx = x(2) - x(1); %
NT = 50000; % 时间步长
TF = 50; %总时长
T = linspace(0,TF,NT); % 时间数组
dT = T(2)-T(1); % 时间步长
%**************************
% % 设置位矢和动量空间的坐标
% 位矢相关
e1 = ones(N,1);
Rx = kron(x',e1); % 2维位矢数组，x坐标
Ry = Rx'; % 2维位矢数组，y坐标
Ra = sqrt(Rx.^2+Ry.^2);% 2维位矢数组，径向坐标
ThR = angle(Rx+1i*Ry);% 2维位矢数组，辐角坐标
% 动量相关
e1 = ones(N,1);
K = linspace(2*pi*(-1/2)/L,2*pi*(N/2-1)/(L*N),N); % 1维动量数组
Kx = kron(K,e1); % 2维动量数组，x
Ky = Kx'; % 2维动量数组，y
% Ka = sqrt(Kx.^2+Ky.^2);
% ThK = angle(Kx+1i*Ky);
% *************************
% 势能设置
G = 1000; % 势能放缩因子
d = 0.46*L*N; %方形势能边界
Hat = @(y) 4*((atan(y/25)+pi/2)/pi).^2; % 平滑的"阶梯"函数
% BOX =@(x,y) G*(Hat(x-d)+ Hat(-x-d))+...
% G*(Hat(y-d)+ Hat(-y-d)); % 边长为 2*d 的方形势阱
% Ub = BOX(Rx,Ry);% 边界势能（方形）
Ub = G*Hat(Ra-d);% 边界势能（圆形）
% thet = 0.3;
% Ub = BOX(Rx*cos(thet)+Ry*sin(thet),...
% -Rx*sin(thet)+Ry*cos(thet));% 边界势能（方形，倾斜）
V11 = exp(-1i*Ub*dT/2); % 势能演化算子，1分量
V22 = V11; %势能演化算子，2分量
% *************************
% 能带参数设置，定义如下：
% H11 = -2*tau*cos(Kx)-2*tau*cos(Ky)-mu
% H22 = 2*tau*cos(Kx)+2*tau*cos(Ky)+mu
% H12 = Delta*(sin(Kx)-1i*sin(Ky))
% H21 = Delta*(sin(Kx)+1i*sin(Ky))
%
Delta = 100;
tau = 100;
mu = 8;
Hkx = Delta*sin(Kx); Hky = Delta*sin(Ky);
Hkz = -2*tau*cos(Kx)-2*tau*cos(Ky)-mu; % 借助Pauli矩阵分解 H
absH = sqrt(Hkx.^2+Hky.^2+Hkz.^2);
%能带（"动能"）给出的时间演化算子
Uk11 = cos(absH*dT)+1i*sin(absH*dT).*Hkz./absH;
Uk22 = cos(absH*dT)-1i*sin(absH*dT).*Hkz./absH;
Uk12 = 1i*sin(absH*dT).*(Hkx-1i*Hky)./absH;
Uk21 = 1i*sin(absH*dT).*(Hkx+1i*Hky)./absH;
% *************************
% 设置 Gauss 波包
ao = 0.01*pi;% 动量分布的标准差
% 中心动量
kxo = 0;
kyo = 2.92;
% 中心位矢
xo = 0;
yo = -100;
% 依据 H(k) 的本征态构造波包
% first component in k space，最后一行由 H(k) 本征矢给出
psi01 = sqrt(1/2)*exp(-((Kx-kxo).^2+(Ky-kyo).^2)/(ao^2))...
.*exp(-1i*(xo*Kx+yo*Ky))...
.*(Hkz-absH); % 为上能带。若为下能带，则取为 (-Hkx+1i*Hky)
% second component，最后一行由 H(k) 本征矢给出
psi02 = sqrt(1/2)*exp(-((Kx-kxo).^2+(Ky-kyo).^2)/(ao^2))...
.*exp(-1i*(xo*Kx+yo*Ky))...
.*(Hkx+1i*Hky); % 为上能带。若为下能带，则取为 (Hkz-absH)
% 位矢空间波函数，未归一
psi0 = [reshape( fftshift(ifft2(fftshift(psi01))),N*N,1 )...
;reshape( fftshift(ifft2( fftshift(psi02) )),N*N,1)];
psi0 = psi0/sqrt(psi0'*psi0); %归一化
%化为2维数组形式
psi21 = reshape(psi0(1:end/2),N,N);
psi22 = reshape(psi0(0.5*end+1:end),[N,N]);
% 绘制能带，及中心动量对应能级的位置
surf(K,K,absH,'EdgeColor','none','FaceAlpha',0.5)
hold on
surf(K,K,-absH,'EdgeColor','none','FaceAlpha',0.5)
hkx = Delta*sin(kxo); hky = Delta*sin(kyo);
hkz = -2*tau*cos(kxo)-2*tau*cos(kyo)-mu;
absh = sqrt(hkx.^2+hky.^2+hkz.^2);% 为上能带。若为下能带，取相反数
scatter3(kxo,kyo,absh,1000,'r.')
hold off
%**************************
% 时间演化以及绘图
fig=figure;
fig.Position = [0 0 2000 1000];
ax = axes(fig);
h = 8e-5;% 绘图的高度范围放缩因子，依据模平方选取
xlim manual
ylim manual
zlim manual
% 自定义颜色函数
CO1=zeros(N,N,3);
CO1(:,:,1) = 0.2*ones(N,N);
CO1(:,:,3) = 0.75*ones(N,N);
CO2 = CO1;
CO2(:,:,3) = 1*ones(N,N);
% 设置显示相位平面的位置
Z0 = -1.5*h*ones(N,N);
for t = 1:NT
%时间演化
%位矢表象，exp(-i U dt/2)
psi21 = V11.*psi21 ;
psi22 = V22.*psi22;
%换为动量表象
psik1 = fftshift(fft2(psi21));
psik2 = fftshift(fft2(psi22));
%exp(-i H(k) dt)
psik10 = psik1;
psik1 = Uk11.*psik1 + Uk12.*psik2;
psik2 = Uk21.*psik10 + Uk22.*psik2 ;
%换为位矢表象
psi21 = full(ifft2(ifftshift(psik1)));
psi22 = full(ifft2(ifftshift(psik2)));
%exp(-i U dt/2)
psi21 = V11.*psi21 ;
psi22 = V22.*psi22 ;
if rem(t-1,20) % 控制图像更新频率
continue
end
rho1 = (abs(psi21).^2); %分量1的模平方
% 画图，着色
CO1(:,:,1) = 0.2*ones(N,N)+10*ones(N,N).*sqrt(rho1/h)./((1+50*rho1/h));
CO1(:,:,2) = 15*ones(N,N).*rho1/h;
CO1(:,:,3) = 0.75*ones(N,N)+ones(N,N).*rho1.^2/h^2;
s01 = surf(ax,x,x,0*(rho1)+0.1*h,CO1,'FaceAlpha',0.4);
shading interp
s01.EdgeColor = 'none';
zlim(ax,[-1.5*h 2.1*h]);
xlim(ax,[-0.5*L*N 0.5*L*N]);
ylim(ax,[-0.5*L*N 0.5*L*N]);
hold on
rho2 = (abs(psi22).^2);% 分量2的模平方
CO2(:,:,1) = 0.2*ones(N,N)+10*ones(N,N).*sqrt(rho2/h)./((1+50*rho2/h));
CO2(:,:,2) = 15*ones(N,N).*rho2/h;
CO2(:,:,3) = 1*ones(N,N)+ones(N,N).*rho2.^2/h^2;
s02 = surf(ax,x,x,0*(-rho2)-0.1*h,1-CO2,'FaceAlpha',0.4);
shading interp
s02.EdgeColor = 'none';
% 文字注释
text(ax,-0.9*L*N,0.4*L*N,2.1*h,...
['$T=',num2str(dT*t,'%.3f'),'$'],...
'Interpreter','latex','FontSize',20);
text(ax,-0.47*L*N,0.48*L*N,1.7*h,...
['$H=\vec{\mathcal{H}}\cdot\vec{\sigma}$,\ ',...
'$\mathcal{H}_x = \Delta \sin(k_x),\ \mathcal{H}_y = \Delta \sin(k_y),\ ',...
'\mathcal{H}_z = -2 \tau(\cos(k_x)+\cos(k_y))-\mu,\ $',newline,...
'$\Delta=',num2str(Delta,'%.0f'),',\ \tau=',num2str(tau,'%.0f'),...
',\ \mu=',num2str(mu,'%.0f'),'$',newline,...
'$k_x(T = 0)\approx',num2str(kxo,'%.1f'),',\ k_y(T = 0)\approx',num2str(kyo,'%.1f'),'$'],...
'Interpreter','latex','FontSize',15);
% 显示势能
imagesc(ax,[-0.5*L*N,0.5*L*N],[-0.5*L*N,0.5*L*N],...
h*Ub*0.01,'AlphaData',0.1);
% 显示相位，仅显示一定阈值（0.01h）以上模平方格点处的相位
phase1 = ...
h*angle(psi21).*heaviside(rho1-0.01*h);
phase2 = ...
h*angle(psi22).*heaviside(rho2-0.01*h);
Ph1 = surf(ax,x,x,-Z0,phase1,'FaceAlpha',0.3);
shading interp
Ph1.EdgeColor = 'none';
Ph2 = surf(ax,x,x,Z0,phase2,'FaceAlpha',0.3);
shading interp
Ph2.EdgeColor = 'none';
hold off
% 其它绘图参数
colormap(hsv)
caxis([-pi*h pi*h])
pbaspect([1 1 0.5]);
campos([850 -850 1.9*h]);
camtarget([0 0 0.24*h])
camva(22)
drawnow
end
sum(sum(rho1))+sum(sum(rho2))
