二维无质量 Dirac 型粒子的相关代码
二维无质量 Dirac 型粒子相关模拟。
L = 20; % Interval Length
N = 600; % 设置 N*N 格点
x = linspace(-L/2,L/2-L/N,N)'; % Coordinate vector
dx = x(2) - x(1); % Coordinate step
NT = 50000; % No. of time steps
TF = 50; T = linspace(0,TF,NT); % Time vector
dT = T(2)-T(1); % Time step
%**************************
% 设置势能
% 可按需调整
[X,Y]=meshgrid(x,x);
V=150;
lth = 2;
Ub = V*heaviside(X).*heaviside(lth-X);
% Ub = V*exp(-X.^2/lth^2);
UVh = exp(-1i*Ub*dT/2);
% 动量相关
v0 = 10; % 无质量Dirac粒子的速度
e1 = ones(N,1);
K = [2*pi*(-N/2:-1)/L,2*pi*(0:N/2-1)/L]; % 动量范围
Kx = kron(K,e1);
Ky = Kx';
Ka = sqrt(Kx.^2+Ky.^2);
ThK = angle(Kx+1i*Ky);
%**************************
% 初态设置
ao = 1.8;% 控制动量空间波包大小
% 动量中心
ko = 12;
Th = 0;
kxo = ko*cos(Th);
kyo = ko*sin(Th);
% 位矢中心
xo = -2.5;
yo = 0;
% first component in k space
psi01 = exp(-((Kx-kxo).^2+100*(Ky-kyo).^2)/(ao^2))...
.*exp(-1i*(xo*Kx+yo*Ky)).*(cos(ThK)-1i*sin(ThK))/sqrt(2);
% second component
psi02 = exp(-((Kx-kxo).^2+100*(Ky-kyo).^2)/(ao^2))...
.*exp(-1i*(xo*Kx+yo*Ky)).*1/sqrt(2);
% 位矢空间波函数，未归一
psi0 = [reshape( fftshift(ifft2(fftshift(psi01))),N*N,1 )...
;reshape( fftshift(ifft2( fftshift(psi02) )),N*N,1)];
psi0 = psi0/sqrt(psi0'*psi0); %归一化
% 初态的双分量位矢空间波函数
psi21 = reshape(psi0(1:end/2),N,N);
psi22 = reshape(psi0(0.5*end+1:end),[N,N]);
%**************************
% 时间演化和画图
%建立figure和axes
fig=figure;
fig.Position = [0 0 2000 1000];
ax = axes(fig);
h = 20e-5; % z轴相关尺度
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
%时间演化，Strang 方法
psi21 = full(UVh.*psi21);
psi22 = full(UVh.*psi22);
psik1 = fftshift(fft2(psi21));
psik2 = fftshift(fft2(psi22));
psik10 = psik1;
psik1 = cos(Ka*v0*dT).*psik1-1i*exp(-1i*ThK).*sin(Ka*v0*dT).*psik2;
psik2 = cos(Ka*v0*dT).*psik2-1i*exp(1i*ThK).*sin(Ka*v0*dT).*psik10;
psi21 = full(ifft2(ifftshift(psik1)));
psi22 = full(ifft2(ifftshift(psik2)));
psi21 = full(UVh.*psi21);
psi22 = full(UVh.*psi22);
%以上是一步（时间间隔为 dT）的时间演化
% 每隔10步画图
if rem(t,10)
if t~=1
continue
end
end
rho1 = (abs(psi21).^2); % 第一个分量的密度
% rgb颜色，常数项是本底颜色，根号项凸显小振幅，一次和平方项给大振幅产生颜色渐变
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,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 0.5*L]);
ylim(ax,[-0.5*L 0.5*L]);
hold on
rho2 = (abs(psi22).^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,-rho2-0.1*h,1-CO2,'FaceAlpha',0.4);
shading interp
s02.EdgeColor = 'none';
text(ax,-0.4*L,0.4*L,2.1*h,...
['$T=',num2str(dT*t,'%.2f'),'$'],...
'Interpreter','latex','FontSize',20); %时间
text(ax,0.35*L,0.35*L,2.05*h,'方势垒高 ','FontSize',20);
text(ax,0.4*L,0.4*L,1.75*h,...
['$V=',num2str(V,'%.1f'),'$'],...
'Interpreter','latex','FontSize',20);
%势能伪色图
imagesc(ax,[-0.5*L,0.5*L],[-0.5*L,0.5*L],...
h*Ub,'AlphaData',0.1);
%显示相位，阈值是0.0001*h，可按需调整
phase1 = ...
h*angle(psi21).*heaviside(rho1-0.0001*h);
phase2 = ...
h*angle(psi22).*heaviside(rho2-0.0001*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([-15 -32 2.2*h])
camtarget([0 0 0.3*h])
camva(20)
%立刻作图
drawnow
end
