c33ae98Companion to 'Coding closed and open quantum systems in MATLAB' (arXiv:1911.04906)Jeremy Magland 1Nspins = 6; % Number of spins
2H0 = zeros(2^Nspins,2^Nspins);
3H1 = zeros(2^Nspins,2^Nspins);
4J = 0;
5alpha = 0.2; % Parameter \alpha
6for j=1:Nspins
7 for i=j+1:Nspins
8 Jij = abs(i-j)^(-alpha);
9 J = J + (Nspins-1)^(-1)*Jij;
10 end
11end
12B = J/0.42; % Magnetic field
13Sx = [0 1;1 0];
14Sz = [1 0;0 -1];
15for i=1:Nspins
16 Szi = getSci(Sz,i,Nspins);
17 H0 = H0 - B*Szi; % Hamiltonian for the magnetic field
18 for j=1:Nspins
19 if i~=j
20 Sxi = getSci(Sx,i,Nspins);
21 Sxj = getSci(Sx,j,Nspins);
22 Vij = abs(i-j)^(-alpha)/J;
23 H1 = H1 - Vij*Sxi*Sxj; % Interaction Hamiltonian
24 end
25 end
26end
27H = H0 + H1; % Total hamiltonian
28xr = [1 1]'/sqrt(2); % Single-particle state |Psi_{-->}>
29xl = [-1 1]'/sqrt(2); % Single-particle state |Psi_{<--}>
30Xr = xr;
31Xl = xl;
32for n=1:Nspins-1
33 Xr = kron(Xr,xr); % Many-body state |Psi_{-->}>
34 Xl = kron(Xl,xl); % Many-body state |Psi_{<--}>
35end
36PSI_0 = Xr; % Initial condition
37ti = 0; % Initial time
38tf = 22; % Final time
39Nt = 10000; % Number of steps
40dt = (tf-ti)/(Nt-1); % Step time dt
41t = ti:dt:tf; % Time vector
42U = expm(-1i*H*dt); % Time propagator operator U(dt)
43Mx = zeros(size(t)); % Average Magnetization <M_x(t)>
44Lambda = zeros(size(t)); % Rate function \Lambda(t)
45SSx = 0;
46for i=1:Nspins
47 Sxi = getSci(Sx,i,Nspins);
48 SSx = SSx + Sxi/Nspins; % Magnetization operator
49end
50for n=1:length(t)
51 if n==1
52 PSI = PSI_0; % Initial wavefunction
53 else
54 PSI = U*PSI; % Wavefunction at time t_n
55 end
56 Pr = abs(Xr'*PSI)^2; % Probability state |Psi_{-->}>
57 Pl = abs(Xl'*PSI)^2; % Probability state |Psi_{<--}>
58 Lambda(n) = min(-Nspins^(-1)*log(Pr),-Nspins^(-1)*log(Pl));
59 Mx(n) = PSI'*SSx*PSI; % Average magnetization along x-axis
60end
62figure()
63plot(B*t,Lambda,'b-','Linewidth',3)
64xlabel('$B t$','Interpreter','LaTex','Fontsize', 30)
65ylabel('$\Lambda(t)$','Interpreter','LaTex','Fontsize', 30)
66set(gca,'fontsize',21)
67xlim([0 5])
69figure()
70box on
71plot(t*B,real(Mx),'r-','Linewidth',2)
72xlabel('$B t$','Interpreter','LaTex','Fontsize', 30)
73ylabel('$\langle M_x\rangle$','Interpreter','LaTex','Fontsize', 30)
74set(gca,'fontsize',21)
75xlim([0 100])