4.5.2.- S=1/2 con pérdidas, ganancias y 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. operador de Linblad de defasaje. SzDagger*)
Sminus={{0,0},{1,0}}; (*operador de Linblad de pérdidas. SplusDagger*)
Splus={{0,1},{0,0}};(*operador de Linblad de ganancias. SminusDagger*)
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 pérdidas ganancias y defasaje para los cuatro componentes de un A genérico*)
RHSpgd11=Simplify[(
I (H2.A-A.H2)
-(p/2) ((Splus .Sminus .A)+(A.Splus .Sminus)-2 (Splus. A.Sminus))
-(g/2) ((Sminus .Splus .A)+(A.Sminus .Splus) -2 (Sminus. A.Splus))
-(d/2) ((Sz .Sz .A)+(A.Sz .Sz) -2 (Sz. A.Sz))
)[[1,1]] ,{p∈Reals, g∈Reals, d∈Reals}];
RHSpgd12=Simplify[(
I (H2.A-A.H2)
-(p/2) ((Splus .Sminus .A)+(A.Splus .Sminus)-2 (Splus. A.Sminus))
-(g/2) ((Sminus .Splus .A)+(A.Sminus .Splus) -2 (Sminus. A.Splus))
-(d/2) ((Sz .Sz .A)+(A.Sz .Sz) -2 (Sz. A.Sz))
)[[1,2]] ,{p∈Reals, g∈Reals, d∈Reals}];
RHSpgd21=Simplify[(
I (H2.A-A.H2)
-(p/2) ((Splus .Sminus .A)+(A.Splus .Sminus)-2 (Splus. A.Sminus))
-(g/2) ((Sminus .Splus .A)+(A.Sminus .Splus) -2 (Sminus. A.Splus))
-(d/2) ((Sz .Sz .A)+(A.Sz .Sz) -2 (Sz. A.Sz))
)[[2,1]] ,{p∈Reals, g∈Reals, d∈Reals}];
RHSpgd22=Simplify[(
I (H2.A-A.H2)
-(p/2) ((Splus .Sminus .A)+(A.Splus .Sminus)-2 (Splus. A.Sminus))
-(g/2) ((Sminus .Splus .A)+(A.Sminus .Splus) -2 (Sminus. A.Splus))
-(d/2) ((Sz .Sz .A)+(A.Sz .Sz) -2 (Sz. A.Sz))
)[[2,2]] ,{p∈Reals, g∈Reals, 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)=" ,RHSpgd11]
Print[Subscript["a'",12],"(t)=" ,RHSpgd12]
Print[Subscript["a'",21],"(t)=" ,RHSpgd21]
Print[Subscript["a'",22],"(t)=" ,RHSpgd22]

(*Resolución de las acuaciones diferenciales acopladas obtenidas en el paso anterior*)
DSolve[{a11'[t]==p(-a11[t]+a22[t]), 
a12'[t]==-1/2a12[t](4d+g+p-2I (Subscript[w, 1])+2I (Subscript[w, 2])), 
a21'[t]==-1/2a21[t] (4d+g+p+2 I (Subscript[w, 1])-2 I (Subscript[w, 2])),
a22'[t]==g(a11[t]-a22[t])},
 {a11[t],a12[t],a21[t],a22[t]} ,t]

A2pgd1={{1,0},{0,1}} ; (*Evolución de la identidad 2x2 bajo operadores de pérdidas ganancias y defasaje, teniendo en cuenta su forma en t=0*)
A2pgd2={{0,(E^(-2d t-((g t)/2)-((p t)/2)+I t (Subscript[w, 1])-I t (Subscript[w, 2])))},{(E^(-2d t-((g t)/2)-((p t)/2)-I t (Subscript[w, 1])+I t (Subscript[w, 2]))),0}} ; (*Evolución de Sx*)
A2pgd3={{0,-I (E^(-2d t-((g t)/2)-((p t)/2)+I t (Subscript[w, 1])-I t (Subscript[w, 2])))},{I (E^(-2d t-((g t)/2)-((p t)/2)-I t (Subscript[w, 1])+I t (Subscript[w, 2]))),0}} ; (*Evolución de Sy*)
A2pgd4={{((g-p+2(E^((-g-p) t)) p)/(g+p)),0},{0,((g-p -2(E^((-g-p) t)) g)/(g+p))}} ; (*Evolución de Sz*)

Print[ Subscript["A",0],"(t)=" ,MatrixForm[A2pgd1]]
Print[ Subscript["A",1],"(t)=" ,MatrixForm[A2pgd2]]
Print[ Subscript["A",2],"(t)=" ,MatrixForm[A2pgd3]]
Print[ Subscript["A",3],"(t)=" ,MatrixForm[A2pgd4]]

A0vector={S0, Sx, Sy, Sz}; (*A_i(0)*)
Apgdtvector ={A2pgd1,A2pgd2, A2pgd3, A2pgd4}; (*A_i(t) bajo el operador de pérdidas ganancias y defasaje*)
Do[Do[Print["<",Subscript["A",n],"(0)",Subscript["A",m],"(t)", "> = ", Tr[Ro2.Apgdtvector[[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, p=0.2, g=0.2, d=0.2, con la proporcion "a" tal que C2pg(t)=0 en algún punto*)

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

(*Correlación A1+a.A3 para w1=3, w2=1, b=1, p=0, g=0, d=0, con la proporción "a" igual al valor que tiene en C2pgd[t_], cuando p=0.2 g=0.2 y d=0.2*)
C2pgtEstable[t_]=((((1)+(((E^(-((2 0.2Pi)/2)-((0.2Pi)/(2 2))-((0.2Pi)/(2 2))-1))+(E^(-((2 0.2Pi)/(2 ))-((0.2Pi)/(2 2))-((0.2Pi)/(2 2))-3)))/((1/((0.2+0.2) ))((E^(-1)) (0.2-2 (E^((-0.2-0.2) (Pi/2))) 0.2-0.2)+(E^(-3))(0.2-0.2+2 (E^((-0.2-0.2) (Pi/2))) 0.2))))(((E^(-3))-(E^(-1)) )/( ((E^(-3))+(E^(-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)-((0.2Pi)/(2 2))-((0.2Pi)/(2 2))-1))+(E^(-((2 0.2Pi)/(2 ))-((0.2Pi)/(2 2))-((0.2Pi)/(2 2))-3)))/((1/((0.2+0.2) ))((E^(-1)) (0.2-2 (E^((-0.2-0.2) (Pi/2))) 0.2-0.2)+(E^(-3))(0.2-0.2+2 (E^((-0.2-0.2) (Pi/2))) 0.2))))(((E^(-3))-(E^(-1)) )/( ((E^(-3))+(E^(-1))))))

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