自发辐射相关代码

关于腔中自发辐射的效率低下的不知道还有没有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))',';']);

  % 粒子数算子

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

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

end

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

%初态

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

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

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

%时间演化

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 激发态

  for i=1:Nf % 各个模式的光子数期望

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

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

  end

  % 光强期望

  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