模拟Rabi振荡的部分代码
用于模拟Rabi振荡的部分代码：
% set photon parameters
w0 =4000; % freq of photon
IM = 1200; % largest phonton number
ptNum = linspace(0,IM,IM+1)';
ptNumM = spdiags(ptNum,0,IM+1,IM+1);% photon number operator
pta = sqrt(linspace(0,IM,IM+1)');
a = spdiags(pta,1,IM+1,IM+1);% photon annihilation operator
ad = a';% photon creation operator
Inum = speye(IM+1,IM+1);% identical operator
% Hamiltonian of free photon in chosen modes
Hp = w0*ptNumM;
% imagesc(full(Hp)),colorbar
% vacuum state
vac = zeros(IM+1,1);
vac(1) = 1; % vacuum state of photon
% 2 level Hamiltonian without coupling
DE = 3500;
H2 = spdiags([0;DE],0,2,2);
% full Hamiltonian without coupling
I2 = speye(2,2);
Ip = speye((IM+1),(IM+1));
H0 = kron(H2,Ip)+kron(I2,Hp);
% coupling
g = 50; % coupling coefficient
sx = spdiags([[1;1] [1;1]],[-1 1],2,2);% Pauli x matrix
HI = g*kron(sx,a+a'); % calculate Interaction term in H
% full Hamiltonian
H = H0+HI;
% spy(H)
% imagesc(full(H),[-100 100]),colorbar,colormap gray
% time evolution
NT = 16000; % No. of time steps
TF = 1; % total time
T = linspace(0,TF,NT); % Time vector
dT = T(2)-T(1); % Time step
% Time displacement operator E=exp(-iHdT/hbar)
E = expm(-1i*H*dT);
% 用 expm() 不必令 dT 很小
% 高维度时，expm() 太慢太费资源太危险
% 替代方案：直接Taylor展开到前几阶，但是 dT 最好比较小
% I = speye(size(H));
% A = -1i*H*dT;
% E = I+A+(A^2/2)+(A^3/6)+(A^4/24)+(A^5/120)+(A^6/720);
% *********************
% 制备光子的初态
%|n>
n1 = 100;
psi01 = vac;
for k=1:n1
psi01 = (a')*psi01/(sqrt(k));
end
%coherent state |alpha>
% 为了简单，取实数 alpha = sqrt(n2)
n2 = 100;
alpha = sqrt(n2);
psi02 = exp(-0.5*n2)*expm(alpha*a')*vac;
% psi0'*a'*a*psi0
% psi0'*psi0
% psi0 =
% *********************
% 设置完整的初态
psi1 = kron([1;0],psi0);
% *********************
% 计算 Bloch 球
sx = spdiags([[1;1] [1;1]],[-1 1],2,2);% Pauli x matrix
sy = spdiags([[1i;1i] [-1i;-1i]],[-1 1],2,2);% Pauli y matrix
sz = spdiags([[1;-1]],0,2,2);% Pauli z matrix
Ip = speye((IM+1),(IM+1));
Sx = kron(sx,Ip);
Sy = kron(sy,Ip);
Sz = kron(sz,Ip);
x = psi1'*Sx*psi1;% 计算赝自旋的期望（*2）即可
y = psi1'*Sy*psi1;
z = psi1'*Sz*psi1;
% 替代方案：直接把psi1的密度矩阵约化到只剩二能级
% 计算效率慢
% P = psi1*psi1';
% P00 = P(1:end/2,1:end/2);
% P01 = P(1:end/2,1+end/2:end);
% P10 = P(1+end/2:end,1:end/2);
% P11 = P(1+end/2:end,1+end/2:end);
% p = [trace(P00),trace(P01);trace(P10),trace(P11)];
% sx = [0 1;1 0];
% sy = [0 -1i;1i 0];
% sz = [1 0;0 -1];
% x = trace(p*sx);
% y = trace(p*sy);
% z = trace(p*sz);
