%% Numerical Solution of Driven Dissipative Jaynes-Cummings Model
% Hari Kumar Yadalam (hari.kumar@icts.res.in) and Manas Kulkarni (manas.kulkarni@icts.res.in)
% Version-1 (6th July 2021)
%
% $\frac{d}{dt}\hat{\rho}(t)=-i\left[\hat{H}_{JC}^{}+\hat{H}_{D}^{},\hat{\rho}(t)\right]-\frac{\kappa}{2}(1+n(\omega_c,T))\left[\left\{a_{}^{\dagger}a_{}^{},\hat{\rho}(t)\right\}-2 a_{}^{} \hat{\rho}(t) a_{}^{\dagger}\right]-\frac{\kappa}{2}n(\omega_c,T)\left[\left\{a_{}^{}a_{}^{\dagger},\hat{\rho}(t)\right\}-2 a_{}^{\dagger} \hat{\rho}(t) a_{}^{}\right]-\frac{\gamma}{2}(1+n(\omega_q,T))\left[\left\{\sigma_{}^{+}\sigma_{}^{-},\hat{\rho}(t)\right\}-2 \sigma_{}^{-} \hat{\rho}(t) \sigma_{}^{+}\right]-\frac{\gamma}{2}n(\omega_q,T)\left[\left\{\sigma_{}^{-}\sigma_{}^{+},\hat{\rho}(t)\right\}-2 \sigma_{}^{+} \hat{\rho}(t) \sigma_{}^{-}\right]]-\frac{\gamma_{\phi}}{4}\left[\left\{\sigma_{}^{z}\sigma_{}^{z},\hat{\rho}(t)\right\}-2 \sigma_{}^{z} \hat{\rho}(t) \sigma_{}^{z}\right]$
% $\hat{H}_{JC}^{}=\tilde{\omega}_{c}^{}a_{}^{\dagger}a_{}^{} + \frac{\tilde{\omega}_{q}^{}}{2} \sigma_{}^{z} + g\left[a_{}^{\dagger}\sigma_{}^{-}+\sigma_{}^{+} a_{}^{}\right]  \ \ \ \ \ \ \ \ \    \hat{H}_{D}^{} = \mu_{}^{}\left[a_{}^{\dagger} + a_{}^{}\right]$
%
%% Specify the System parameters

omegac=1.0; % Cavity frequency - Drive frequency
omegaq=1.0; % Qubit frequency - Drive frequency
g=0.1; % Cavity-Qubit coupling
kappac=0.01; % Cavity dissipation rate
gammaq=0.01; % Qubit dissipation rate
gammaphi=0.01; % Qubit dephasing rate
mu=0.01; % Cavity drive strength
Tres=0.1; % Reservoir temperature
nresc=1/(exp(omegac/Tres)-1); % Average photon number of cavity if it is at thermal equilibrium with temperature Tres
nresq=1/(exp(omegaq/Tres)-1); % Average magnetization of qubit if it is at thermal equilibrium with temperature Tres

%
%% Specify the Photon Number Cutoff

Ncut=25; % Photon cutoff

%
%% Construct Photon and Qubit Operators in the individual Hilbert spaces

a=zeros([Ncut+1,Ncut+1]);
ip=eye([Ncut+1,Ncut+1]);
for k=1:Ncut
    a(k,k+1)=sqrt(k);
end
adag=a';
sz=[1,0;0,-1];
sm=[0,0;1,0];
iq=eye([2,2]);
sp=sm';

%
%% Lift the Operators onto the Composite Hilbert Space

Sz=kron(ip,sz);
Sm=kron(ip,sm);
Sp=kron(ip,sp);
A=kron(a,iq);
Adag=kron(adag,iq);
Id=kron(ip,iq);
clear a;
clear adag;
clear ip;
clear sz;
clear sm;
clear sp;
clear iq;

%
%% Construct Hamiltonian

Hjc=omegac.*(Adag*A)+(omegaq/2).*Sz+g.*(Adag*Sm+Sp*A);
Hd=mu.*(Adag+A);

%
%% Construct the vectorized Liouvillian

LvN=-1i.*(kron(Id,Hjc+Hd)-kron((Hjc+Hd).',Id))-((kappac/2)*(1+nresc)).*(kron(Id,Adag*A)+kron((Adag*A).',Id)-2.*kron(Adag.',A))-((kappac/2)*nresc).*(kron(Id,A*Adag)+kron((A*Adag).',Id)-2.*kron(A.',Adag))-((gammaq/2)*(1+nresq)).*(kron(Id,Sp*Sm)+kron((Sp*Sm).',Id)-2.*kron(Sp.',Sm))-((gammaq/2)*nresq).*(kron(Id,Sm*Sp)+kron((Sm*Sp).',Id)-2.*kron(Sm.',Sp))-(gammaphi/4).*(kron(Id,Sz*Sz)+kron((Sz*Sz).',Id)-2.*kron(Sz.',Sz));
clear Hjc;
clear Hd;
clear Id;

%
%% Specify the parameters for Dynamics

Nsteps=1000; % Number of time steps
deltaT=(1/kappac)/100; % Time step length
Tmax=100.*deltaT; % Maximum time
TIME=linspace(0,Tmax,Nsteps+1);
NPHOTONS=zeros([Nsteps+1,1]);
SZ=zeros([Nsteps+1,1]);
clear Tmax;

%
%% Specify Initial State

Psi0p=zeros([Ncut+1,1]);
Psi0p(1,1)=1;
Psi0q=[0;1];
Psi0=kron(Psi0p,Psi0q); % INITIAL STATE IS A PRODUCT STATE
Rho0=Psi0*Psi0'; % INITIAL STATE IS A PURE STATE
Rho0=reshape(Rho0,[(2*(Ncut+1))^2,1]); % VECTORIZING DENSITY MATRIX
clear Psi0p;
clear Psi0q;
clear Psi0;

%
%% Dynamics

RhoT=Rho0; % TIME ZERO DENSITY MATRIX
clear Rho0;
RhoT=(1/trace(reshape(RhoT,[2*(Ncut+1),2*(Ncut+1)]))).*RhoT;
NPHOTONS(1)=trace(Adag*A*reshape(RhoT,[2*(Ncut+1),2*(Ncut+1)])); % DEVECTORIZING DENSITY MATRIX AND COMPUTE AVERAGE
SZ(1)=trace(Sz*reshape(RhoT,[2*(Ncut+1),2*(Ncut+1)])); % DEVECTORIZING DENSITY MATRIX AND COMPUTE AVERAGE
GdeltaT=expm(deltaT.*LvN); % PROPOGATOR FOR EACH TIME STEP
clear deltaT;
clear LvN;
for k=1:Nsteps
    RhoT=GdeltaT*RhoT;
    RhoT=(1/trace(reshape(RhoT,[2*(Ncut+1),2*(Ncut+1)]))).*RhoT;
    NPHOTONS(k+1)=trace(Adag*A*reshape(RhoT,[2*(Ncut+1),2*(Ncut+1)]));
    SZ(k+1)=trace(Sz*reshape(RhoT,[2*(Ncut+1),2*(Ncut+1)]));
end
clear GdeltaT;
clear RhoT;
clear A;
clear Adag;
clear Sz;
clear Sm;
clear Sp;
clear k;
clear Nsteps;

%
%% Plot the results

plot(kappac.*TIME,real(NPHOTONS),'r-','LineWidth',3)
xlabel('\kappa t')
ylabel('< a^{T}(t) a(t) >')
axis("tight")
plot(kappac.*TIME,real(SZ),'r-','LineWidth',3)
xlabel('\kappa t')
ylabel('< \sigma^z(t) >')
axis("tight")

%

