c33ae98Companion to 'Coding closed and open quantum systems in MATLAB' (arXiv:1911.04906)Jeremy Magland 1Omega = 1; % Raby frequency from |1> to |2>
2gamma0 = 0.2*Omega; % Decay rate
3dim = 2; % Dimension Hilbert space two-level system
4Is = eye(dim); % Identity matrix Hilbert space
5v1 = Is(:,1); % Excited state for the atom
6v2 = Is(:,2); % Ground state fro the atom
7s11 = v1*v1'; % Atom operator sigma_{11}
8s22 = v2*v2'; % Atom operator sigma_{22}
9sp = v1*v2'; % Atom operator sigma_{+}
10sm = v2*v1'; % Atom operator sigma_{-}
11HL = -0.5*Omega*(sp+sm); % Hamiltonian of the two-level system
12N = 0; % Mean number of photons at zero temperature
13Lrad = gamma0*(N+1)*(kron(conj(sm),sm)-0.5*kron(Is,sm'*sm)-0.5*kron(sm.'*conj(sm),Is)) + ...
14 gamma0*N*(kron(conj(sp),sp)-0.5*kron(Is,sp'*sp)-0.5*kron(sp.'*conj(sp),Is));
15L = -1i*kron(Is,HL)+1i*kron(HL.',Is)+Lrad; % Lindblad operator
16TOL = 1e-6;
17[R_sort,L_sort,lambda_sort] = sortingEigenvalues(dim,TOL,L);
18psi_0 = v2; % Ground state
19rho_0 = psi_0*psi_0'; % Initial density matrix
20Nt = 100000; % Number of steps for time
21ti = 0; % Initial time
22tf = 50/Omega; % Final time
23dt = (tf-ti)/(Nt-1); % Step time dt
24t = ti:dt:tf; % Time vector
25rho11 = zeros(size(t)); % Matrix elements \rho_{ij} = <i|\rho|j>
26sigmap = zeros(size(t));
27for n=1:Nt % General solution
28 rho = zeros(dim,dim);
29 for k=1:length(lambda_sort)
30 Lk = L_sort{k};
31 Rk = R_sort{k};
32 ck = trace(rho_0*Lk);
33 rho = rho + ck*exp(lambda_sort(k)*t(n))*Rk;
34 end
35 rho11(n) = rho(1,1);
36 sigmap(n) = trace(rho*sp);
37end
38% Exact Solution at zero temperature
39% p_e: population excited state, sp = <\sigma_x>
40mu = 1i*sqrt((gamma0/4)^2-Omega^2);
41pe_exact = Omega^2/(gamma0^2+2*Omega^2)*(1-exp(-3*gamma0*t/4).*(cos(mu*t)+3*gamma0/4/mu*sin(mu*t)));
42if gamma0~=0
43 sp_exact = -1i*Omega*gamma0/(gamma0^2+2*Omega^2)*(1-exp(-3*gamma0*t/4).*(cos(mu*t)+(gamma0/4/mu)-Omega^2/gamma0/mu*sin(mu*t)));
44else
45 sp_exact = -1i/2*sin(Omega*t);
46end
47figure
48box on
49hold on
50plot(t*Omega,real(rho11),'r-','Linewidth',2)
51plot(t*Omega,pe_exact,'b--','Linewidth',2)
52hold off
53xlabel('$\Omega t$','Interpreter','LaTex','Fontsize', 30)
54ylabel('$p_e(t)$','Interpreter','LaTex','Fontsize', 30)
55legend({'$\mbox{Numerical}$','$\mbox{Exact}$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
56set(gca,'fontsize',21)
57xlim([0 50])
59figure
60box on
61hold on
62plot(t*Omega,imag(sigmap),'r-','Linewidth',2)
63plot(t*Omega,imag(sp_exact),'b--','Linewidth',2)
64hold off
65xlabel('$\Omega t$','Interpreter','LaTex','Fontsize', 30)
66ylabel('$\mbox{Im}\langle \hat{\sigma}_+ \rangle$','Interpreter','LaTex','Fontsize', 30)
67legend({'$\mbox{Numerical}$','$\mbox{Exact}$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
68set(gca,'fontsize',21)
69xlim([0 50])