模拟Bloch波包用的部分代码
% 用于制备Bloch波包、观察Bloch波包运动的matlab代码
% 渣代码能力，轻喷
% 代码参考了 https://arxiv.org/pdf/0704.1622.pdf，有很多部分并非原创（英文注释基本上是原版代码里面复制过来的）
% 这个文献是在 https://www.zhihu.com/answer/278044905 里面翻到的
L = 80; % Interval Length
N = 4000; % No of points
x = linspace(-L/2,L/2,N)'; % Coordinate vector，粒子被限制在-L/2至L/2范围内，相当于一个无限深势阱
dx = x(2) - x(1); % Coordinate step
% Parameters for making intial wavefunction psi
%设置波包的中心动量
ko1 = 5.2;
ko2 = 4.7;
ko3 = 4;
a = 7; % 设置Gauss波包的初始宽度
% 设置波包的初始中心位置
xo1 = -20;
xo2 = -20;
xo3 = -20;
V=10; % 用于调整势垒高度
b = -0.1; % 设置外力大小和方向
% 由于空间自由度被离散化，可以准确描述的动量范围有上界。动量太大的后果其实也是一种 Bloch 振荡
if abs(ko1)>(0.5*pi*N/L)
error('too large k');
end
if abs(ko2)>(0.5*pi*N/L)
error('too large k');
end
if abs(ko3)>(0.5*pi*N/L)
error('too large k');
end
U1 = V*heaviside(x).*abs(rem( floor(x) ,2)); % 半边自由，半边 K-P 周期势能，可以用制备Bloch波包；如果要全局都是K-P,把heaviside删掉
% U2=U1+b*x; %有外力
% Finite-difference representation of Laplacian and Hamiltonian
e = ones(N,1); Lap = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2;
H1 = -(1/2)*Lap + spdiags(U1,0,N,N);
% H2 = -(1/2)*Lap + spdiags(U2,0,N,N);
% Parameters for computing psi(x,T) at different times 0 < T < TF
NT = 1600; % No. of time steps
TF = 16; T = linspace(0,TF,NT); % Time vector
dT = T(2)-T(1); % Time step
hbar = 1;
% Time displacement operator E=exp(-iHdT/hbar)
E1 = expm(-1i*full(H1)*dT/hbar); % time desplacement operator
% E2 = expm(-1i*full(H2)*dT/hbar);
% 初始波函数
psi1 = exp(-(x-xo1).^2/(2*a.^2)).*exp(1i*ko1*x);
psi1 = psi1/sqrt(psi1'*psi1*dx); % normalize the psi(x,0)
psi2 = exp(-(x-xo2).^2/(2*a.^2)).*exp(1i*ko2*x);
psi2 = psi2/sqrt(psi2'*psi2*dx);
psi3 = exp(-(x-xo3).^2/(2*a.^2)).*exp(1i*ko3*x);
psi3 = psi3/sqrt(psi3'*psi3*dx);
%对于已经制备好的Bloch波包psi0
% psi1 = psi0;
% psi1 = psi1/sqrt(psi1'*psi1*dx);
% psi2=psi1;
%***************************************************************
% Simulate rho(x,T) and plot for each T
%***************************************************************
h=0.3;
fig=figure;
fig.Position = [0 0 2000 1000];
for t = 1:NT % time index for loop
% calculate probability density rho(x,T)，相同的H
psi1 = E1*psi1; % calculate new psi from old psi
rho1 = conj(psi1).*psi1; % rho(x,T)
psi2 = E1*psi2; % calculate new psi from old psi
rho2 = conj(psi2).*psi2; % rho(x,T)
% psi3 = E1*psi3; % calculate new psi from old psi
% rho13 = conj(psi3).*psi3; % rho(x,T)
% 如果要用不同的H
% psi2 = E2*psi2;
% rho2 = conj(psi2).*psi2;
% psi3 = E*psi3;
% rho3 = conj(psi3).*psi3;
% f1 = fft(psi1);
% f2 = fft(psi2);
% f3 = fft(psi3);
subplot(211);
plot(x,(h/2)*U1/(2*V));% 画势能
hold on
% plot(x,rho1,'b',x,rho2,'k',x,rho3,'m'); % plot rho(x,T) vs. x
% plot(x,rho1,'b',x,real(psi1).^2,'c--',x,rho2,'k',x,real(psi2).^2,'b--',x,rho3,'r',x,real(psi3).^2,'g--');
plot(x,rho1,'b',x,real(psi1).^2,'c--');
axis([-L/2 L/2 0 h]); % set x,y axis parameters for plotting
xlabel('x', 'FontSize', 24);
hold off
subplot(212);
plot(x,(h/2)*(U1)/(2*V)+0.1);
% plot(x,(h/2)*(U2)/(2*V)+0.1);
hold on
plot(x,rho2,'k',x,real(psi2).^2,'b--');
axis([-L/2 L/2 0 h]); % set x,y axis parameters for plotting
xlabel('x', 'FontSize', 24);
hold off
pause(0.001);
end
