无质量2维Dirac粒子的朗道能级数值模拟
在一定的约定下，可以把哈密顿算子写作
其中 ，矢势 A 使用旋转规范的形式
。
没有出现动量和位矢直接相乘的项，因此可以沿用之前的Time Splitting方法，分别在位矢表象和动量表象处理哈密顿算子中相应的部分。
但是也可以用差分构造离散形式的哈密顿算子，可以帮助计算波函数的能量期望等值。
离散化哈密顿算子：
L = 15; % Interval Length
N = 500; % 设置 N*N 格点. 有大动量时格点不能太稀疏
x = linspace(-L/2,L/2-L/N,N)'; % Coordinate vector
dx = x(2) - x(1); % Coordinate step
% operators in 1D:
e = ones(N,1);
Dif1 = spdiags([-e 0*e e],[-1 0 1],N,N)/(2*dx);% 硬边界
I1 = speye(N,N);
% operators in 2D:
PDx = kron(Dif1,I1); % partial diff of x
PDy = kron(I1,Dif1);% partial diff of y
Px = -1i*PDx;
Py = -1i*PDy;
xy = spdiags(x,0,N,N);
X = kron(xy,I1); % x position operator
Y = kron(I1,xy); % y position operator
B = 4;% 磁场强度
Ax = 0.5*B*(-Y);% vector potential
Ay = 0.5*B*(X);
v0 = 10;
e2 = [1;1];
sigmax = spdiags([e2,e2],[1,-1],2,2);
sigmay = spdiags([-1i*e2,1i*e2],[1,-1],2,2);
I2 = spdiags(e2,0,2,2);
% Hamiltonian
H = v0*(kron(sigmax,Px-Ax)+kron(sigmay,Py-Ay));
构造朗道能级的波函数（当然，也可以是其他的初态函数，可以自己调整），并用离散化的H计算能量等：
e1 = ones(N,1);
xf = kron(x,e1);
yf = kron(e1,x);
NLL = -8; % Landau Level number，整数
psi01 = ((xf-1i*yf).^abs(NLL)).*exp(-0.25*B*(xf.^2+yf.^2));
psi02 = -sign(NLL)*1i*sqrt(2*abs(NLL)/B)*((xf-1i*yf).^(abs(NLL)-1)).*exp(-0.25*B*(xf.^2+yf.^2));
psi0 = [psi01;psi02];
psi0 = psi0/sqrt(psi0'*psi0); %归一化
EM=real(psi0'*H*psi0); % 能量期望
ES=sqrt(real(psi0'*H*H*psi0-EM.^2)); % 能量标准差
EN0 = v0*sqrt(2*B*abs(NLL))*sign(NLL);%理论上应该得到的能量
主要的演化计算还是依靠后面的分割方法：
% 取出波函数的两个分量，分别化为N*N形式
psi21 = reshape(psi0(1:end/2),N,N);
psi22 = reshape(psi0(0.5*end+1:end),[N,N]);
NT = 50000; % No. of time steps
TF = 50; T = linspace(0,TF,NT); % Time vector
dT = T(2)-T(1); % Time step
%磁场相关
B = 4;
e1 = ones(N,1);
Rx = kron(x',e1);
Ry = Rx';
Ra = sqrt(Rx.^2+Ry.^2);
ThR = angle(Rx+1i*Ry);
% 动量相关
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);
% 演化，Strang 方法
for t = 1:NT
% 第一步磁场相关演化
psi210 = psi21;
psi21 = cos(0.25*v0*B*dT*Ra).*psi21-1i*sin(0.25*v0*B*dT*Ra).*1i.*exp(1i*(-ThR)).*psi22;
psi22 = cos(0.25*v0*B*dT*Ra).*psi22-1i*sin(0.25*v0*B*dT*Ra).*(-1i).*exp(1i*(ThR)).*psi210;
% 动能相关演化
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)));
% 第二步磁场相关
psi210 = psi21;
psi21 = cos(0.25*v0*B*dT*Ra).*psi21-1i*sin(0.25*v0*B*dT*Ra).*1i.*exp(1i*(-ThR)).*psi22;
psi22 = cos(0.25*v0*B*dT*Ra).*psi22-1i*sin(0.25*v0*B*dT*Ra).*(-1i).*exp(1i*(ThR)).*psi210;
%密度分布和相位等
rho1 = (abs(psi21).^2);
rho2 = (abs(psi22).^2);
h = 30e-5;
phase1 = ...
h*angle(psi21).*heaviside(rho1-0.0001*h);
phase2 = ...
h*angle(psi22).*heaviside(rho2-0.0001*h);
%按需增加内容
end
