LoadSpinChainMatrices[2,"a"];

H = ((-h)/2) σz//mf;

ρ = MatrixExp[-β H]//mf;

Z = Tr[ρ]//Simplify

ρ = ρ/Z//Simplify//mf;

Tr[ρ]//Simplify

Tr[ρ.ρ]//Simplify

H = -(h1/2) sz[1]+(h2/2)sz[2]//mf;

ρ = MatrixExp[β H]//mf;

Z = Tr[ρ]//Simplify

ρ = ρ/Z//mf;

PTr[ρ,{1}]//Simplify//mf;

PTr[ρ,{2}]//Simplify//mf;

H = -(h1/2) sz[1]+(h2/2)sz[2]+ J sz[1].sz[2]//mf;

Sin[10.0]

Sin@10.0

Log[10,1.0]

ρ = ((MatrixExp[-β H])/(Tr[MatrixExp[-β H]]))//Simplify//mf;

ρA = PTr[ρ,{2}]//Simplify//mf;

ρB = PTr[ρ,{1}]//Simplify//mf;

ρrecombined = kron[ρA,ρB]//Simplify//mf;

ρ-ρrecombined//Simplify//mf;

ρ-ρrecombined/.J->0//Simplify//mf;

H = g (sp[1].sm[2]+sm[1].sp[2])//mf;

U = MatrixExp[-I t H]//Simplify//mf;

ρt=U.ρ.U\[ConjugateTranspose]//cf//mf;

PTr[ρt,{1}]//Simplify//mf;

LoadSpinChainMatrices[L=6];

H = Sum[sz[i].sz[i+1],{i,1,L-1}];
ρ = MatrixExp[- H];
ρ = ρ/Tr[ρ]//mf;

Dimensions[ρ]

PTr[ρ,{1,2,3,4,5}]//mf;

LoadBosonicOperators[nmax = 20]

LoadedMatrices

(H = kron[a\[ConjugateTranspose].a,σ0] + kron[𝕀,σz] + kron[a+a\[ConjugateTranspose],σx])//MatrixPlot

ρ = MatrixExp[-H];
(ρ = ρ/Tr[ρ])//MatrixPlot

Tr[ρ.kron[a\[ConjugateTranspose].a,σ0]]

(ρqho = PTr[ρ,{2},{nmax,2}]);

Tr[ρqho.a\[ConjugateTranspose].a]