可以用來模擬外加勢場的 QWZ 模型。

% 定義晶格和時(shí)間

L = 1; % 晶格常數(shù)

N = 750; % 設(shè)置 N*N 格點(diǎn)

x = linspace(-N*L/2,(N-2)*L/2,N)'; % 1維位矢數(shù)組

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

NT = 50000; % 時(shí)間步長

TF = 50; %總時(shí)長

T = linspace(0,TF,NT); % 時(shí)間數(shù)組

dT = T(2)-T(1); % 時(shí)間步長

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

% % 設(shè)置位矢和動(dòng)量空間的坐標(biāo)

% 位矢相關(guān)

e1 = ones(N,1);

Rx = kron(x',e1); % 2維位矢數(shù)組，x坐標(biāo)

Ry = Rx'; % 2維位矢數(shù)組，y坐標(biāo)

Ra = sqrt(Rx.^2+Ry.^2);% 2維位矢數(shù)組，徑向坐標(biāo)

ThR = angle(Rx+1i*Ry);% 2維位矢數(shù)組，輻角坐標(biāo)

% 動(dòng)量相關(guān)

e1 = ones(N,1);

K = linspace(2*pi*(-1/2)/L,2*pi*(N/2-1)/(L*N),N); % 1維動(dòng)量數(shù)組

Kx = kron(K,e1); % 2維動(dòng)量數(shù)組，x

Ky = Kx'; % 2維動(dòng)量數(shù)組，y

% Ka = sqrt(Kx.^2+Ky.^2);

% ThK = angle(Kx+1i*Ky);

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

% 勢能設(shè)置

G = 1000; % 勢能放縮因子

d = 0.46*L*N; %方形勢能邊界

Hat = @(y) 4*((atan(y/25)+pi/2)/pi).^2; % 平滑的"階梯"函數(shù)

% 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分量

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

% 能帶參數(shù)設(shè)置，定義如下：

% 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);

%能帶（"動(dòng)能"）給出的時(shí)間演化算子

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;

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

% 設(shè)置 Gauss 波包

ao = 0.01*pi;% 動(dòng)量分布的標(biāo)準(zhǔn)差

% 中心動(dòng)量

kxo = 0;

kyo = 2.92;

% 中心位矢

xo = 0;

yo = -100;

% 依據(jù) H(k) 的本征態(tài)構(gòu)造波包

% 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)

% 位矢空間波函數(shù)，未歸一

psi0 = [reshape( fftshift(ifft2(fftshift(psi01))),N*N,1 )...

? ?;reshape( fftshift(ifft2( fftshift(psi02) )),N*N,1)];

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

%化為2維數(shù)組形式

psi21 = reshape(psi0(1:end/2),N,N);

psi22 = reshape(psi0(0.5*end+1:end),[N,N]);

% 繪制能帶，及中心動(dòng)量對應(yīng)能級的位置

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);% 為上能帶。若為下能帶，取相反數(shù)

scatter3(kxo,kyo,absh,1000,'r.')

hold off

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

% 時(shí)間演化以及繪圖

fig=figure;

fig.Position = [0 0 2000 1000];

ax = axes(fig);

h = 8e-5;% 繪圖的高度范圍放縮因子，依據(jù)模平方選取

xlim manual

ylim manual

zlim manual

% 自定義顏色函數(shù)

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);

% 設(shè)置顯示相位平面的位置

Z0 = -1.5*h*ones(N,N);

for t = 1:NT

? ?%時(shí)間演化

? ?%位矢表象，exp(-i U dt/2)

? ?psi21 = V11.*psi21 ;

? ?psi22 = V22.*psi22;

? ?%換為動(dòng)量表象

? ?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）以上模平方格點(diǎn)處的相位

? ? ? ?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

? ?% 其它繪圖參數(shù)

? ?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))