1dim = 4; % Dimension of the total Hilbert space
2Is = eye(dim); % Identity matrix for the total Hilbert space
3J = 1; % Value of the coupling term J
4B = 0.1*J; % Value of the magnetic field B
5Sx = [0 1;1 0]; % S_x operator for one spin
6Sz = [1 0;0 -1]; % S_z operator for one spin
7I = eye(2); % Identity matrix for one spin 1/2
8H = -J*kron(Sx,Sx)-B*(kron(Sx,I)+kron(I,Sx)); % Hamiltonian of the system
9L_H = -1i*kron(Is,H)+1i*kron(H.',Is); % Lindblad operator L_H
10gamma_1 = 0.1*B; % Decay rate gamma_1
11gamma_2 = 0.5*B; % Decay rate gamma_2
12S_minus =[0 0; 1 0]; % Lowering operator of the particle 1
13L1 = kron(S_minus,I); % Lowering operator of the particle 1 in the total Hilbert space
14DL_1 = gamma_1*(kron(conj(L1),L1)-0.5*kron(Is,L1'*L1)-0.5*kron(L1.'*conj(L1),Is));
15L2 = kron(I,S_minus); % Lowering operator of the particle 2 in the total Hilbert space
16DL_2 = gamma_2*(kron(conj(L2),L2)-0.5*kron(Is,L2'*L2)-0.5*kron(L2.'*conj(L2),Is));
17L_diss = DL_1 + DL_2; % Lindbladian L_diss
18L = L_H + L_diss; % Total Lindblad operator
19TOL = 1e-10;
20[R_sort,L_sort,lambda_sort] = sortingEigenvalues(dim,TOL,L);
22down = [0 1]'; % Quantum state down = [0 1]
23Psi_0 = kron(down,down); % Initial wavefunction
24rho_0 = Psi_0*Psi_0'; % Initial density matrix
25Nt = 3000; % Number of steps for time
26T = 2*pi/(2*B) ; % Period of time
27ti = 0; % Intial time
28tf = 4*T; % Final time
29dt = (tf-ti)/(Nt-1); % Step for time
30t = ti:dt:tf; % Time vector
31Mz = zeros(size(t)); % Initial average magnetization
32SSz = (kron(Sz,I)+kron(I,Sz))/2; % Operator S1^z+S2^z
33for n=1:length(t) % Iteration to find general sulution of rho(t)
34 rho = zeros(dim,dim);
35 for k=1:length(lambda_sort)
36 Lk = L_sort{k};
37 Rk = R_sort{k};
38 ck = trace(rho_0*Lk);
39 rho = rho + ck*exp(lambda_sort(k)*t(n))*Rk; % General sulution for rho(t)
40 end
41 Mz(n) = trace(SSz*rho); % General sulution for Mz(t)
42end
43figure()
44hold on
45plot(t/T,real(Mz),'r-','LineWidth',3)
46plot(t/T, exp(real(lambda_sort(end))*t),'k--','LineWidth',2)
47plot(t/T,-exp(real(lambda_sort(end))*t),'k--','LineWidth',2)
48hold off
49xlabel('$Bt$','Interpreter','LaTex','Fontsize', 30)
50ylabel('$\langle M_z \rangle$','Interpreter','LaTex','Fontsize', 30)
51set(gca,'fontsize',21)