1sz = [1 0; % Pauli matrix s_z
2 0 -1];
3Nw = 5000; % Number of points for \omega
4wi = 0.01; % Initial frequency
5wf = 5; % Final frequency
6dw = (wf-wi)/(Nw-1); % Step frequency d\omega
7w = wi:dw:wf; % Frequency vector
8alpha = 5; % Strength spectral density function J1(w)
9s = 2.5; % Ohmic-parameter s
10wc = 0.1; % Cut-off frequency
11J1 = alpha*wc^(1-s)*w.^s.*exp(-w/wc); % Spectral density funcion J1(w)
12J0 = 0.2; % Strength spectral density function J2(w)
13w0 = 2; % Resonant frequency
14Gamma = 0.1; % Width of the spectral density function J2(w)
15J2 = J0*(Gamma/2)./((w-w0).^2+(Gamma/2)^2).*w.^s./(w/w0+1).^2; % Spectral density funcion J2(w)
16figure
17box on
18hold on
19plot(w,J1,'r-','Linewidth',2)
20plot(w,J2,'b-','Linewidth',2)
21hold off
22xlabel('$\omega$','Interpreter','LaTex','Fontsize', 30)
23ylabel('$J(\omega)$','Interpreter','LaTex','Fontsize', 30)
24legend({'$J_1(\omega)$','$J_2(\omega)$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
25set(gca,'fontsize',21)
26xlim([0 5])
27Nt = 5000; % Number of points for time
28ti = 0; % Initial time
29tf = 100; % Final time
30dt = (tf-ti)/(Nt-1); % Step time dt
31t = ti:dt:tf; t= t'; % Time vector
32Psi_0 = [1 1]'/sqrt(2); % Initial wavefunction
33rho1 = Psi_0*Psi_0'; % Initial density matrix
34rho2 = rho1;
35p11 = zeros(size(t));
36p22 = zeros(size(t));
37p12 = zeros(size(t));
38p21 = zeros(size(t));
39wa = ones(size(t))*w; % Auxiliar frecuency vector
40J1a = ones(size(t))*J1; % Auxiliar J1 vector
41J2a = ones(size(t))*J2; % Auxiliar J2 vector
42ta = t*ones(size(w)); % Auxiliar time vector
43T = 0.001*w0; % Temperature
44gamma1 = sum(J1a./wa.*sin(wa.*ta).*coth(wa/T/2),2)*dw; % Rate gamma_1(t)
45gamma_1_teo = alpha*wc*gamma(s)*sin(s*atan(wc*t))./(1+(wc*t).^2).^(s/2);
46gamma2 = sum(J2a./wa.*sin(wa.*ta).*coth(wa/T/2),2)*dw; % Rate gamma_2(t)
47figure()
48box on
49hold on
50plot(t,gamma1,'r-','Linewidth',2)
51plot(t,gamma_1_teo,'k--','Linewidth',2)
52plot(t,gamma2,'b-','Linewidth',2)
53hold off
54xlabel('$t$','Interpreter','LaTex','Fontsize', 30)
55ylabel('$\gamma(t)$','Interpreter','LaTex','Fontsize', 30)
56legend({'$\gamma_1(t)$','$\gamma_1^{\rm teo}(t)$','$\gamma_2(t)$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
57set(gca,'fontsize',21)
58xlim([0 100])
59nk = 15; % nk steps per interval dt
60C1 = zeros(size(t));
61C2 = zeros(size(t));
62N1 = zeros(size(t));
63N2 = zeros(size(t));
64for n = 1:length(t)
65 for k = 1:nk
66 L1 = gamma1(n)*(sz*rho1*sz -rho1); % Lindbladian for gamma_1
67 L2 = gamma2(n)*(sz*rho2*sz -rho2); % Lindbladian for gamma_2
68 rho1_pred = rho1 + L1*dt/nk; % Predictor \rho_1
69 rho2_pred = rho2 + L2*dt/nk; % Predictor \rho_2
70 rho1_m = 0.5*(rho1+rho1_pred); % Mean \rho_1
71 rho2_m = 0.5*(rho1+rho2_pred); % Mean \rho_1
72 L1 = gamma1(n)*(sz*rho1_m*sz -rho1_m); % Lindbladian using mean \rho_1
73 L2 = gamma2(n)*(sz*rho2_m*sz -rho2_m); % Lindbladian using mean \rho_2
74 rho1 = rho1 + L1*dt/nk; % Density matrix \rho_1
75 rho2 = rho2 + L2*dt/nk; % Density matrix \rho_2
76 end
77 C1(n) = 2*abs(rho1(1,2)); % Coherence C1(t)
78 C2(n) = 2*abs(rho2(1,2)); % Coherence C2(t)
79 f1 = (abs(gamma1(1:n))-gamma1(1:n));
80 N1(n) = sum(f1)*dt; % Degree of non-Markovianity for gamma_1
81 f2 = (abs(gamma2(1:n))-gamma2(1:n));
82 N2(n) = sum(f2)*dt; % Degree of non-Markovianity for gamma_2
83end
84figure()
85box on
86hold on
87plot(t,C1,'r-','Linewidth',2)
88plot(t,C2,'b-','Linewidth',2)
89hold off
90xlabel('$t$','Interpreter','LaTex','Fontsize', 30)
91ylabel('$C(t)$','Interpreter','LaTex','Fontsize', 30)
92legend({'$C_1(t)$','$C_2(\omega)$'},'Interpreter','latex','Fontsize', 21,'Location','northeast')
93set(gca,'fontsize',21)
94xlim([0 50])
96figure()
97box on
98hold on
99plot(t,N1,'r-','Linewidth',2)
100plot(t,N2,'b-','Linewidth',2)
101hold off
102xlabel('$t$','Interpreter','LaTex','Fontsize', 30)
103ylabel('$\mathcal{N}_{\gamma}(t)$','Interpreter','LaTex','Fontsize', 30)
104legend({'$\mathcal{N}_{\gamma_1}(t)$','$\mathcal{N}_{\gamma_2}(t)$'},'Interpreter','latex','Fontsize', 21,'Location','best')
105set(gca,'fontsize',21)
106xlim([0 100])