tic

%
%% Numerical Spectrum of the Hamiltonian for the Closed Dicke Model
% Hari Kumar Yadalam (hari.kumar@icts.res.in), Mahaveer Prasad (mahaveer.prasad@icts.res.in) and Manas Kulkarni (manas.kulkarni@icts.res.in)
% Version-1 (9th July 2021)
%
%% \hat{H}_{D}^{}= \omega_{c}^{} a_{}^{\dagger}a_{}^{} + \omega_{z}^{} S_{}^{z} + \frac{\lambda}{\sqrt{N}}\left[\left(a_{}^{\dagger}+a_{}^{}\right)\left(S_{}^{+} +S_{}^{-}\right)\right]
%% S_{}^{z}=b_{}^{\dagger}b_{}^{}-\frac{N}{2} \ \ \ S_{}^{+} = b_{}^{\dagger}\sqrt{N-b_{}^{\dagger}b_{}^{}} \ \ \  S_{}^{-} = \sqrt{N-b_{}^{\dagger}b_{}^{}} b_{}^{}
%% \lim_{N \to \infty}^{}\lambda_{c}^{}=\frac{\sqrt{\omega_{c}^{}\omega_{z}^{}}}{2}
%
%% Specify the System Parameters

omegac=1.0; % Energy of the photon mode
omegaz=1.0; % Energy of the Spin-N/2 energy
lambdacz=0.8; % Photon-Spin coupling strength

%
%% Specify the Spin size and Photon Cutoff

N=20; % 'N' for Spin-N/2
Ncut=500; % Photon number cutoff

%
%% Specify the Parity Sector

PARITY="EVEN"; % "EVEN" or "ODD" parity sector

%
%% Generate the Basis for the Parity Sector of the Hilbert Space

BASIS=[]; % Initialize the BASIS to empty array
counter=0; % Counter for the BASIS generation loop
switch PARITY % BASIS generation loop
    case "EVEN"
        for n=0:Ncut
            for s=0:N
                if mod((n+s),2)==0
                    counter=counter+1;
                    BASIS(counter,:)=[n,s];
                end
            end
        end
    case "ODD"
        for n=0:Ncut
            for s=0:N
                if mod((n+s),2)~=0
                    counter=counter+1;
                    BASIS(counter,:)=[n,s];
                end
            end
        end
end
clear n;
clear s;
clear counter;

%
%% Generate the Hamiltonian acting on the Parity Sector of the Hilbert Space

H=zeros([length(BASIS),length(BASIS)]); % Initialize the Hamiltonian to a null matrix
for basisindex2=1:length(BASIS) % Generate the matrix elements of the Hamiltonian and store them
    n2=BASIS(basisindex2,1);
    s2=BASIS(basisindex2,2);
    % Diagonal elements of Hamiltonain (Ignore the constant piece -omega0 S/2)
    %Coming from omega a^\dagger a + omega0 J^z
    basisindex1=basisindex2;
    H(basisindex1,basisindex2)=omegac*n2+omegaz*(s2);
    % Off-diagonal elements of Hamiltonian
    %Coming from (lambda/sqrt(S)) a^\dagger J^+
    n1=n2+1;
    s1=s2+1;
    basisindex1=ismember(BASIS,[n1,s1],'rows');
    H(basisindex1,basisindex2)=(lambdacz/sqrt(N))*(sqrt(n2+1)*sqrt(N-s2)*sqrt(s2+1));
    %Coming from (lambda/sqrt(S)) a J^-
    n1=n2-1;
    s1=s2-1;
    basisindex1=ismember(BASIS,[n1,s1],'rows');
    H(basisindex1,basisindex2)=(lambdacz/sqrt(N))*(sqrt(n2)*sqrt(N-s1)*sqrt(s2));
    %Coming from (lambda/sqrt(S)) a^\dagger J^-
    n1=n2+1;
    s1=s2-1;
    basisindex1=ismember(BASIS,[n1,s1],'rows');
    H(basisindex1,basisindex2)=(lambdacz/sqrt(N))*(sqrt(n2+1)*sqrt(N-s1)*sqrt(s2));
    %Coming from (lambda/sqrt(S)) a J^+
    n1=n2-1;
    s1=s2+1;
    basisindex1=ismember(BASIS,[n1,s1],'rows');
    H(basisindex1,basisindex2)=(lambdacz/sqrt(N))*(sqrt(n2)*sqrt(N-s2)*sqrt(s2+1));
end
clear BASIS;
clear n1;
clear s1;
clear n2;
clear s2;
clear basisindex1;
clear basisindex2;

%
%% Compute the Spectrum of the Hamiltonian in the Parity Sector

EVALS=eig(H); % Compute the spectrum
EVALS=sort(EVALS); % Sort the spectrum
clear H;

%
%% Plot the Spectrum and the Density of States

plot(EVALS,'r+') % Plot the spectrum
xlabel('k')
ylabel('E_{k}')
EVALS=EVALS(EVALS>50);
EVALS=EVALS(EVALS<350); % Extract flat part of the spectrum
histogram(EVALS,100,'Normalization','pdf') % Plot the density of states
xlabel('E')
ylabel('\rho(E)')

%
%% Compute the consecutive level spacing

EVALSPACE=diff(EVALS); % Compute the consecutive level spacing
EVALSPACE=sort(EVALSPACE)./mean(EVALSPACE); % Normalize to the mean

%
%% Plot the histogram of the consecutive level spacing

histogram(EVALSPACE,40,'Normalization','pdf') % Plot the histogram of the consecutive level spacing
hold on
s=linspace(0,max(EVALSPACE),100); % Points for generating data for Poisson and Wigner surmise
POISSON=exp(-s); % Poisson distribution for comparision
GOE=((pi/2).*s).*exp(-(pi/4).*(s.*s)); % Wigner surmise (GOE) for comparision
clear Npoints;
plot(s,POISSON,'r-','LineWidth',3) % Plot the Poisson distribution
hold on
plot(s,GOE,'b-','LineWidth',3) % Plot the Wigner surmise for GOE
hold off
xlabel('s')
ylabel('P(s)')
axis('tight');
legend('Dicke model','Poisson','Wigner (GOE)');

%

toc