Feshbach共振Hamiltonian的构造和初态的设置

模拟Feshbach的视频中用到的部分代码：

% 设置格点化的空间坐标

L = 50; % Interval Length

N = 2000; % No of points

x = linspace(0.5,0.5+L,N)'; % Coordinate vector

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

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

%设置 Hamiltonian

e = ones(N,1);

e2 = [1;0];% open channel对应的（二能级）态矢量

e3 = [0;1];% close channel对应的（二能级）态矢量

e4 = [1;1];

V = 10;% 控制相互作用的强度

Lap = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2; % 1维Laplacian

Uclose = 1e3./(x.^12)-1e3./(x.^6)+ 107.965; % close channel 的势能

Uopen = (1e3./(x.^12)-1.5e3./(x.^6));% open channel 的势能

C = V*(pi/2-atan(x-1.5));% 在左侧加上channel间的耦合

G =0.5;% 控制边界上的衰减

U1 = -(1i*G).*(atan(0.5*(x-(L)))+pi/2).*(1./(1+(3e-2)*(L-x).^2));

% 在右侧缓缓加上虚势能，让反射的低能量的波包也可以被移出系统

% plot(x,Uopen,x,Uclose)

% axis([0.5 10 -150 105])

% Hamiltonian

H1 = -(1/2)*kron(Lap,speye(2)) +...      % 动能

    kron(spdiags(Uclose,0,N,N),spdiags(e3,0,2,2))+... %close channel potential

    kron(spdiags(Uopen,0,N,N),spdiags(e2,0,2,2))+... %open channel potential

    kron(spdiags(C,0,N,N),spdiags([e4 e4],[-1 1],2,2))+... % coupling between channels

    kron(spdiags(U1,0,N,N),speye(2)); % 虚势能

% Parameters for computing psi(x,T) at different times 0 < T < TF

NT = 1000; % No. of time steps

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

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

hbar = 1;

% Time displacement operator E=exp(-iHdT/hbar)

E = expm(-1i*full(H1)*dT/hbar); % time desplacement operator

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

%设置初态

ko1 = -0.5; % Peak momentum

a = 4; % package width

xo=15;% peak center

% 离散化的空间坐标对最大动量有限制，不然会出现Bloch振荡这类现象

if ko1>(pi*N/L)

    error('too large ko1');

end

S1 = [1;0];

S1 = S1/sqrt(S1'*S1);% 设置初态channel

% 设置空间波函数为 平面波*Gauss轮廓

psi1 = exp(-(x-xo).^2/(2*a.^2)).*exp(1i*x*ko1);

psi1 = psi1/sqrt(psi1'*psi1*dx); % normalize the psi(x,0)

% 把两个通道的自由度和空间自由度的态矢量张量积起来

psi1 = kron(psi1,S1);