(* Write down the Hamiltonians he, hs, v. We will use Fs (= E in the paper) as the energy of the system, and Fe = (E_0 in the paper) as the energy of the environment. *)
a = {{0,1},{0,0}};
astar = ConjugateTranspose[a];
hs = Fs(astar.a);

b={{0,1},{0,0}};
bstar = ConjugateTranspose[b];
he = Fe(bstar.b);
h0 = KroneckerProduct[he,IdentityMatrix[dimhs]] + KroneckerProduct[IdentityMatrix[dimhe],hs] ;

(* We define both interaction potentials, with the rotating wave approximation and without. Note that Mathematica's convention for taking tensor products is the opposite of ours, so we write the arguments in the opposite order. That is, if we want A tensor B, we write KroneckerProduct[B,A]. *)
vfull=( λ/2)* KroneckerProduct[b + ConjugateTranspose[b],a+astar] ;
vrot= ( λ/2)*( KroneckerProduct[b,astar] + KroneckerProduct[ConjugateTranspose[b],a]);
(* Now we select which potential: *)
v=vfull;

(* now put our definitions together to get h *)
dimhs=First[Dimensions[hs]];
dimhe=First[Dimensions[he]];
h = h0+v;
(* We'll display the matrix of h to check. *)
MatrixForm[h]

(* compute time evolution operator. It seems to run faster by computing each U of and Ustar, instead of setting Ustar = ConjugateTranspose[U]. *)
U =FullSimplify[ MatrixExp[I τ h],{Fe >0,Fs>0,β>0,τ>0,λ>0}];
Ustar =FullSimplify[ MatrixExp[-I τ h],{Fe >0,Fs>0,β>0,τ>0,λ>0}];

(* the state of the environment *)
rhobetaenv = MatrixExp[-β he]/Tr[MatrixExp[-β he]];

(* Partial trace from Prof. Cubitt: http://www.dr-qubit.org/matlab.php *)
PT[p_,a_,b_,sys_]:=Module[{s},s=Partition[p,{b,b}];
If[sys==2,s=Transpose[s,{3,4,1,2}]];
s=Tr[s,Plus,2]]

(* Now we wish to compute L. We'll initialize a variable LMatrixT, written with a T because our computation will give us the transpose of L. *)
LMatrixT=.;

(* First we write a basis for the Hilbert space B(H_S), bounded operators on the system's Hilbert space, i.e. dim h_s by dim h_s matrices: *)
For[i=1, i<dimhs+1,i++,
For[j=1,j<dimhs+1,j++,
A=Normal[SparseArray[{i,j}->1,dimhs]];

(* Now we act our operator on this basis element  A *)
LofA =PT[U.KroneckerProduct[rhobetaenv,A].Ustar,dimhe,dimhs,1];

(* This yields a dim h_s by dim h_s matrix. We want to write down L as a (dim h_s * dim h_e) x ( dim h_s * dim h_e) matrix, so we'll flatten this L(A) into a column vector, and simplify it: *)
LofAvector =FullSimplify[Flatten[LofA,1],{Fe >0,Fs>0,β>0,τ>0,λ>0}];


(* Okay, now we'll append this vector onto LMatrixT. If this is the first basis element, we'll set LMatrixT to be this result. Otherwise, we'll append it as the next row. Since we append as rows, we need to take the transpose at the end so L(A_11) is the first column of the matrix, etc. *)
LMatrixT = If[i==1&&j==1,{LofAvector},
Append[LMatrixT,LofAvector]];

]
]
LMatrix=Transpose[LMatrixT];
(* Now we have our operator L written as a matrix in the basis defined above. We can print it: *)
Print[MatrixForm[LMatrix]];

(* Now we can find the eigenvalues and eigenvectors of LMatrix. The eigenvalues will allow us to check our spectral assumptions, and the eigenvectors will give us the invariant state. This command takes a long time to run! It ran overnight on my desktop before finishing. So I'll comment it out and put in the results I get.*)

(* eigsys = FullSimplify[Eigensystem[LMatrix],{Fe >0,Fs>0,β>0,τ>0,λ>0}]; *)

eigsys={{(1/(8 ((((Fe-Fs))^2)+(λ^2)) ((((Fe+Fs))^2)+(λ^2))))(E^(-I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) ((((Fe+Fs))^2) (λ^2)+(λ^4)+(E^(2 I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (λ^2) ((((Fe+Fs))^2)+(λ^2))-2 (E^((1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (Sqrt[(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))]) (λ^2) Abs[-2 ((Fe^2)+(Fs^2)+(λ^2))+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]+((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]]+2 (E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (4 ((((Fe^2)-(Fs^2)))^2)+6 ((Fe^2)+(Fs^2)) (λ^2)+2 (λ^4)+(λ^2) ((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])),(1/(8 ((((Fe-Fs))^2)+(λ^2)) ((((Fe+Fs))^2)+(λ^2))))(E^(-I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) ((((Fe+Fs))^2) (λ^2)+(λ^4)+(E^(2 I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (λ^2) ((((Fe+Fs))^2)+(λ^2))+2 (E^((1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (Sqrt[(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))]) (λ^2) Abs[-2 ((Fe^2)+(Fs^2)+(λ^2))+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]+((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]]+2 (E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (4 ((((Fe^2)-(Fs^2)))^2)+6 ((Fe^2)+(Fs^2)) (λ^2)+2 (λ^4)+(λ^2) ((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])),(1/(4 (Sqrt[((((Fe-Fs))^2)+(λ^2)) ((((Fe+Fs))^2)+(λ^2))])))(E^(-(1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (2 (1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (Sqrt[(Fe^4)+2 (Fe^2) (-(Fs^2)+(λ^2))+((((Fs^2)+(λ^2)))^2)]) Cos[(1/2) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]-2 I (-1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (Fe-Fs) (Fe+Fs) Sin[(1/2) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]-2 √((E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (-2 ((((Fe^2)-(Fs^2)))^2)-2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)-(λ^2) (-4 Fe Fs+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]+Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ] (-(λ^2) (4 Fe Fs+(λ^2))+(2 ((((Fe^2)-(Fs^2)))^2)+2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])+2 (Fe-Fs) (Fe+Fs) (Sqrt[((((Fe^2)-(Fs^2)))^2)+2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)]) Sin[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ] Sin[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]))),(1/(4 (Sqrt[((((Fe-Fs))^2)+(λ^2)) ((((Fe+Fs))^2)+(λ^2))])))(E^(-(1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (2 (1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (Sqrt[(Fe^4)+2 (Fe^2) (-(Fs^2)+(λ^2))+((((Fs^2)+(λ^2)))^2)]) Cos[(1/2) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]+2 (-I (-1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (Fe-Fs) (Fe+Fs) Sin[(1/2) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]+√((E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (-2 ((((Fe^2)-(Fs^2)))^2)-2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)-(λ^2) (-4 Fe Fs+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]+Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ] (-(λ^2) (4 Fe Fs+(λ^2))+(2 ((((Fe^2)-(Fs^2)))^2)+2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])+2 (Fe-Fs) (Fe+Fs) (Sqrt[((((Fe^2)-(Fs^2)))^2)+2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)]) Sin[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ] Sin[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]))))},{{((E^(-(1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) ((1+(E^(Fe β))) (λ^2) (Sqrt[(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) ((((Fe-Fs))^2)+(λ^2)) ((((Fe+Fs))^2)+(λ^2))]) Abs[-2 ((Fe^2)+(Fs^2)+(λ^2))+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]+((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]]+(E^((1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (-1+(E^(Fe β))) (λ^2) (Sqrt[(Fe^4)+2 (Fe^2) (-(Fs^2)+(λ^2))+((((Fs^2)+(λ^2)))^2)]) (-4 Fe Fs+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]-((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])))/(2 (λ^2) (Sqrt[(Fe^4)+2 (Fe^2) (-(Fs^2)+(λ^2))+((((Fs^2)+(λ^2)))^2)]) (((((Fe+Fs))^2)+(λ^2)) (-1+Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ])+(E^(Fe β)) ((((Fe-Fs))^2)+(λ^2)) (-1+Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]))),0,0,1},{((E^(-(1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (-(1+(E^(Fe β))) (λ^2) (Sqrt[(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) ((((Fe-Fs))^2)+(λ^2)) ((((Fe+Fs))^2)+(λ^2))]) Abs[-2 ((Fe^2)+(Fs^2)+(λ^2))+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]+((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]]+(E^((1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (-1+(E^(Fe β))) (λ^2) (Sqrt[(Fe^4)+2 (Fe^2) (-(Fs^2)+(λ^2))+((((Fs^2)+(λ^2)))^2)]) (-4 Fe Fs+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]-((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])))/(2 (λ^2) (Sqrt[(Fe^4)+2 (Fe^2) (-(Fs^2)+(λ^2))+((((Fs^2)+(λ^2)))^2)]) (((((Fe+Fs))^2)+(λ^2)) (-1+Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ])+(E^(Fe β)) ((((Fe-Fs))^2)+(λ^2)) (-1+Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]))),0,0,1},{0,(1/((-1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (λ^2)))((1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (Fe+Fs) (Sqrt[(((Fe-Fs))^2)+(λ^2)])+I (-1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (Fe-Fs) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) Cot[(1/2) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]-I Csc[(1/2) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ] √((E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (-2 ((((Fe^2)-(Fs^2)))^2)-2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)-(λ^2) (-4 Fe Fs+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]+Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ] (-(λ^2) (4 Fe Fs+(λ^2))+(2 ((((Fe^2)-(Fs^2)))^2)+2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])+2 (Fe-Fs) (Fe+Fs) (Sqrt[((((Fe^2)-(Fs^2)))^2)+2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)]) Sin[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ] Sin[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]))),1,0},{0,(1/((-1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (λ^2)))((1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (Fe+Fs) (Sqrt[(((Fe-Fs))^2)+(λ^2)])+I (-1+(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))) (Fe-Fs) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) Cot[(1/2) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]+I Csc[(1/2) (Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ] √((E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (-2 ((((Fe^2)-(Fs^2)))^2)-2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)-(λ^2) (-4 Fe Fs+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]+Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ] (-(λ^2) (4 Fe Fs+(λ^2))+(2 ((((Fe^2)-(Fs^2)))^2)+2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])+2 (Fe-Fs) (Fe+Fs) (Sqrt[((((Fe^2)-(Fs^2)))^2)+2 ((Fe^2)+(Fs^2)) (λ^2)+(λ^4)]) Sin[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ] Sin[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]))),1,0}}};

(* We remark the eigenvalues are independent of β. We can see this by asking
Print[eigsys[[1,j]] for each j=1,2,3,4, and finding no β's.*)

(* From here, we need to identify the invariant state as the eigenvector associated to the constant eigenvalue 1, and then normalize by the trace. By plotting each eigenvalue (substituting in numbers for each parameter, we notice the eigenvalue eigsys[[1,2]] is a candidate for this constant eigenvalue. We expect that it we may have problems (i.e. it is not constantly 1) in some exceptional set. Let us investigate this eigenvalue. *)

Print[eigsys[[1,2]]]

(* We think this should be 1, so we see what stops Mathematica from simplifying further. We see we have an absolute value. Let's consider the argument of the absolute value:*)
(* A = -2 ((Fe^2)+(Fs^2)+(λ^2))+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]+((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]
We estimate Cosine above by 1:
A≤ -2 ((Fe^2)+(Fs^2)+(λ^2))+((((Fe+Fs))^2)+(λ^2))+((((Fe-Fs))^2)+(λ^2)) =-2 Fe^2 -2 Fs^2 + Fe^2  + Fs^2 + 2Fe*Fs + Fe^2 + Fs^2 - 2Fe*Fs = 0 
Okay, so fairly easy to see this term is always negative. Not sure why Mathematica couldn't spot that. So let's plug in the negative version in:*)

eigsys[[1,2]]=(1/(8 ((((Fe-Fs))^2)+(λ^2)) ((((Fe+Fs))^2)+(λ^2))))(E^(-I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) ((((Fe+Fs))^2) (λ^2)+(λ^4)+(E^(2 I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (λ^2) ((((Fe+Fs))^2)+(λ^2))+2 (E^((1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (Sqrt[(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))]) (λ^2)*(-1)* (-2 ((Fe^2)+(Fs^2)+(λ^2))+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]+((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])+2 (E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (4 ((((Fe^2)-(Fs^2)))^2)+6 ((Fe^2)+(Fs^2)) (λ^2)+2 (λ^4)+(λ^2) ((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]));

FullSimplify[eigsys[[1,2]],{Fe >0,Fs>0,β>0,τ>0,λ>0}]

(* Huh, still not there yet. We notice Mathematica wrote  (Sqrt[(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))]), and if we took the principal root, (Sqrt[(E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))]) = (E^(I (1/2)(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)), we get simpliciation: *)

Simplify[(1/(8 ((((Fe-Fs))^2)+(λ^2)) ((((Fe+Fs))^2)+(λ^2))))(E^(-I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) ((((Fe+Fs))^2) (λ^2)+(λ^4)+(E^(2 I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (λ^2) ((((Fe+Fs))^2)+(λ^2))-2  (E^((1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (E^((1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ))  (λ^2) (-2 ((Fe^2)+(Fs^2)+(λ^2))+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]+((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])+2 (E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (4 ((((Fe^2)-(Fs^2)))^2)+6 ((Fe^2)+(Fs^2)) (λ^2)+2 (λ^4)+(λ^2) ((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]))]

(* Great, constantly 1! *)
eigsys[[1,2]] = 1;

(* Now let's consider the associated eigenvector:*)
eigsys[[2,2]]

(* We make the same simplifications as for the eigenvalue: *)
eigsys[[2,2]] = FullSimplify[{((E^(-(1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) ((-1-(E^(Fe β))) (λ^2) (E^(I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ/2)) Sqrt[((((Fe-Fs))^2)+(λ^2)) ((((Fe+Fs))^2)+(λ^2))] (-1)*(-2 ((Fe^2)+(Fs^2)+(λ^2))+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]+((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])+(E^((1/2) I (Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ)) (-1+(E^(Fe β))) (λ^2) (Sqrt[(Fe^4)+2 (Fe^2) (-(Fs^2)+(λ^2))+((((Fs^2)+(λ^2)))^2)]) (-4 Fe Fs+((((Fe+Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ]-((((Fe-Fs))^2)+(λ^2)) Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ])))/(2 (λ^2) (Sqrt[(Fe^4)+2 (Fe^2) (-(Fs^2)+(λ^2))+((((Fs^2)+(λ^2)))^2)]) (((((Fe+Fs))^2)+(λ^2)) (-1+Cos[(Sqrt[(((Fe-Fs))^2)+(λ^2)]) τ])+(E^(Fe β)) ((((Fe-Fs))^2)+(λ^2)) (-1+Cos[(Sqrt[(((Fe+Fs))^2)+(λ^2)]) τ]))),0,0,1},{Fe >0,Fs>0,β>0,τ>0,λ>0}]

(* And then we normalize to be trace 1: *)
rhoinv =FullSimplify[ eigsys[[2,2]] / Tr[eigsys[[2,2]]],{Fe >0,Fs>0,β>0,τ>0,λ>0}]

(* We can check this is the same as what I wrote in the paper: *)
rhoinvpaper=With[{ν=Sqrt[(Fe-Fs)^2 + λ^2], η=Sqrt[(Fe+Fs)^2 + λ^2]},{{(ν^2 (1-Cos[η τ])+E^(Fe β) η^2 (1-Cos[ν τ]))/((1+E^(Fe β)) (ν^2 (1-Cos[η τ])+η^2 (1-Cos[ν τ]))),0},{0,(E^(Fe β) ν^2 (1-Cos[η τ])+η^2 (1-Cos[ν τ]))/((1+E^(Fe β)) (ν^2 (1-Cos[η τ])+η^2 (1-Cos[ν τ])))}}]

FullSimplify[rhoinv -Flatten[rhoinvpaper]]

(* Note we had to flatten rhoinvpaper so it is in the vector representation we have here. Good, all zeros! *)

(* This is good! We can check that this is indeed an invariant state: trace 1, positive, and invariant under L. *)
Print[FullSimplify[Tr[rhoinv],{Fe >0,Fs>0,β>0,τ>0,λ>0}]];

eigrhoinv = FullSimplify[Eigenvalues[ArrayReshape[rhoinv,{2,2}]],{Fe >0,Fs>0,β>0,τ>0,λ>0}];
Print[eigrhoinv];

(* We want to check that these are positive, so we have a faithful state. *)
(* For the first eigenvalue, we notice the denominator is A from the absolute value before; we know A≤0, so we take the minus sign outside to flip it to non-negative. Then the numerator is already non-negative, as cosine is at most 1, so it is the sum of two non-negative terms. We have two potential problems: the denominator could be zero, or the numerator could be zero. For the numerator, the first term is zero when Cos[ντ]=1, and the second is zero when Cos[ητ]=1. So for the numerator to be zero, we need ντ = π/2 + nπ, and ητ = π/2 + mπ; this gives us a countable set to exclude from our parameter space. By our considerations of A before,we see these are the same conditions for the denominator to vanish. We inspect the second eigenvalue and note that both the numerator and denominator vanish under the same conditions, by the same logic. *)

Print[FullSimplify[LMatrix.rhoinv - rhoinv,{Fe >0,Fs>0,β>0,τ>0,λ>0}]]

(* We obtain {0,0,0,0} as we want.  *)

(* Now we have a faithful invariant state, given our countable set of assumptions. Note also that for any parameters {Fe,Fs,β,λ,τ}, we can plot the eigenvalues of LMatrix, and of LMatrix easily, since we have symbolic expressions for them: *)

With[{rr= eigsys[[1]]},Manipulate[Show[ListPlot[{Re[#],Im[#]}&/@rr,AxesOrigin->{0,0},PlotRange->{{-1.2,1.2},{-1.2,1.2}},ImagePadding->55,AspectRatio->1,Frame->True,FrameLabel->{{Im,None},{Re,"Eigenvalues of LMatrix"}},BaseStyle->{FontSize->Scaled[0.05]},AxesStyle->AbsoluteThickness[.5],ImageSize->Scaled[.3],PlotStyle->Directive[Red,PointSize[.02]]],Graphics[{AbsoluteThickness[.5],Circle[{0,0},1]}]],{{λ,2},0,20},{{Fe,.8},0,2},{{Fs,.9},0,2},{{τ,.5},0,1}]]

With[{rr= eigrhoinv},Manipulate[Show[ListPlot[{Re[#],Im[#]}&/@rr,AxesOrigin->{0,0},PlotRange->{{-1.2,1.2},{-1.2,1.2}},ImagePadding->55,AspectRatio->1,Frame->True,FrameLabel->{{Im,None},{Re,"Eigenvalues of rhoinv"}},BaseStyle->{FontSize->Scaled[0.05]},AxesStyle->AbsoluteThickness[.5],ImageSize->Scaled[.3],PlotStyle->Directive[Red,PointSize[.02]]],Graphics[{AbsoluteThickness[.5],Circle[{0,0},1]}]],{{λ,2},0,20},{{Fe,.8},0,2},{{Fs,.9},0,2},{{τ,.5},0,1},{{β,20},0.01,50}]]

(* Okay, that allows us to easily confirm our spectral hypothesis. Now, let us proceed to compute the entropy to lowest order in the SCL. To this end, we first need to compute Mk. This computation depends on the spectral properties of h0. Recall h0 is degenerate iff Fe = Fs: *)

Print[MatrixForm[h0]]
(* So we need to divide our computation of Mk into the degenerate case (Mkdeg) and the non-degenerate case (Mk). First the non-degenerate case. We find the spectral data of h0: *)

h0esys = Eigensystem[h0];

h0eigproj={};
AppendTo[h0eigproj, Outer[Times,h0esys[[2]][[1]],h0esys[[2]][[1]]]];
AppendTo[h0eigproj, Outer[Times,h0esys[[2]][[2]],h0esys[[2]][[2]]]];
AppendTo[h0eigproj, Outer[Times,h0esys[[2]][[3]],h0esys[[2]][[3]]]];
AppendTo[h0eigproj, Outer[Times,h0esys[[2]][[4]],h0esys[[2]][[4]]]];
(* Now h0eigproj contains the spectral projectors of h0, assuming Fe ≠ Fs *)

(* We wish to compute M_k, which is given by a commutator of the zeroth order term of the invariant state tensor the environment state with a sum of terms, as shown in the paper. We'll calculate that sum here. *) 
sum=0;
For[i=1,i≤ dimhe*dimhs,i++,
For[j=1,j≤dimhe*dimhs,j++,
sum=If[i≠j,
sum+h0eigproj[[i]].v.h0eigproj[[j]]*(Exp[-I*τ*(h0esys[[1]][[i]]-h0esys[[1]][[j]])] -1)/(h0esys[[1]][[i]]-h0esys[[1]][[j]]),sum+
h0eigproj[[i]].v.h0eigproj[[i]]*(-Iτ)];
]
]

(* Now the joint state of interest: *)
omega = KroneckerProduct[ArrayReshape[rhobetaenv,{2,2}],ArrayReshape[rhoinv,{2,2}]];

(* We check if there is a first order term in λ: *)
Print[Series[omega,{λ,0,1}]]

(* There is not, so we just need the zeroth order term: *)
(* omega0=FullSimplify[Limit[omega,λ->0],{Fe>0,Fs>0,τ>0,β>0,λ>0}]; *)
omega0=Limit[omega,λ->0];

Mk = (omega0.sum - sum.omega0)/λ;

(* This is the first order correction to the state in the non-degenerate case *)

(* Now the degenerate case. We repeat the same steps. *)

h0eigprojdeg={};
AppendTo[h0eigprojdeg,h0eigproj[[1]]];
AppendTo[h0eigprojdeg,h0eigproj[[2]]+h0eigproj[[3]]];
AppendTo[h0eigprojdeg,h0eigproj[[4]]];

h0eigdeg={};
AppendTo[h0eigdeg,h0esys[[1]][[1]]];
AppendTo[h0eigdeg,h0esys[[1]][[2]]];
AppendTo[h0eigdeg,h0esys[[1]][[4]]];
h0eigdeg = h0eigdeg/.{Fe->Fs};

(* Now we find the sum of the terms to commute with *)
sumdeg=0;
For[i=1,i<4,i++,
For[j=1,j<4,j++,
sumdeg=If[i≠j,
sumdeg+h0eigprojdeg[[i]].v.h0eigprojdeg[[j]]*(Exp[-I*τ*(h0eigdeg[[i]]-h0eigdeg[[j]])] -1)/(h0eigdeg[[i]]-h0eigdeg[[j]]),sumdeg+
h0eigprojdeg[[i]].v.h0eigprojdeg[[j]]*(-I*τ)];
]
]

(* Still no first order term:  *)
Print[Series[omega/.{Fe-> Fs},{λ,0,1}]]

(* So we compute the zeroth order term *)
(* omega0deg=FullSimplify[Limit[(omega/.{Fe->Fs}),λ->0],{Fe>0,Fs>0,τ>0,β>0,λ>0}]; *)
omega0deg=Limit[(omega/.{Fe->Fs}),λ->0];

(* Mkdeg= FullSimplify[(sumdeg.omega0deg-omega0deg.sumdeg)/λ,{Fe>0,Fs>0,τ>0,β>0,λ>0}]; *)
Mkdeg=(omega0deg.sumdeg-sumdeg.omega0deg)/λ;

(* And there we have it. *)

(* Now we wish to use Mk to compute sigma. First the non-degenerate case. We need the spectral data of omega0. We note that it is diagonal, so we already know the eigenprojectors*)

omega0proj = h0eigproj;

(* And we use the eigenvalues are omega0[[i]][[i]]. *)
entropyUnsimp=0;
For[i=1,i<5,i++,
For[j=i,j<5,j++,
entropyUnsimp=If[i≠j,
entropyUnsimp+Tr[Mk.omega0proj[[j]].Mk.omega0proj[[i]]]*((Log[omega0[[i]][[i]]]-Log[omega0[[j]][[j]]])/(omega0[[i]][[i]]-omega0[[j]][[j]])),entropyUnsimp+Tr[Mk.omega0proj[[i]].Mk.omega0proj[[i]]] / (2 omega0[[i]][[i]])];
]
]

(* Now we simplify: *)
 entropy= FullSimplify[entropyUnsimp,{Fe >0,Fs>0,β>0,τ>0 }];

Print[entropy]

(* Great! Now the degenerate case. Note the eigenprojectors for omega0 in the degenerate case are still the eigenprojectors for h0 in the non-degenerate case; omega0 does not become degenerate.*)

omega0degproj = h0eigproj;

Print[MatrixForm[omega0deg]]

entropydegUnsimp=0;
For[i=1,i<5,i++,
For[j=i,j<5,j++,
entropydegUnsimp=If[i≠j,
entropydegUnsimp+Tr[Mkdeg.omega0degproj[[j]].Mkdeg.omega0degproj[[i]]]*((Log[omega0deg[[i]][[i]]]-Log[omega0deg[[j]][[j]]])/(omega0deg[[i]][[i]]-omega0deg[[j]][[j]])),entropydegUnsimp+Tr[Mkdeg.omega0degproj[[i]].Mkdeg.omega0degproj[[i]]] / (2 omega0deg[[i]][[i]])];
]
]

entropydeg= FullSimplify[entropydegUnsimp,{Fe >0,Fs>0,β>0,τ>0, λ>0}];
Print[entropydeg]

(* We note that the non-degenerate entropy converges to the degenerate entropy in the limit Fe->Fs. *)

difference = FullSimplify[entropydeg - Limit[entropy, Fe->Fs],{Fe >0,Fs>0,β>0,τ>0,λ>0}]