% Atomic inversion with a coherent field: collapse and revival (Sec. V)
% Atom excited, field in a coherent state. The atomic inversion <sigma_z(t)>
% collapses and later revives -- a hallmark of field quantization.
%
% Note: d = 50 and a long evolution make this the heaviest example here;
% give it a few seconds.

clear;          % Clear memory
clc;            % Clear the command window/screen
d = 50;         % dimension of the cavity field
hbar = 1;
W0 = 1;         % atomic frequency
Wf = 1;         % cavity field frequency
g = 0.1;        % coupling constant
A = diag(sqrt(1:d-1), 1);   % Annihilation operator
Ad = A';                    % Creation operator
Sz = [1,0;0,-1];            % sigma_z
Splus = [0,1;0,0];          % sigma_+
Sminus = [0,0;1,0];         % sigma_-
gs = [0;1];                 % ground state
es = [1;0];                 % excited state
I_a = eye(2);               % identity operator for the atom
I_f = eye(d);               % identity operator for the field
Hatom = (1/2)*hbar*W0*kron(Sz,I_f);
Hfield = hbar*Wf*kron(I_a,Ad*A);
Hint = hbar*g*(kron(Splus,A)+kron(Sminus,Ad));
H = Hatom + Hfield + Hint;  % JC Hamiltonian
alpha = 3;                  % coherent state amplitude
Coh = 0;
for x = 0:d-1
  Coh = Coh + exp(-norm(alpha)^2/2)*alpha^x/sqrt(prod(1:x))*I_f(:,x+1);
end
Psi = kron(es,Coh);         % initial state: atom |e>, field coherent
dt = 0.1;                   % time step
U = expm(-j*H*dt);          % unitary time-evolution operator
T = 0:dt:500;               % total evolution time
for t = 1:length(T)
  W(t) = Psi'*kron(Sz,I_f)*Psi;   % <sigma_z> (atomic inversion)
  Psi = U*Psi;
  Psi = Psi/norm(Psi);
end
plot(T, W)
title('Atomic inversion <\sigma_z(t)>: collapse and revival (\alpha = 3)')
xlabel('t'); ylabel('<\sigma_z>')