4.3.2.- S=1/2 con pérdidas y ganancias- 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*)
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 y ganancias para los cuatro componentes de un A genérico*)
RHSpg11=Simplify[(I (H2.A-A.H2)-(p/2) ((Splus .Sminus .A)[[1,1]]+(A.Splus .Sminus)[[1,1]] -2 (Splus. A.Sminus)))[[1,1]]-(g/2) ((Sminus .Splus .A)[[1,1]]+(A.Sminus .Splus)[[1,1]] -2 (Sminus. A.Splus))[[1,1]] ,{p∈Reals, g∈Reals}];
RHSpg12=Simplify[(I (H2.A-A.H2)-(p/2) ((Splus .Sminus .A)[[1,2]]+(A.Splus .Sminus)[[1,2]] -2 (Splus. A.Sminus)))[[1,2]]-(g/2) ((Sminus .Splus .A)[[1,2]]+(A.Sminus .Splus)[[1,2]] -2 (Sminus. A.Splus))[[1,2]] ,{p∈Reals, g∈Reals}];
RHSpg21=Simplify[(I (H2.A-A.H2)-(p/2) ((Splus .Sminus .A)[[2,1]]+(A.Splus .Sminus)[[2,1]] -2 (Splus. A.Sminus)))[[2,1]]-(g/2) ((Sminus .Splus .A)[[2,1]]+(A.Sminus .Splus)[[2,1]] -2 (Sminus. A.Splus))[[2,1]] ,{p∈Reals, g∈Reals}];
RHSpg22=Simplify[(I (H2.A-A.H2)-(p/2) ((Splus .Sminus .A)[[2,2]]+(A.Splus .Sminus)[[2,2]] -2 (Splus. A.Sminus)))[[2,2]]-(g/2) ((Sminus .Splus .A)[[2,2]]+(A.Sminus .Splus)[[2,2]] -2 (Sminus. A.Splus))[[2,2]] ,{p∈Reals, g∈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)=" ,RHSpg11]
Print[Subscript["a'",12],"(t)=" ,RHSpg12]
Print[Subscript["a'",21],"(t)=" ,RHSpg21]
Print[Subscript["a'",22],"(t)=" ,RHSpg22]

(*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](g+p-2I (Subscript[w, 1])+2I (Subscript[w, 2])), 
a21'[t]==-1/2a21[t] (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]

A2pg1={{1,0},{0,1}} ; (*Evolución de la identidad 2x2 bajo operadores de pérdidas y ganancias, teniendo en cuenta su forma en t=0*)
A2pg2={{0,(E^(-((g t)/2)-((p t)/2)+I t (Subscript[w, 1])-I t (Subscript[w, 2])))},{(E^(-((g t)/2)-((p t)/2)-I t (Subscript[w, 1])+I t (Subscript[w, 2]))),0}} ; (*Evolución de Sx*)
A2pg3={{0,-I (E^(-((g t)/2)-((p t)/2)+I t (Subscript[w, 1])-I t (Subscript[w, 2])))},{I (E^(-((g t)/2)-((p t)/2)-I t (Subscript[w, 1])+I t (Subscript[w, 2]))),0}} ; (*Evolución de Sy*)
A2pg4={{((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[A2pg1]]
Print[ Subscript["A",1],"(t)=" ,MatrixForm[A2pg2]]
Print[ Subscript["A",2],"(t)=" ,MatrixForm[A2pg3]]
Print[ Subscript["A",3],"(t)=" ,MatrixForm[A2pg4]]

A0vector={S0, Sx, Sy, Sz}; (*A_i(0)*)
Apgtvector ={A2pg1,A2pg2, A2pg3, A2pg4}; (*A_i(t) bajo el operador de pérdidas*)
Do[Do[Print["<",Subscript["A",n],"(0)",Subscript["A",m],"(t)", "> = ", Tr[Ro2.Apgtvector[[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, con la proporcion "a" tal que C2pg(t)=0 en algún punto*)

C2pg[t_]=((((((E^(-1))/((E^(-3))+(E^(-1))))+((E^(-3))/((E^(-3))+(E^(-1)))))+(((E^(-((0.2Pi)/(2 2))-((0.2Pi)/(2 2))-1))+(E^(-((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^(-((0.2t)/2)-((0.2t)/2)-3I t-1+I t))/((E^(-3))+(E^(-1))))+((E^(-((0.2t)/2)-((0.2t)/2)-3+3I t-I t ))/((E^(-3))+(E^(-1)))))+(((E^(-((0.2Pi)/(2 2))-((0.2Pi)/(2 2))-1))+(E^(-((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, con la proporción "a" igual al valor que tiene en C2pg[t_], cuando p=0.2 y g=0.2*)
C2pgEstable[t_]=((((((E^(-3))/(((E^(-3))+(E^(-1)))))+((E^(-1))/(((E^(-3))+(E^(-1))))))+(0.2+0.2)(E^(((-(0.2+0.2) ((Pi )/2))/2) ))((E^(-3))+(E^(-1)))(((-(E^(-1))(0.2-0.2-2 0.2 (E^(-(0.2+0.2) ((Pi )/2) )))+(E^(-3))(0.2-0.2+2 0.2 (E^(-(0.2+0.2) (Pi/2) )))))^(-1))))^(-1))((((E^(3I t-3-I t))/(((E^(-3))+(E^(-1)))))+((E^(I t-1-3I t))/(((E^(-3))+(E^(-1))))))+(0.2+0.2)(E^(((-(0.2+0.2) ((Pi )/2))/2) ))((E^(-3))+(E^(-1)))(((-(E^(-1))(0.2-0.2-2 0.2 (E^(-(0.2+0.2) ((Pi )/2) )))+(E^(-3))(0.2-0.2+2 0.2 (E^(-(0.2+0.2) (Pi/2) )))))^(-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=0*)
ParametricPlot[{{Re[C2pgEstable[t]],Im[C2pgEstable[t]]},{Re[C2pg[t]],Im[C2pg[t]]}},{t,0, 30Pi},AxesLabel->{Style["Re (OverTilde[C])(t)" ,Bold,16],Style["Im (OverTilde[C])(t)",Bold,16]}, FormatType->StandardForm]