% Coupled cavities: energy exchange between two field modes (Sec. VI.A)
% One photon starts in cavity 1, cavity 2 is empty. The photon is
% periodically transferred between the two cavities.
%   P10(t) = cos^2(Jt),  P01(t) = sin^2(Jt)

clear;          % Clear memory
clc;            % Clear the command window/screen
d = 10;         % dimension of the cavity field
hbar = 1;
W1 = 1;         % resonance frequency of the first cavity
W2 = 1;         % resonance frequency of the second cavity
J = 0.1;        % inter-cavity coupling constant
A = diag(sqrt(1:d-1), 1);   % Annihilation operator
Ad = A';                    % Creation operator
I_f = eye(d);               % identity operator for the fields
H = hbar*W1*kron(Ad*A,I_f) + hbar*W2*kron(I_f,Ad*A) + hbar*J*(kron(Ad,A)+kron(A,Ad));
S10 = kron(I_f(:,2),I_f(:,1));  % cavity 1 has |1>, cavity 2 has |0>
S01 = kron(I_f(:,1),I_f(:,2));  % cavity 1 has |0>, cavity 2 has |1>
Psi = S10;                  % initial state
dt = 0.1;                   % time step
U = expm(-i*H*dt/hbar);     % unitary time-evolution operator
T = 0:dt:50;                % total evolution time
for t = 1:length(T)
  P10(t) = norm(S10'*Psi)^2;    % prob. photon in cavity 1
  P01(t) = norm(S01'*Psi)^2;    % prob. photon in cavity 2
  Psi = U*Psi;                  % time evolution
  Psi = Psi/norm(Psi);          % normalization
end
plot(T, P10, 'r', T, P01, 'k')      % cavity 1 red, cavity 2 black
title('Photon exchange between coupled cavities (J = 0.1)')
xlabel('t'); ylabel('Probability'); legend('P_{10}', 'P_{01}')