4.4.2- S=1/2 con defasaje - Heisenberg

S0={{1,0},{0,1}} ; (*Identidad 2x2*)
Sx={{0,1},{1,0}} ; (*Sigma x*)
Sy={{0,-I},{I,0}} ; (*Sigma y*)
Sz={{1,0},{0,-1}} ; (*Sigma z*)
L=Sz;(*operador de Linblad de defasaje*)
Ldg=ConjugateTranspose[L]; (*Ldagger*)
A={{(Subscript[a, 11]),(Subscript[a, 12])},{(Subscript[a, 21]),(Subscript[a, 22])}}; 
H2= {{(Subscript[w, 1]), 0},{0,  (Subscript[w, 2])}}; (*Hamiltoniano 2-dimensional*)
Ro2 = 1/Z {{Exp[-b (Subscript[w, 1])], 0},{0, Exp[-b (Subscript[w, 2])]}}; (*Matriz de densidad 2x2, constante en esta imagen*)


(*Calcular el Right Hand Side de la ecuación de Lindblad con defasaje para los cuatro componentes de un A genérico*)
RHSd11=Simplify[(I (H2.A-A.H2)-(d/2) ((Ldg .L .A)[[1,1]]+(A.Ldg .L)[[1,1]] -2 (Ldg. A.L)))[[1,1]] ,d∈Reals];
RHSd12=Simplify[(I (H2.A-A.H2)-(d/2) ((Ldg .L .A)[[1,2]]+(A.Ldg .L)[[1,2]] -2 (Ldg. A.L)))[[1,2]] ,d∈Reals];
RHSd21=Simplify[(I (H2.A-A.H2)-(d/2) ((Ldg .L .A)[[2,1]]+(A.Ldg .L)[[2,1]] -2 (Ldg. A.L)))[[2,1]] ,d∈Reals];
RHSd22=Simplify[(I (H2.A-A.H2)-(d/2) ((Ldg .L .A)[[2,2]]+(A.Ldg .L)[[2,2]] -2 (Ldg. A.L)))[[2,2]] ,d∈Reals];

(*Left Hand Side = Right Hand Side. Las ecuaciones que surgen para cada componente de A de la ecuación de Lindblad*)
Print[Subscript["a'",11],"(t)=" ,RHSd11]
Print[Subscript["a'",12],"(t)=" ,RHSd12]
Print[Subscript["a'",21],"(t)=" ,RHSd21]
Print[Subscript["a'",22],"(t)=" ,RHSd22]

(*Resolución de las acuaciones diferenciales desacopladas obtenidas en el paso anterior*)
DSolve[a11'[t]==0, a11[t], t]
DSolve[a12'[t]==-a12[t](2d-I (Subscript[w, 1])+I (Subscript[w, 2])), a12[t], t]
DSolve[a21'[t]==-a21[t](2d+I (Subscript[w, 1])-I (Subscript[w, 2])), a21[t], t]
DSolve[a22'[t]==0, a22[t], t]

A2d1={{1,0},{0,1}} ; (*Evolución de la identidad 2x2 bajo operador de pérdidas, teniendo en cuenta su forma en t=0*)
A2d2={{0,(E^(-2d t+I t (Subscript[w, 1])-I t (Subscript[w, 2])))},{(E^(-2d t-I t (Subscript[w, 1])+I t (Subscript[w, 2]))),0}} ; (*Evolución de Sx*)
A2d3={{0,-I (E^(-2d t+I t (Subscript[w, 1])-I t (Subscript[w, 2])))},{I (E^(-2d t-I t (Subscript[w, 1])+I t (Subscript[w, 2]))),0}} ; (*Evolución de Sy*)
A2d4={{1,0},{0,-1}} ; (*Evolución de Sz*)

Print[ Subscript["A",0],"(t)=" ,MatrixForm[A2d1]]
Print[ Subscript["A",1],"(t)=" ,MatrixForm[A2d2]]
Print[ Subscript["A",2],"(t)=" ,MatrixForm[A2d3]]
Print[ Subscript["A",3],"(t)=" ,MatrixForm[A2d4]]

A0vector={S0, Sx, Sy, Sz}; (*A_i(0)*)
Adtvector ={A2d1,A2d2, A2d3, A2d4}; (*A_i(t) bajo el operador de defasaje*)
Do[Do[Print["<",Subscript["A",n],"(0)",Subscript["A",m],"(t)", "> = ", Tr[Ro2.Adtvector[[n+1]].A0vector[[m+1]]]],{n, 0,3}],{m,0,3}] (* los índices del vector en mathematica van del 1 al 4, de ahí el n+1 y el m+1 *)

(*Valores de las autocorrelaciones C_{nm}*)

(*Correlación normalizada A1+a.A3 para w1=3, w2=1, b=1 y d=0.2, con la proporcion "a" tal que C2d(t)=0 en algún punto*)

C2d[t_]=(((1+(E^(-((2 0.2Pi)/2) ))(1)))^(-1))((((E^(- 2 0.2t-3I t-1+I t))/((E^(-3))+(E^(-1))))+((E^(-2 0.2t-3+3I t-I t ))/((E^(-3))+(E^(-1)))))+(E^(-((2 0.2Pi)/2) ))(1))

(*Correlación A1+a.A3 para w1=3, w2=1, b=1 y d=0, con la proporción "a" igual al valor que tiene en C2d[t_], cuando d=0.2*)
C2dEstable[t_]=(((1+(E^(-((2 0.2Pi)/2) ))(1)))^(-1))((((E^(-3I t-1+I t))/((E^(-3))+(E^(-1))))+((E^(-3+3I t-I t ))/((E^(-3))+(E^(-1)))))+(E^(-((2 0.2Pi)/2) ))(1))

(*En amarillo, gráfica con los valores ajustados para que C2d(t)=0 en algún punto. En azul la órbita que conseguiríamos si, manteniendo el resto de valores inalterados, hiciesemos d=0*)
ParametricPlot[{{Re[C2dEstable[t]],Im[C2dEstable[t]]},{Re[C2d[t]],Im[C2d[t]]}},{t,0, 30Pi},AxesLabel->{Style["Re (OverTilde[C])(t)" ,Bold,16],Style["Im (OverTilde[C])(t)",Bold,16]}, FormatType->StandardForm]