關(guān)于腔中自發(fā)輻射的效率低下的不知道還有沒有bug的渣渣代碼：

% set photon degrees

L = 100*pi/100; % length of cavity

xc = linspace(-L/2,L/2,400)';

Dw = pi/L;

w0 =51*Dw; % centre of freq, odd integer*pi/L

Nf = 2*2+1; % number of freq taken into account

IM = 2; % largest phonton number of single freq

ptFr = (w0-Dw*(Nf-1):2*Dw:w0+Dw*(Nf-1))'; % phonton frequencies of interest

Dim_p = (IM+1)^Nf;

pta = sqrt(linspace(0,IM,IM+1)');

a_s = spdiags(pta,1,IM+1,IM+1);% photon annihilation operator for single freq

a = sparse(Dim_p,Dim_p*Nf);

for j=1:Nf

?I1 = speye((IM+1).^(j-1));

?I2 = speye((IM+1).^(Nf-j));

?a0 = kron(I1,a_s);

?a0 = kron(a0,I2);

?a(:,(j-1)*Dim_p+1:(j)*Dim_p) = a0;

end

% a(:,(m-1)*(IM+1).^Nf+1:(m)*(IM+1).^Nf) is

% the annihilation operator for the m'th mode

% Hamiltonian and vacuum state

% of free photon in chosen modes

Hp = sparse(Dim_p,Dim_p);

e1 = zeros(IM+1,1);

e1(1) = 1;

vac = 1;

for k=1:Nf

? vac = kron(vac, e1); % vacuum state of photon

? a1 = a(:,(k-1)*Dim_p+1:(k)*Dim_p); % a of k'th freq

? Hp = ?Hp+ptFr(k)*(a1'*a1); % calculate Hamiltonian of free photon

end

% 2 level Hamiltonian without coupling

DE = w0;

H2 = spdiags([0;DE],0,2,2);

% full-system Hamiltonian without coupling

I2 = speye(2,2);

Ip = speye(Dim_p,Dim_p);

H0 = kron(H2,Ip)+kron(I2,Hp);

% coupling

g = 0.4/7; % coupling coefficient

sx = spdiags([[1;1] [1;1]],[-1 1],2,2);% Pauli x matrix

HI = sparse(2*Dim_p,2*Dim_p);

for j=1:Nf

? a1 = a(:,(j-1)*Dim_p+1:(j)*Dim_p);

? HI = HI+g*sqrt(ptFr(j))*kron(sx,a1+a1'); % calculate Interaction term in H

end

% full Hamiltonian

H = H0+HI;

% time evolution

NT = 5001; % No. of time steps

TF = 4; % 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*full(H)*dT);

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

% 可能用到的算子

sz = spdiags([1;-1],0,2,2);% Pauli z matrix

Ip = speye(Dim_p,Dim_p); % 光子空間的恒等算子

Sz = kron(sz,Ip);

Itt = kron(I2,Ip); % 全空間的恒等算子

for i=1:Nf

? % 湮滅算子

? eval(['a',num2str(i),'=',...

? ? ? 'kron(I2,a(:,(i-1)*Dim_p+1:(i)*Dim_p))',';']);

? % 粒子數(shù)算子

? eval(['n',num2str(i),'=',...

? ? ? 'a',num2str(i),'''*a',num2str(i),';']);

end

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

%初態(tài)

psi1 = kron([0;1],vac);

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

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

%時(shí)間演化

Nx = length(xc);

In = zeros(Nx,1);

Pop1 = zeros(NT,1);

Pop1(1) = TwoLvPop(psi1,Sz);

for i=1:Nf

? ?eval(['Num',num2str(i),'=',...

? ? ? ?'zeros(NT,1) ;']);

? ?eval(['Num',num2str(i),'(1) =real(',...

? ? ? 'psi1''*n',num2str(i), '*psi1 );']);

end

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

? psi1 = E*psi1; % time evolution

? % 按需修改后面的東西

? Pop1(t+1) = TwoLvPop(psi1,Sz); % population of 激發(fā)態(tài)

? for i=1:Nf % 各個(gè)模式的光子數(shù)期望

? ? ? eval(['Num',num2str(i),'(t+1) =real(',...

? ? ? ? ? 'psi1''*n',num2str(i), '*psi1 );']);

? end

? % 光強(qiáng)期望

? parfor j=1:Nx

? ? ? In(j) = ElcAmp(psi1,xc(j),a,IM,ptFr);

? end

end

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

function z = TwoLvPop(psi1,Sz)

? z = psi1'*Sz*psi1;

? z = (1-z)/2;

end

function ElecI = ElcAmp(S1,x1,a,IM,ptFr)

? Nf = length(ptFr);

? dimp = (IM+1).^Nf;

? El = sparse(dimp*2,dimp*2);

? I2 = speye(2,2);

? for k=1:Nf

? ? ? wk = ptFr(k);

? ? ? a1 = kron(I2,a(:,(k-1)*dimp+1:(k)*dimp));

? ? ? El = El+sqrt(wk)...

? ? ? ? ? *cos(wk*x1)...

? ? ? ? ? *(a1);

? end

? ElecI = real((S1'*(El'*El)*S1));

?

end