fϕθ=Sin[θ]^2(3+(Cos[θ]^2-Sin[θ]^2) - 2 (Cos[ϕ]^2-Sin[ϕ]^2)Sin[θ]^2)/16
%/.{Cos[θ]->z,Sin[θ]->Sqrt[(1-z^2)],Cos[ϕ]->x/Sqrt[(1-z^2)],Sin[ϕ]->y/Sqrt[(1-z^2)]};
FullSimplify[%];
FullSimplify[%,Assumptions->x^2+y^2+z^2==1]

D[Sqrt[Sin[θ]^2(1-Cos[ϕ]^2Sin[θ]^2)/4],ϕ]/.{ϕ->90*(2 Pi/360),θ-> 90*(2Pi/360)}
D[Sqrt[Sin[θ]^2(1-Cos[ϕ]^2Sin[θ]^2)/4],θ]/.{ϕ->90*(2 Pi/360),θ-> 90*(2Pi/360)}
D[Sqrt[Sin[θ]^2(1-Cos[ϕ]^2Sin[θ]^2)/4],ϕ]/.{ϕ->180*(2 Pi/360),θ-> 45*(2Pi/360)}
D[Sqrt[Sin[θ]^2(1-Cos[ϕ]^2Sin[θ]^2)/4],θ]/.{ϕ->180*(2 Pi/360),θ-> 45*(2Pi/360)}

LittleH[d_,e_,β_,δ_]=blockHamiltonian[[23]];
%//MatrixForm
EpParallel[β_,δ_]=FullSimplify[(Eigenvalues[LittleH[d,e,β,δ]/.{ϕ->0,θ->0}]-HTest[0,0,β,δ][[7]][[7]])/.{d->0.1,e->-0.1/4}][[2]];
EmParallel[β_,δ_] =FullSimplify[(Eigenvalues[LittleH[d,e,β,δ]/.{ϕ->0,θ->0}]-HTest[0,0,β,δ][[7]][[7]])/.{d->0.1,e->-0.1/4}][[1]];
EpBent[β_,δ_]=FullSimplify[(Eigenvalues[LittleH[d,e,β,δ]/.{ϕ->Pi/2,θ->Pi/2}]-HTest[Pi/2,Pi/2,β,δ][[7]][[7]])/.{d->0.1,e->-0.1/4}][[2]];
EmBent[β_,δ_] =FullSimplify[(Eigenvalues[LittleH[d,e,β,δ]/.{ϕ->Pi/2,θ->Pi/2}]-HTest[Pi/2,Pi/2,β,δ][[7]][[7]])/.{d->0.1,e->-0.1/4}][[1]];
Ep=FullSimplify[(Eigenvalues[LittleH[d,e,β,δ]])][[2]]
Em=FullSimplify[(Eigenvalues[LittleH[d,e,β,δ]])][[1]]

Eigensystem[LittleH[d,e,β,δ]/.{ϕ->Pi/2,θ->Pi/2}]/.{d->0.1,e->-0.1/3}//MatrixForm
Eigensystem[LittleH[d,e,β,δ]/.{ϕ->Pi/2,θ->Pi/2}]/.{d->0.1,e->-0.1/4}//MatrixForm

EpBentList=Re[Table[EpBent[β,δ],{δ,0.9,1.1,0.01},{β,0,Pi,0.01}]];
EmBentList=Re[Table[EmBent[β,δ],{δ,0.9,1.1,0.01},{β,0,Pi,0.01}]];
ListPlot3D[{EpBentList,EmBentList}]

ListPlot3D[EmBentList-Table[HTest[Pi/2,Pi/2,β,δ][[7]][[7]]/.{d->0.1,e->-0.1/4},{δ,0.9,1.1,0.01},{β,0,Pi,0.01}]]
ListPlot3D[EpBentList-Table[HTest[Pi/2,Pi/2,β,δ][[7]][[7]]/.{d->0.1,e->-0.1/4},{δ,0.9,1.1,0.01},{β,0,Pi,0.01}]]
ListPlot3D[EpBentList-EmBentList]

Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/Data Files/21_1n1_avoidedCrossing_3Ddata_Bparallely_Ep.txt",Transpose[EpBentList],"Table"];
Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/Data Files/21_1n1_avoidedCrossing_3Ddata_Bparallely_Em.txt",Transpose[EmBentList],"Table"];

EpParallelList=Table[EpParallel[β,δ],{δ,0.9,1.1,0.01},{β,0,Pi,0.01}];
EmParallelList=Table[EmParallel[β,δ],{δ,0.9,1.1,0.01},{β,0,Pi,0.01}];
ListPlot3D[{EpParallelList,EmParallelList}]

Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/Data Files/21_1n1_avoidedCrossing_3Ddata_Bparallelz_Ep.txt",Transpose[EpParallelList],"Table"];
Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/Data Files/21_1n1_avoidedCrossing_3Ddata_Bparallelz_Em.txt",Transpose[EmParallelList],"Table"];

U=Transpose[{{1,0,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0,0},{0,0,1,0,0,0,0,0,0},{0,0,0,I/Sqrt[2],1/Sqrt[2],0,0,0,0},{0,0,0,-I/Sqrt[2],1/Sqrt[2],0,0,0,0},{0,0,0,0,0,1,0,0,0},{0,0,0,0,0,0,1,0,0},{0,0,0,0,0,0,0,1,0},{0,0,0,0,0,0,0,0,1}}];
%//MatrixForm
FullSimplify[Chop[Transpose[U].HTest[0,0,0,0.7].U/.{d->0.1,e->0}]]//MatrixForm

Abs[{0,0,0,0,0,0,1,0,0}.(sAx+sBx).{0,0,0,I/Sqrt[2],0,1/Sqrt[2],0,0,0}]^2

Abs[{0,0,0,0,0,0,1,0,0}.(sAy+sBy).{0,0,-I/Sqrt[2],0,1/Sqrt[2],0,0,0,0}]^2

Abs[HTest[ϕ,θ,Pi/2,δ][[1]][[7]]/.{d->0.1,e->-0.1/30}]^2

FullSimplify[HTest[Pi/2,Pi/2,β,δ]/.{d->0.1,e->-0.1/3}]//MatrixForm

ContourPlot[Abs[HTest[ϕ,θ,Pi/2,δ][[1]][[7]]/.{d->0.1,e->-0.1/30}]^2,{ϕ,0,2Pi},{θ,0,Pi},PlotLegends->Automatic]//MatrixForm

HTest[Pi/2,Pi/2,β,δ]//MatrixForm

ΔEpm[δ_]=FullSimplify[(Ep-Em),Assumptions->{β ∈ Reals,θ ∈ Reals,d>e,δ ∈ Reals,d ∈ Reals, e ∈ Reals}]/.{Cos[2θ]->Cos[θ]^2-Sin[θ]^2,Cos[2ϕ]->Cos[ϕ]^2-Sin[ϕ]^2,Cos[2β]->Cos[β]^2-Sin[β]^2,Cos[4θ]->TrigExpand[Cos[4θ]]}
ΔEpmAdd[δ_]=FullSimplify[(Ep+Em)-2HTest[ϕ,θ,β,δ][[7]][[7]],Assumptions->{β ∈ Reals,θ ∈ Reals,d>e,δ ∈ Reals,d ∈ Reals, e ∈ Reals}]/.{Cos[2θ]->Cos[θ]^2-Sin[θ]^2,Cos[2ϕ]->Cos[ϕ]^2-Sin[ϕ]^2,Cos[2β]->Cos[β]^2-Sin[β]^2,Cos[4θ]->TrigExpand[Cos[4θ]]}
ΔEp0[δ_]=FullSimplify[(Ep-HTest[ Pi/2,Pi/2,β,δ][[7]][[7]]),Assumptions->{β ∈ Reals,θ ∈ Reals,d>e,δ ∈ Reals,d ∈ Reals, e ∈ Reals}]/.{Cos[2θ]->Cos[θ]^2-Sin[θ]^2,Cos[2ϕ]->Cos[ϕ]^2-Sin[ϕ]^2,Cos[2β]->Cos[β]^2-Sin[β]^2,Cos[4θ]->TrigExpand[Cos[4θ]]}
ΔEm0[δ_]=FullSimplify[(Em-HTest[Pi/2,Pi/2,β,δ][[7]][[7]]),Assumptions->{β ∈ Reals,θ ∈ Reals,d>e,δ ∈ Reals,d ∈ Reals, e ∈ Reals}]/.{Cos[2θ]->Cos[θ]^2-Sin[θ]^2,Cos[2ϕ]->Cos[ϕ]^2-Sin[ϕ]^2,Cos[2β]->Cos[β]^2-Sin[β]^2,Cos[4θ]->TrigExpand[Cos[4θ]]}

FullSimplify[FullSimplify[ΔEpmAdd[δ]/(2(d-e))//.{Cos[θ]->z,Sin[θ]->Sqrt[1-z^2],Cos[ϕ]->x/Sqrt[1-z^2],Sin[ϕ]->y/Sqrt[1-z^2],x^2->1-y^2-z^2},Assumptions->{β ∈ Reals,θ ∈ Reals,d>e,δ ∈ Reals,d ∈ Reals, e ∈ Reals}]/.{y->Sin[ϕ]Sin[θ],z->Cos[θ]}]

FullSimplify[%*(d-e)==1 +(d+3e)/12(1-3Sin[θ]^2Sin[ϕ]^2)-(d-e)Cos[β](Cos[2 θ] + Sin[θ]^2Sin[ϕ]^2)/4]

FullSimplify[-Cos[2θ]-Sin[θ]^2Sin[ϕ]^2==Sin[θ]^2(1+Cos[ϕ]^2)-1]

ΔEm0xyz[δ_] =FullSimplify[ΔEm0[δ]//.{Cos[θ]->z,Sin[θ]->Sqrt[1-z^2],Cos[ϕ]->x/Sqrt[1-z^2],Sin[ϕ]->y/Sqrt[1-z^2],x^2->1-y^2-z^2},Assumptions->{β ∈ Reals,θ ∈ Reals,d>e,δ ∈ Reals,d ∈ Reals, e ∈ Reals}]

FullSimplify[ΔEm0[δ]==1+(d+3e)/12(1-3Sin[θ]^2Sin[ϕ]^2)+(d-e)Cos[β](Sin[θ]^2(1+Cos[ϕ]^2)-1)/4 - (d-e)/2Sqrt[Sin[β]^2Sin[θ]^2(1-Sin[θ]^2Cos[ϕ]^2)+4(δ-1)^2/(d-e)^2],Assumptions->{β ∈ Reals,θ ∈ Reals,d>e,δ ∈ Reals,d ∈ Reals, e ∈ Reals}]

FullSimplify[ΔEp0[δ]==1+(d+3e)/12(1-3Sin[θ]^2Sin[ϕ]^2)+(d-e)Cos[β](Sin[θ]^2(1+Cos[ϕ]^2)-1)/4 + (d-e)Sqrt[Sin[β]^2Sin[θ]^2(1-Sin[θ]^2Cos[ϕ]^2)/4+(δ-1)^2/(d-e)^2],Assumptions->{β ∈ Reals,θ ∈ Reals,d>e,δ ∈ Reals,d ∈ Reals, e ∈ Reals}]

(*Calculate the 1st derivative of the new energy curves to find stationary points. To simplify use δ=1 where the on diagonal components of the 2x2 matrix are degenerate*)
FullSimplify[D[{Ep[d,e,β,δ]-Em[d,e,β,δ]},β]]
%/.δ->1

(* Calculate the 2nd derivative of the energy splitting to characterize the curvature *)
FullSimplify[D[D[Ep[d,e,β,ϵ]-Em[d,e,β,ϵ],ϵ],ϵ]/.β->Pi/2]

Look at the Hessian to get a feel for the curvature of the space

Hess=FullSimplify[Det[{{D[D[Ep[d,e,β,δ]-Em[d,e,β,δ],δ],δ],D[D[Ep[d,e,β,δ]-Em[d,e,β,δ],δ],β]},{D[D[Ep[d,e,β,δ]-Em[d,e,β,δ],β],δ],D[D[Ep[d,e,β,δ]-Em[d,e,β,δ],β],β]}}]/.{θ->0,ϕ->0}]

Plot3D[Hess/.{d->0.1,e->0.1/3},{δ,0.9,1.1},{β,0,Pi},AxesLabel->{Style["B/J",20,Black],Style["β",20,Black],Style["",20,Black]},PlotRange->All,Boxed->False,ImageSize->{800,400}]
Export["C:\\Users\\Sina\\Google Drive\\Research\\Singlet Fission Project\\21_1n1_DetHessian_3Dimage.png",%,ImageResolution->500];

Make a 2D plot

diagEnergies[d_,e_,β_,δ_] = Diagonal[H[ϕ,θ,β,δ]];
params={d,e,β,δ}/.{d->0.1,e-> 0.1/3,β->Pi/2,δ->δ};
Plot[Evaluate[{diagEnergies[params[[1]],params[[2]],params[[3]],params[[4]]][[4]],diagEnergies[params[[1]],params[[2]],params[[3]],params[[4]]][[6]],Ep[params[[1]],params[[2]],params[[3]],params[[4]]],Em[params[[1]],params[[2]],params[[3]],params[[4]]]}/.{θ->Pi/2,ϕ->Pi/2}],{δ,0,10},PlotLabel->Style["",Black,28],PlotRange->{{0.9,1.1},{-0.1,0.2}},PlotLegends->Placed[Evaluate[{plotLegend[[4]],plotLegend[[6]],Style["Antisymmetric adiabat",20],Style["Symmetric adiabat",20]}],{1,0.6}],Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy",Black,28],None},{Style["B",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {900,500},AspectRatio->Full,PlotStyle->{(*{Green, Dashed,Thick},{Blue,Thick},{Blend[{Green,Blue},0.4],Dashed,Thick},{Blend[{Green,Blue},0.6],Thick}*){Black,Thick},{Gray,Thick}}]
(*Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/clockTransitionEx.png",%,ImageResolution->500]*)

eigs1

export=True;
eigs1=Eigenvalues[LittleH[0.1,0.1/3,β,δ]]/.{ϕ->Pi/2,θ->Pi/2,β->Pi};
eigs2=Eigenvalues[LittleH[0.1,0.1/3,β,δ]]/.{ϕ->Pi/2,θ->Pi/2,β->Pi/2};
Plot[{eigs1[[1]],eigs1[[2]],eigs2[[1]],eigs2[[2]]},{δ,0.9,1.1},PlotLabel->Style["",Black,28],PlotRange->{{0.9,1.1},{-0.10,0.15}}(*,PlotLegends->Placed[plotLegend,{1.1,0.6}]*), Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy/J",Black,28],None},{Style["Bz/J",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {600,400},AspectRatio->Full,PlotStyle->{{Green,Thick},{Green,Dashed,Thick},{Blue,Thick},{Blue,Thick,Dashed}}(*,Epilog->transitionPoints*)]
If[export,Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/21_1n1_avoidedCrossing_imagev2.png",%,ImageResolution->500],Print["Image not exported"]]

Make a 3D plot

Plot3D[Evaluate[{EpParallel[0.1,-0.1/3,β,δ],EmParallel[0.1,0.1/3,β,δ]}],{δ,0.9,1.1},{β,0,2Pi}]

Plot3D[Evaluate[Eigenvalues[LittleH[0.1,-0.1/3,β,δ]/.{ϕ->Pi/2,θ->Pi/2}]],{δ,0.9,1.1},{β,0,2Pi}]

Plot3D[Evaluate[Eigenvalues[LittleH[0.1,-0.1/3,Pi,δ]]/.{ϕ->ϕt,θ->Pi/2}],{δ,0.9,1.1},{ϕt,0,Pi}]

Plot3D[Evaluate[{Ep[0.1,0.1/3,β,ϵ],Em[0.1,0.1/3,β,ϵ]}/.{θ->0,ϕ->0}],{ϵ,0.9,1.1},{β,0,Pi},BaseStyle->{16,Black,FontFamily->"Avenir"},LabelStyle->Directive[Blend[{Black,White},0.3]],AxesLabel->{Style["(Subscript[B, z])/J",20,Black,FontFamily->"Big Caslon"],Style["β",20,Black,FontFamily->"Big Caslon"],Style["E/J",20,Black,FontFamily->"Big Caslon"]},PlotLabel->Style["B || z",26,Black,FontFamily->"Times"],
PlotStyle->{Directive[RGBColor[0.61,0.17,0.58],Specularity[White,50]],Directive[RGBColor[0.,0.62,1.],Specularity[White,50]]},Boxed->False,ImageSize->{800,400},PlotTheme->"ZMesh",Mesh->15,MeshStyle->{Blend[{Black,White},0.9]}(*,BoundaryStyle->None*),PerformanceGoal->"Quality",Ticks->{Automatic,{0,π/4,π/2,3π/4,π}},FaceGrids->{{{-1,0,0},{{0,Pi/4,Pi/2,3Pi/4,Pi},{-0.05,0,0.05,0.1,0.143}}},{{0,-1,0},{{0.9,0.95,1,1.05,1.1},{-0.05,0,0.05,0.1}}},{{0,0,-1},{{0.9,0.95,1,1.05,1.1},{0,Pi/4,Pi/2,3Pi/4,Pi}}}},FaceGridsStyle->Directive[Black],Background->White,Lighting->"Standard",ViewVertical->{0,0,1},ViewPoint->{10,4,2.5},PlotPoints->50]
Export["/Users/sina/Google Drive/Research/Singlet Fission Project/22_10_avoidedCrossing_3Dimage_Bparallely.png",%,ImageResolution->500]
Plot3D[Evaluate[{Ep[0.1,0.0,β,ϵ],Em[0.1,0.0,β,ϵ]}/.{θ->Pi/2,ϕ->Pi/2}],{ϵ,0.9,1.1},{β,0,Pi},BaseStyle->{16,Black,FontFamily->"Avenir"},LabelStyle->Directive[Blend[{Black,White},0.3]],AxesLabel->{Style["(Subscript[B, z])/J",20,Black,FontFamily->"Big Caslon"],Style["β",20,Black,FontFamily->"Big Caslon"],Style["E/J",20,Black,FontFamily->"Big Caslon"]},PlotLabel->Style["B || z",26,Black,FontFamily->"Times"],
PlotStyle->{Directive[RGBColor[0.61,0.17,0.58],Specularity[White,50]],Directive[RGBColor[0.,0.62,1.],Specularity[White,50]]},Boxed->False,ImageSize->{800,400},PlotTheme->"ZMesh",Mesh->15,MeshStyle->{Blend[{Black,White},0.9]}(*,BoundaryStyle->None*),PerformanceGoal->"Quality",Ticks->{Automatic,{0,π/4,π/2,3π/4,π}},FaceGrids->{{{-1,0,0},{{0,Pi/4,Pi/2,3Pi/4,Pi},{-0.05,0,0.05,0.1,0.143}}},{{0,-1,0},{{0.9,0.95,1,1.05,1.1},{-0.05,0,0.05,0.1}}},{{0,0,-1},{{0.9,0.95,1,1.05,1.1},{0,Pi/4,Pi/2,3Pi/4,Pi}}}},FaceGridsStyle->Directive[Black],Background->White,Lighting->"Standard",ViewVertical->{0,0,1},ViewPoint->{10,4,2.5},PlotPoints->50]
Export["/Users/sina/Google Drive/Research/Singlet Fission Project/22_10_avoidedCrossing_3Dimage_Bparallely.png",%,ImageResolution->500](*
Export["C:\\Users\\Sina\\Google Drive\\Research\\Singlet Fission Project\\21_1n1_avoidedCrossing_3Dimage.png",%,ImageResolution->500];*)

Plot3D[Evaluate[{Ep[0.1,0.1/3,β,δ],Em[0.1,0.1/3,β,δ]}/.{δ->1,ϕ->0}],{θ,0,Pi/2},{β,0,θ},BaseStyle->{16,Black,FontFamily->"Avenir"},LabelStyle->Directive[Blend[{Black,White},0.3]],AxesLabel->{Style["θ",20,Black,FontFamily->"Big Caslon"],Style["β",20,Black,FontFamily->"Big Caslon"],Style["E/J",20,Black,FontFamily->"Big Caslon"]},PlotLabel->Style["Energy Eigenvalues \n (Subscript[B, z])/J = 1, ϕ = 0",26,Black,FontFamily->"Times"],PlotStyle->{Directive[RGBColor[0.61,0.17,0.58],Specularity[White,50]],Directive[RGBColor[0.,0.62,1.],Specularity[White,50]]},Boxed->False,ImageSize->{800,400},PlotTheme->"ZMesh",Mesh->15,MeshStyle->{Blend[{Black,White},0.9]}(*,BoundaryStyle->None*),PerformanceGoal->"Quality",Ticks->{Automatic,{0,π/4,π/2,3π/4,π}},FaceGrids->{{-1,0,0},{0,-1,0},{0,0,-1}},FaceGridsStyle->Directive[Black],Background->White,Lighting->"Standard",ViewVertical->{0,0,1},ViewPoint->{4,10,2.5},PlotPoints->50]

(* Plot the difference of the energy eigenvalues wrt θ,β at δ=1, ϕ=0*)
Plot3D[Evaluate[(Ep[0.1,0.1/3,β,δ]-Em[0.1,0.1/3,β,δ])/.{δ->1,ϕ->0}],{θ,0,Pi/2},{β,0,Pi},BaseStyle->{16,Black,FontFamily->"Avenir"},LabelStyle->Directive[Blend[{Black,White},0.3]],AxesLabel->{Style["θ",20,Black,FontFamily->"Big Caslon"],Style["β",20,Black,FontFamily->"Big Caslon"],Style["E/J",20,Black,FontFamily->"Big Caslon"]},PlotLabel->Style["Difference of Energy Eigenvalues \n (Subscript[B, z])/J = 1, ϕ = 0",26,Black,FontFamily->"Times"],PlotStyle->{Directive[RGBColor[0.61,0.17,0.58]],Directive[RGBColor[0.,0.62,1.]]},Boxed->False,ImageSize->{800,400},PlotTheme->"ZMesh",Mesh->15,MeshStyle->{Blend[{Black,White},0.9]}(*,BoundaryStyle->None*),PerformanceGoal->"Quality",Ticks->{Automatic,{0,π/4,π/2,3π/4,π}},FaceGrids->{{-1,0,0},{0,1,0},{0,0,-1}},FaceGridsStyle->Directive[Black],Background->White,Lighting->"Standard",ViewVertical->{0,0,1},PlotPoints->50]

(* Plot the difference of the energy eigenvalues wrt θ,β at δ=1, ϕ=π/2*)Plot3D[Evaluate[(Ep[0.1,0.1/3,β,δ]-Em[0.1,0.1/3,β,δ])/.{δ->1,ϕ->Pi/2}],{θ,0,Pi/2},{β,0,Pi},BaseStyle->{16,Black,FontFamily->"Avenir"},LabelStyle->Directive[Blend[{Black,White},0.3]],AxesLabel->{Style["θ",20,Black,FontFamily->"Big Caslon"],Style["β",20,Black,FontFamily->"Big Caslon"],Style["E/J",20,Black,FontFamily->"Big Caslon"]},PlotLabel->Style["Difference of Energy Eigenvalues \n (Subscript[B, z])/J = 1, ϕ = π/2",26,Black,FontFamily->"Times"],PlotStyle->{Directive[RGBColor[0.61,0.17,0.58]],Directive[RGBColor[0.,0.62,1.]]},Boxed->False,ImageSize->{800,400},PlotTheme->"ZMesh",Mesh->15,MeshStyle->{Blend[{Black,White},0.9]}(*,BoundaryStyle->None*),PerformanceGoal->"Quality",Ticks->{Automatic,{0,π/4,π/2,3π/4,π}},FaceGrids->{{-1,0,0},{0,1,0},{0,0,-1}},FaceGridsStyle->Directive[Black],Background->White,Lighting->"Standard",ViewVertical->{0,0,1},PlotPoints->50]

Manipulate[Plot3D[Evaluate[(Ep[0.1,0.1/3,β,δ]-Em[0.1,0.1/3,β,δ])/.{δ->1,ϕ->ϕtrue}],{θ,0,Pi/2},{β,0,Pi},BaseStyle->{16,Black,FontFamily->"Avenir"},LabelStyle->Directive[Blend[{Black,White},0.3]],AxesLabel->{Style["θ",20,Black,FontFamily->"Big Caslon"],Style["β",20,Black,FontFamily->"Big Caslon"],Style["E/J",20,Black,FontFamily->"Big Caslon"]},PlotLabel->Style["Difference of Energy Eigenvalues \n (Subscript[B, z])/J = 1, ϕ = π/2",26,Black,FontFamily->"Times"],PlotStyle->{Directive[RGBColor[0.61,0.17,0.58]],Directive[RGBColor[0.,0.62,1.]]},Boxed->False,ImageSize->{800,400},PlotTheme->"ZMesh",Mesh->15,MeshStyle->{Blend[{Black,White},0.9]}(*,BoundaryStyle->None*),PerformanceGoal->"Quality",Ticks->{Automatic,{0,π/4,π/2,3π/4,π}},FaceGrids->{{-1,0,0},{0,1,0},{0,0,-1}},FaceGridsStyle->Directive[Black],Background->White,Lighting->"Standard",ViewVertical->{0,0,1},PlotPoints->50],{{ϕtrue,0,"ϕ"},0,Pi}]

(* Plot the difference of the energy eigenvalues wrt θ,ϕ at δ=1, β=π/2*)Plot3D[Evaluate[(Ep[0.1,0.1/3,β,δ]-Em[0.1,0.1/3,β,δ])/.{δ->1,β->Pi/2}],{θ,0,Pi},{ϕ,0,Pi},BaseStyle->{16,Black,FontFamily->"Avenir"},LabelStyle->Directive[Blend[{Black,White},0.3]],AxesLabel->{Style["θ",20,Black,FontFamily->"Big Caslon"],Style["ϕ",20,Black,FontFamily->"Big Caslon"],Style["E/J",20,Black,FontFamily->"Big Caslon"]},PlotLabel->Style["Difference of Energy Eigenvalues \n (Subscript[B, z])/J = 1, β = π/2",26,Black,FontFamily->"Times"],PlotStyle->{Directive[RGBColor[0.61,0.17,0.58]],Directive[RGBColor[0.,0.62,1.]]},Boxed->False,ImageSize->{800,400},PlotTheme->"ZMesh",Mesh->15,MeshStyle->{Blend[{Black,White},0.9]}(*,BoundaryStyle->None*),PerformanceGoal->"Quality",Ticks->{Automatic,{0,π/4,π/2,3π/4,π}},FaceGrids->{{-1,0,0},{0,1,0},{0,0,-1}},FaceGridsStyle->Directive[Black],Background->White,Lighting->"Standard",ViewVertical->{0,0,1},PlotPoints->50]

derivative=FullSimplify[D[(Ep[0.1,0.1/3,β,δ]-Em[0.1,0.1/3,β,δ]),β]/.δ->1]

list=NumeratorDenominator[derivative]

ContourPlot3D[list[[1]]==0,{ϕ,0,Pi/2},{θ,0,Pi/2},{β,0,Pi},BaseStyle->{18,Black,FontFamily->"Avenir"},LabelStyle->Directive[Blend[{Black,White},0.3]],AxesLabel->{Style["θ",24,Black,FontFamily->"Avenir"],Style["β",24,Black,FontFamily->"Avenir"],Style["ϕ",24,Black,FontFamily->"Avenir"]},Boxed->False,ImageSize->{800,400},PlotTheme->"ZMesh",Mesh->15,MeshStyle->{Blend[{Black,White},0.9]}(*,BoundaryStyle->None*),PerformanceGoal->"Quality",Ticks->{Automatic,{0,π/4,π/2,3π/4,π}},FaceGrids->{{-1,0,0},{0,-1,0},{0,0,-1}},FaceGridsStyle->Directive[Black],Background->White,Lighting->"Standard",ViewVertical->{0,0,1},ViewPoint->{10,4,2.5},PlotPoints->50]

ContourPlot3D[list[[2]]==0,{ϕ,0,Pi/2},{θ,0,Pi/2},{β,0,Pi},BaseStyle->{18,Black,FontFamily->"Avenir"},LabelStyle->Directive[Blend[{Black,White},0.3]],AxesLabel->{Style["θ",24,Black,FontFamily->"Avenir"],Style["β",24,Black,FontFamily->"Avenir"],Style["ϕ",24,Black,FontFamily->"Avenir"]},Boxed->False,ImageSize->{800,400},PlotTheme->"ZMesh",Mesh->15,MeshStyle->{Blend[{Black,White},0.9]}(*,BoundaryStyle->None*),PerformanceGoal->"Quality",Ticks->{Automatic,{0,π/4,π/2,3π/4,π}},FaceGrids->{{-1,0,0},{0,-1,0},{0,0,-1}},FaceGridsStyle->Directive[Black],Background->White,Lighting->"Standard",ViewVertical->{0,0,1},ViewPoint->{10,4,2.5},PlotPoints->50]

list[[2]]

stepSize=0.01;
dataEm[ϕvar_,θvar_]:=Table[Table[EmParallel[0.1,0.1/3,β,δ]/.{ϕ->ϕvar,θ->θvar},{δ,0.9,1.1,stepSize}],{β,0,Pi,stepSize}];
dataEp[ϕvar_,θvar_]:=Table[Table[EpParallel[0.1,0.1/3,β,δ]/.{ϕ->ϕvar,θ->θvar},{δ,0.9,1.1,stepSize}],{β,0,Pi,stepSize}];
Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/Data Files/21_1n1_avoidedCrossing_3Ddata_Bparallely_Em.txt",dataEm[Pi/2,Pi/2],"Table"];
Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/Data Files/21_1n1_avoidedCrossing_3Ddata_Bparallely_Ep.txt",dataEp[Pi/2,Pi/2],"Table"];
(*Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/Data Files/21_1n1_avoidedCrossing_3Ddata_Bparallelz_Em.txt",dataEm[0,0],"Table"];
Export["/Users/sina/Google Drive/Research/Singlet Fission Project/21_1n1_Analysis/Data Files/21_1n1_avoidedCrossing_3Ddata_Bparallelz_Ep.txt",dataEp[0,0],"Table"];*)

dataEm[ϕ,θ]

export=False;
Plot3D[Evaluate[Ep[0.1,0.1/3,β,ϵ]-Em[0.1,0.1/3,β,ϵ]/.{θ->Pi/2,ϕ->Pi/2}],{ϵ,0.9,1.1},{β,0,Pi},AxesLabel->{Style["B/J",20,Black],Style["β",20,Black],Style["E/J",20,Black]},Boxed->False,ImageSize->{800,400}]
If[export,Export["C:\\Users\\Sina\\Google Drive\\Research\\Singlet Fission Project\\21_1n1_avoidedCrossing_energyDiff_3Dimage.png",%,ImageResolution->500];,Print["Image not saved"]]

Analyze the allowed transitions

(* DIAGONALIZATION OF HAMILTONIAN *)
(* Re-arrange Hamiltonian so |1,-1⟩ and |2,1⟩ are next to each other *)
Uswap = IdentityMatrix[9];
Uswap[[5;;6]]={Uswap[[6]],Uswap[[5]]}
Uswap//MatrixForm
Hswapped=Uswap.H[ϕ,θ,β,δ].Transpose[Uswap];
%//MatrixForm

eigs=Eigenvectors[LittleH[d,e,β,δ]];
soln={FullSimplify[Normalize[eigs[[1]]],Assumptions->{δ∈ Reals,e ∈ Reals,d ∈ Reals,β ∈ Reals,δ>0,β>0,β≤Pi,d>e}],FullSimplify[Normalize[eigs[[2]]],Assumptions->{δ∈ Reals,e ∈ Reals,d ∈ Reals,β ∈ Reals,δ>0,β>0,β≤Pi,d>e}]}

soln/.δ->1

Plot[Evaluate[{soln[[1]][[1]],-soln[[2]][[1]]}/.{d->0.1,e->0.1/3,β->Pi/2,θ->0,ϕ->0}],{δ,0.5,1.5}]

(*Define matrix to take |1,-1⟩ and |2,1⟩ to linear combinations of each other *)
Udiag=IdentityMatrix[9];
Udiag[[4;;5]]={{0,0,0,soln[[1]][[1]],soln[[2]][[1]],0,0,0,0},{0,0,0,soln[[1]][[2]],soln[[2]][[2]],0,0,0,0}}/.{β->Pi/2};
(*Udiag[[4;;5]]=1/Sqrt[2]{{0,0,0,1,-1,0,0,0,0},{0,0,0,1,1,0,0,0,0}};*)
Udiag//MatrixForm

transitionZFS=FullSimplify[Udiag.Hswapped.Transpose[Udiag]/.{δ->1,β->Pi/2}];
%//MatrixForm

A transition due to an externally applied perpendicular magnetic field will be proportional to  |⟨α|(Subscript[S, x])|β⟩((|)^2) a la Fermi’s golden rule.

vecList

vecs={{1,0,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0,0},{0,0,1,0,0,0,0,0,0},{0,0,0,0,1,0,0,0,0},{0,0,0,0,0,0,1,0,0},{0,0,0,0,0,0,0,1,0},{0,0,0,0,0,0,0,0,1},{0,0,0,1/Sqrt[2],0,1/Sqrt[2],0,0,0}}
Table[Abs[vecs[[i]].(UtotSpin.(sAx+sBx).Transpose[UtotSpin]).{0,0,0,1/Sqrt[2],0,1/Sqrt[2],0,0,0}]^2,{i,1,8}]

(UtotSpin.(sAx+sBx).Transpose[UtotSpin])//MatrixForm
(UtotSpin.(sAx+sBx).Transpose[UtotSpin])*(UtotSpin.(sAx+sBx).Transpose[UtotSpin])//MatrixForm
FullSimplify[MatrixPower[(UtotSpin.(sAx+sBx).Transpose[UtotSpin]),2]]//MatrixForm

Possible transitions from |2,0⟩?

vecs={{1,0,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0,0},{0,0,1,0,0,0,0,0,0},{0,0,0,1,0,0,0,0,0},{0,0,0,0,1,0,0,0,0},{0,0,0,0,0,1,0,0,0},{0,0,0,0,0,0,0,1,0},{0,0,0,0,0,0,0,0,1}}
Table[Abs[vecs[[i]].(UtotSpin.(sAz+sBz).Transpose[UtotSpin]).{0,0,0,0,0,0,0,0,0}]^2,{i,1,8}]

(UtotSpin.(sAz+sBz).Transpose[UtotSpin])//MatrixForm

Abs[{0,0,0,0,0,0,0,1,0}.(UtotSpin.(sAz+sBz).Transpose[UtotSpin]).{0,0,0,0,0,0,0,1,0}]^2

FullSimplify[-I(UtotSpin.(sAy+sBy).Transpose[UtotSpin])]//MatrixForm
FullSimplify[(UtotSpin.(sAx+sBx).Transpose[UtotSpin])]//MatrixForm

tensor=sy
T00 = -Tr[tensor]/Sqrt[3]
T10=-I(tensor[[1]][[2]]-tensor[[2]][[1]])/Sqrt[2]
T1p1 = -(tensor[[3]][[1]]-tensor[[1]][[3]]+I(tensor[[3]][[2]]-tensor[[2]][[3]]))/2
T1m1 = -(tensor[[3]][[1]]-tensor[[1]][[3]]-I(tensor[[3]][[2]]-tensor[[2]][[3]]))/2
T20 = (3tensor[[3]][[3]]-Tr[tensor])/Sqrt[6]
T2p1 = -(tensor[[1]][[3]]+tensor[[3]][[1]]+I(tensor[[2]][[3]]+tensor[[3]][[2]]))/2
T2m1 = +(tensor[[1]][[3]]+tensor[[3]][[1]]-I(tensor[[2]][[3]]+tensor[[3]][[2]]))/2
T2p2 = (tensor[[1]][[1]]-tensor[[2]][[2]]+I(tensor[[1]][[2]]+tensor[[2]][[1]]))/2
T2m2 = (tensor[[1]][[1]]-tensor[[2]][[2]]-I(tensor[[1]][[2]]+tensor[[2]][[1]]))/2

≥