4.2.2.- S=1/2 con 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*)
Splus={{0,1},{0,0}}; (*Sigma minus*)
L=Splus;(*operador de Linblad de ganancias*)
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 ganancias para los cuatro componentes de un A genérico*)
RHSg11=Simplify[(I (H2.A-A.H2)-(g/2) ((Ldg .L .A)[[1,1]]+(A.Ldg .L)[[1,1]] -2 (Ldg. A.L)))[[1,1]] ,g∈Reals];
RHSg12=Simplify[(I (H2.A-A.H2)-(g/2) ((Ldg .L .A)[[1,2]]+(A.Ldg .L)[[1,2]] -2 (Ldg. A.L)))[[1,2]] ,g∈Reals];
RHSg21=Simplify[(I (H2.A-A.H2)-(g/2) ((Ldg .L .A)[[2,1]]+(A.Ldg .L)[[2,1]] -2 (Ldg. A.L)))[[2,1]] ,g∈Reals];
RHSg22=Simplify[(I (H2.A-A.H2)-(g/2) ((Ldg .L .A)[[2,2]]+(A.Ldg .L)[[2,2]] -2 (Ldg. A.L)))[[2,2]] ,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)=" ,RHSg11]
Print[Subscript["a'",12],"(t)=" ,RHSg12]
Print[Subscript["a'",21],"(t)=" ,RHSg21]
Print[Subscript["a'",22],"(t)=" ,RHSg22]

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

A2g1={{1,0},{0,1}} ; (*Evolución de la identidad 2x2 bajo operador de pérdidas, teniendo en cuenta su forma en t=0*)
A2g2={{0,(E^(-((g t)/2)+I t (Subscript[w, 1])-I t (Subscript[w, 2])))},{(E^(-((g t)/2)-I t (Subscript[w, 1])+I t (Subscript[w, 2]))),0}} ; (*Evolución de Sx*)
A2g3={{0,-I (E^(-((g t)/2)+I t (Subscript[w, 1])-I t (Subscript[w, 2])))},{I (E^(-((g t)/2)-I t (Subscript[w, 1])+I t (Subscript[w, 2]))),0}} ; (*Evolución de Sy*)
A2g4={{1,0},{0,1-2(E^(-g t))}} ; (*Evolución de Sz*)

Print[ Subscript["A",0],"(t)=" ,MatrixForm[A2g1]]
Print[ Subscript["A",1],"(t)=" ,MatrixForm[A2g2]]
Print[ Subscript["A",2],"(t)=" ,MatrixForm[A2g3]]
Print[ Subscript["A",3],"(t)=" ,MatrixForm[A2g4]]

A0vector={S0, Sx, Sy, Sz}; (*A_i(0)*)
Agtvector ={A2g1,A2g2, A2g3, A2g4}; (*A_i(t) bajo el operador de ganancias*)
Do[Do[Print["<",Subscript["A",n],"(0)",Subscript["A",m],"(t)", "> = ", Tr[Ro2.Agtvector[[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 g=0.2, con la proporcion "a" tal que C2g(t)=0 en algún punto*)

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

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

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