% 用于制備Bloch波包、觀察Bloch波包運(yùn)動(dòng)的matlab代碼

% 渣代碼能力，輕噴

% 代碼參考了 https://arxiv.org/pdf/0704.1622.pdf，有很多部分并非原創(chuàng)（英文注釋基本上是原版代碼里面復(fù)制過來的）

% 這個(gè)文獻(xiàn)是在 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范圍內(nèi)，相當(dāng)于一個(gè)無限深勢(shì)阱

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

% Parameters for making intial wavefunction psi

%設(shè)置波包的中心動(dòng)量

?ko1 = 5.2;?

?ko2 = 4.7;

?ko3 = 4;

?

a = 7; % 設(shè)置Gauss波包的初始寬度

% 設(shè)置波包的初始中心位置

xo1 = -20;

xo2 = -20;

xo3 = -20;

V=10; % 用于調(diào)整勢(shì)壘高度

b = -0.1; % 設(shè)置外力大小和方向

?% 由于空間自由度被離散化，可以準(zhǔn)確描述的動(dòng)量范圍有上界。動(dòng)量太大的后果其實(shí)也是一種 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 周期勢(shì)能，可以用制備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);

% 初始波函數(shù)

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

%對(duì)于已經(jīng)制備好的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));% 畫勢(shì)能

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