用于模擬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 很小

% 高維度時(shí)，expm() 太慢太費(fèi)資源太危險(xiǎn)

% 替代方案：直接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);

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

% 制備光子的初態(tài)

%|n>

n1 = 100;

psi01 = vac;

for k=1:n1

? ? psi01 = (a')*psi01/(sqrt(k));

end

%coherent state |alpha>

% 為了簡(jiǎn)單，取實(shí)數(shù) alpha = sqrt(n2)

n2 = 100;

alpha = sqrt(n2);

psi02 = exp(-0.5*n2)*expm(alpha*a')*vac;

% psi0'*a'*a*psi0

% psi0'*psi0

% psi0 =?

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

% 設(shè)置完整的初態(tài)

psi1 = kron([1;0],psi0);?

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

% 計(jì)算 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;% 計(jì)算贗自旋的期望（*2）即可

y = psi1'*Sy*psi1;

z = psi1'*Sz*psi1;

% 替代方案：直接把psi1的密度矩陣約化到只剩二能級(jí)

% 計(jì)算效率慢

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