% Atom-field interaction: vacuum Rabi oscillations (Sec. V)
% Jaynes-Cummings Hamiltonian; the atom starts excited with n photons in the
% cavity and periodically exchanges a photon with the field.
%   Pe(t) = cos^2(g t sqrt(n+1)),  Pg(t) = sin^2(g t sqrt(n+1))

clear;          % Clear memory
clc;            % Clear the command window/screen
d = 10;         % 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);             % atomic Hamiltonian
Hfield = hbar*Wf*kron(I_a,Ad*A);                % field Hamiltonian
Hint = hbar*g*(kron(Splus,A)+kron(Sminus,Ad));  % interaction Hamiltonian
H = Hatom + Hfield + Hint;  % JC Hamiltonian
n = 4;                      % initial number of photons in the cavity
en = kron(es,I_f(:,n+1));   % atom excited, field has n photons   |e,n>
gn = kron(gs,I_f(:,n+2));   % atom ground, field has n+1 photons  |g,n+1>
Psi = en;                   % initial state
dt = 0.1;                   % time step
U = expm(-i*H*dt/hbar);     % unitary time-evolution operator
T = 0:dt:30;                % total evolution time
for t = 1:length(T)
  Pe(t) = norm(en'*Psi)^2;  % probability of |e,n>
  Pg(t) = norm(gn'*Psi)^2;  % probability of |g,n+1>
  Psi = U*Psi;              % time-evolved state
  Psi = Psi/norm(Psi);      % renormalize
end
plot(T, Pe, 'r', T, Pg, 'k')        % Pe in red, Pg in black
title('Vacuum Rabi oscillations (n = 4, g = 0.1)')
xlabel('t'); ylabel('Probability'); legend('P_e', 'P_g')