無質(zhì)量2維Dirac粒子的朗道能級數(shù)值模擬
在一定的約定下，可以把哈密頓算子寫作
其中?，矢勢 A 使用旋轉(zhuǎn)規(guī)范的形式
。
沒有出現(xiàn)動量和位矢直接相乘的項，因此可以沿用之前的Time Splitting方法，分別在位矢表象和動量表象處理哈密頓算子中相應(yīng)的部分。
但是也可以用差分構(gòu)造離散形式的哈密頓算子，可以幫助計算波函數(shù)的能量期望等值。
離散化哈密頓算子：
L = 15; % Interval Length
N = 500; % 設(shè)置 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));
構(gòu)造朗道能級的波函數(shù)（當然，也可以是其他的初態(tài)函數(shù)，可以自己調(diào)整），并用離散化的H計算能量等：
e1 = ones(N,1);?
xf = kron(x,e1);
yf = kron(e1,x);
NLL = -8; % Landau Level number，整數(shù)
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);%理論上應(yīng)該得到的能量
主要的演化計算還是依靠后面的分割方法：
% 取出波函數(shù)的兩個分量，分別化為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
%磁場相關(guān)
B = 4;
e1 = ones(N,1);?
Rx = kron(x',e1);
Ry = Rx';
Ra = sqrt(Rx.^2+Ry.^2);
ThR = angle(Rx+1i*Ry);
% 動量相關(guān)
v0 = 10; % 無質(zhì)量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
????% 第一步磁場相關(guān)演化
? ? 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;
????% 動能相關(guān)演化
? ? 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)));
? ??% 第二步磁場相關(guān)
? ? 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);
????%按需增加內(nèi)容
end
