三能級與光子相關。

• Hamiltonian 以及演化算子

% set photon parameters

w0 = 5500; % freq of photon 0

w1 = 3500; % freq of photon 0

IM = 40; % largest phonton number

% operators for a single photon mode

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

a0 = kron(a,Inum); % annihilation operator for w0

a1 = kron(Inum,a);% annihilation operator for w1

% Hamiltonian of free photon in chosen modes

Hp = w0*kron(ptNumM,Inum)+w1*kron(Inum,ptNumM);

% vacuum state

vac00 = zeros(IM+1,1);

vac00(1) = 1;

vac = kron(vac00,vac00);% vacuum state of photons

% 3 level Hamiltonian without coupling

% level 1 2 3

H3 = spdiags([0;2000;7000],0,3,3);

% full space's Hamiltonian without coupling

I3 = speye(3,3);

Ip = speye((IM+1)^2,(IM+1)^2);

H0 = kron(H3,Ip)+kron(I3,Hp);

% coupling

g0 = 50; % coupling coefficient between photon 0 and state 1,3

g1 = 50; % coupling coefficient between photon 1 and state 2,3

sig0 = sparse(3,3);

sig0(1,3) = 1;sig0(3,1) = 1;

sig1 = sparse(3,3);

sig1(2,3) = 1;sig1(3,2) = 1;

HI = g0*kron(sig0,a0+a0')+g1*kron(sig1,a1+a1'); % 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 最好比較小

• 初態(tài)（僅給出了張量積/直積態(tài)的情形）

% photons

n0 = 25;

n1 = 0;

psi01 = vac;

psi02 = vac;

NorA = 0; % num (0) or coherent (1)

if ~NorA? ?

? ? for k=1:n0

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

? ? end

? ? for k=1:n1

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

? ? end

else? ??

? ? alpha0 = sqrt(n0);

? ? psi01 = exp(-0.5*n0)*expm(alpha0*a0')*vac;

? ? alpha1 = sqrt(n1);

? ? psi01 = exp(-0.5*n1)*expm(alpha1*a1')*psi01;

end

% 3 levels

S = [1;0;0];

S = S/sqrt(S'*S);

% full state vec

psi0 = kron(S,psi01);

• 演化、相關量的計算及畫圖的簡單示例

FO = [ones(6,1);zeros(12,1)];

SO = [zeros(6,1);ones(6,1);zeros(6,1)];

TO = [zeros(12,1);ones(6,1)];

for t = 1:NT? % time index for loop

? ? psi0 = E*psi0; % 演化

????

????rho3 = Rho3(psi0,IM); % 約化密度矩陣，函數(shù)見后面

? ? phase = angle(rho3); % 提取相位

? ? Purity = trace(rho3*rho3); % 粗略反映混態(tài)的程度

? ? B = bar3(ax1,abs(rho3)); % 畫圖

%著色

? ? B(1).CData = kron([1,1,1,1],...

? ? ? ? phase(1,1)*FO+phase(1,2)*SO+phase(1,3)*TO);

? ? B(2).CData = kron([1,1,1,1],...

? ? ? ? phase(2,1)*FO+phase(2,2)*SO+phase(2,3)*TO);

? ? B(3).CData = kron([1,1,1,1],...

? ? ? ? phase(3,1)*FO+phase(3,2)*SO+phase(3,3)*TO);

? ? set(gca,'CLim',[-pi,pi]);

????colormap(ax1,hsv);

? ? set(gca,'zlim',[-0.01,1.1]);

? ? daspect(gca,[1 1 0.4]);

? ??

????drawnow

end

function rho3=Rho3(psi,IM) % 計算 3 level? 約化密度矩陣

Ip = speye((IM+1)^2,(IM+1)^2);

rho3 = zeros(3,3);

for j=1:3

? ? for i=1:j

? ? ? ? sig_ij = sparse(3,3);

? ? ? ? sig_ij(i,j) = 1;

? ? ? ? rho3(i,j) = psi'*(kron(sig_ij,Ip))*psi;

? ? ? ? rho3(j,i) = rho3(i,j)';

? ? end

end

end