% Zero time-delay second-order coherence function g^(2)(0) (Sec. IV.C)
%   g2(0) = <a'^2 a^2> / <a'a>^2
% Number state |n>: 1 - 1/n (sub-Poissonian); coherent: 1; thermal: 2.

clear;          % Clear memory
clc;            % Clear the command window/screen
d = 25;         % dimension of the field
I = eye(d);     % identity matrix
A = diag(sqrt(1:d-1), 1);   % Annihilation operator
Ad = A';                    % Creation operator
AdA = Ad*A;                 % Number operator
Ad2A2 = Ad*Ad*A*A;          % a'^2 a^2

%% Number state |4>
Ket4 = I(:,5);                          % four-photon state |4>
AdANumber = Ket4'*AdA*Ket4;
Ad2A2Number = Ket4'*Ad2A2*Ket4;
g2Number = Ad2A2Number/(AdANumber)^2    % g2(0) for number state (= 1 - 1/4 = 0.75)

%% Coherent state
alpha = sqrt(3);    % amplitude of the coherent state
Coh = 0;            % initialization
for x = 0:d-1
  Coh = Coh + exp(-norm(alpha)^2/2)*alpha^x/sqrt(prod(1:x))*I(:,x+1);
end
AdACoherent = Coh'*AdA*Coh;
Ad2A2Coherent = Coh'*Ad2A2*Coh;
g2Coherent = Ad2A2Coherent/(AdACoherent)^2   % g2(0) for coherent state (= 1)

%%% Thermal state
nth = 0.85;     % average number of photons in the thermal state
RhoTh = 0;
for x = 0:d-1
  RhoTh = RhoTh + nth^(x)/(1+nth)^(x+1)*I(:,x+1)*I(:,x+1)';
end
AdAThermal = trace(AdA*RhoTh);
Ad2A2Thermal = trace(Ad2A2*RhoTh);
g2Thermal = Ad2A2Thermal/(AdAThermal)^2       % g2(0) for thermal state (= 2)