朗道能级相关

% **************************

L = 15; % Interval Length

N = 200; % 设置 N*N 格点

% x = linspace(-L/2,L/2-L/(N),N)'; % Coordinate vector， 周期性边界时

x = linspace(-L/2,L/2,N)'; % Coordinate vector，硬边界时

dx = x(2) - x(1); % Coordinate step

% operators in 1D:

e = ones(N,1); 

% Dif1 = spdiags([-e 0*e e -e e],[-1 0 1 N-1 -N+1],N,N)/(2*dx);% 有周期性边界

% Lap1 = spdiags([e -2*e e e e],[-1 0 1 N-1 -N+1],N,N)/dx^2;

Dif1 = spdiags([-e 0*e e],[-1 0 1],N,N)/(2*dx);% 硬边界

Lap1 = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2; 

I1 = speye(N,N);

% operators in 2D:

Lap2 = kron(Lap1,I1)+kron(I1,Lap1); % 2D Laplacian

PDx = kron(Dif1,I1); % partial diff of x

PDy = kron(I1,Dif1);% partial diff of y

xy = spdiags(x,0,N,N);

X = kron(xy,I1); % x position operator

Y = kron(I1,xy); % y position operator

R2 = kron(xy.^2,I1)+kron(I1,xy.^2); % R^2=x^2+y^2

L3 = -1i*(X*PDy-Y*PDx);% angular momentum, z向

B = 4; % 磁感应强度

H = 0.5*(-Lap2+B*L3+0.25*B*B*R2); % 哈密顿算子

NT = 50000; % No. of time steps

TF = 50; T = linspace(0,TF,NT); % Time vector

dT = T(2)-T(1); % Time step

% E = expm(-1i*H*dT); % Evolution

I = speye(size(H));

A = -1i*H*dT;

E = I;

for m=1:11

    E = E+(A^m/factorial(m));

end

% spy(E)

% **************************

e1 = ones(N,1);

xf = kron(x,e1);

yf = kron(e1,x);

NLL = 4;% Landau 能级量子数，与z角动量量子数

a = 3;% 控制波包大小

d = 0.5*L;% 平移距离，平移波函数时使用

% 以原点为中心的旋转对称的朗道能级波函数

psi0 = ((xf+1i*yf).^NLL).*exp(-0.25*B*(xf.^2+yf.^2)); 

% 中心按照 d 在 x 方向平移的朗道能级波函数

% psi0 = exp(-1i*0.5*B*yf*d).*(((xf-d)+1i*yf).^NLL).*exp(-0.25*B*((xf-d).^2+yf.^2));

% psi0 = psi0.*exp(-((xf-0.3*L).^2+(yf).^2)/a); % 用于截取波包

psi0 = psi0/sqrt(psi0'*psi0); %归一化

EM=psi0'*H*psi0 % 能量期望

LM=psi0'*L3*psi0; % Lz 期望

% 理论预测的能量本征值是 B*(2*NLL+1)/2

ES=sqrt(real(psi0'*H*H*psi0-EM.^2)) % 能量标准差

LS=sqrt(real(psi0'*L3*L3*psi0-LM.^2)); % Lz 标准差

% 计算密度

rho = abs(psi0).^2; 

rho1 = reshape(rho,[N,N]); % 变成二维数组，方便画图

% 在密度足够大的范围内计算相位，阈值设为0.05*h

h=1e-3;

 phase = reshape(angle(psi0).*heaviside(abs(psi0).^2-0.05*h),[N,N]);