submitted version (141128)
***VERSION SD3.10.2 IS REQUIRED***

Needs["SpinDynamica`"]
SetSpinSystem[2];
SetBasis[SingletTripletBasis[]]
SetOperatorBasis[CartesianProductOperatorBasis[]]

general

Symbols

Format[Ω1]=(Subscript[Ω, 1])
Format[Ω2]=(Subscript[Ω, 2])
Format[ωΔ]=(Subscript[ω, "Δ"])
Format[ωΣ]=(Subscript[ω, "Σ"])
Format[ωJ]=(Subscript[ω, "J"])
Format[θST]=(Subscript[θ, "ST"])

Format[ωnut]=(Subscript[ω, "nut"])

Format[Δωnut]=(Subscript[Δω, "nut"])

Format[θe]=(Subscript[θ, "e"])
Format[ωe]=(Subscript[ω, "e"])

Format[ϕB]=(Subscript[ϕ, "B"])
Format[ϕC]=(Subscript[ϕ, "C"])

Format[τ1]=(Subscript[τ, "1"])
Format[τ2]=(Subscript[τ, "2"])

Format[ϵrf]=(Subscript[ϵ, "rf"])

Format[κp]=(κ^("'"))

Hamiltonians

HJ

HJ=ωJ opI[1].opI[2]

HΣ

HΣ=ωΣ (opI[1,"z"]+opI[2,"z"])/2

HΔ

HΔ=ωΔ(opI[1,"z"]-opI[2,"z"])/2

INADEQUATE

Hamiltonian

factor out overall offset

H=(1/2)ωΔ (opI[1,"z"]- opI[2,"z"])+ ωJ opI[1].opI[2]

ρini

ρini=opI["z"]

calculate ρf

ExpressOperator[RotationSuperoperator[{π/2,"y"}][ρini]]

ρf=Fold[
Operator[MatrixRepresentation[#2[#1]]]&,
ρini,
{RotationSuperoperator[{π/2,"y"}],
Exp[-I τ1 CommutationSuperoperator[H]],
RotationSuperoperator[{π,"x"}],
Exp[-I τ1 CommutationSuperoperator[H]],
RotationSuperoperator[{π/2,"y"}]
}
]

ρf=Fold[
Operator[Together@Simplify@PowerExpand@Normal[MatrixRepresentation[#2[#1]]/.ωJ^2+ωΔ^2->Ω^2]]&,
ρini,
{RotationSuperoperator[{π/2,"y"}],
Exp[-I τ1 CommutationSuperoperator[H]],
RotationSuperoperator[{π,"x"}],
Exp[-I τ1 CommutationSuperoperator[H]],
RotationSuperoperator[{π/2,"y"}]
}
]

MatrixRepresentation[ρf]

DQ amplitude

ρf

aDQ["INADQ"]=Simplify[Simplify[ExpToTrig@MatrixRepresentation[ρf][[2,4]]]/.{Ω^2->ωΔ^2+ωJ^2}//.{ωJ->Ω Cos[θST],ωΔ-> Ω Sin[θST]}]

as in eq.(2) [vs 251028]

DQF amplitude

aDQF["INADQ"]=ComplexExpand[aDQ["INADQ"]*Conjugate[aDQ["INADQ"]]]//Simplify

INADEQUATE density plot vs θ and τ, as in the paper

aDQF[“INADQ”]

aDQF["INADQ"]

INADEQUATE density plot vs θ and τ, as in the paper

aDQF["INADQ"]

aDQF["INADQ"]//.{τ1->Jτ1/J,J->ωJ/(2π),ωJ->Ω Cos[θST],θST->N[θSTdeg °]}

DensityPlot[
Evaluate[aDQF["INADQ"]//.{τ1->Jτ1/J,J->ωJ/(2π),ωJ->Ω Cos[θST],θST->N[θSTdeg °]}],
{θSTdeg,0,90},{Jτ1,0,2},
FrameLabel->{{"J τ1",None},{"(Subscript[θ, ST]) / °",None}},
LabelStyle->16,
PlotLegends->Automatic,
PlotPoints->100,
PlotLabel->"INADQ",
ColorFunction->(GrayLevel[1-#]&)
]

DQ excitation schemes

geoDQ

ρ1

ρ1=opI["z"]

ρ2

ρ2=opR[{π/2,"y"}].ρ1.Adjoint@opR[{π/2,"y"}]
ExpressOperator[%]

MatrixRepresentation[ρ2]==MatrixRepresentation[2^(1/2)(SingleTransitionOperator[{2,3},"x"]+SingleTransitionOperator[{3,4},"x"])]

ρ3

ρ3=Exp[-I π SingleTransitionOperator[{1,3},"z"]].ρ2.Adjoint@Exp[-I π SingleTransitionOperator[{1,3},"z"]]
ExpressOperator[%]

ρ4

ρ4=opR[{π/2,"y"}].ρ3.Adjoint@opR[{π/2,"y"}]
ExpressOperator[%]

MatrixRepresentation[%]==MatrixRepresentation[2SingleTransitionOperator[{2,4},"y"]]

MatrixRepresentation[ρ4][[2,4]]

MatrixRepresentation[ρ4][[4,2]]

spinorDQ

ρ1

ρ1=opI["z"]

ρ2

ρ2=opR[{π/2,"y"}].ρ1.Adjoint@opR[{π/2,"y"}]
ExpressOperator[%]

MatrixRepresentation[ρ2]==MatrixRepresentation[2^(1/2)(SingleTransitionOperator[{2,3},"x"]+SingleTransitionOperator[{3,4},"x"])]

ρ3

C12=Exp[-I 2π SingleTransitionOperator[{1,2},"x"]]

ρ3=C12.ρ2.Adjoint@C12
ExpressOperator[%]

ρ4

ρ4=opR[{π/2,"x"}].ρ3.Adjoint@opR[{π/2,"x"}]
ExpressOperator[%]

MatrixRepresentation[%]==MatrixRepresentation[2SingleTransitionOperator[{2,4},"y"]]

MatrixRepresentation[ρ4][[2,4]]

MatrixRepresentation[ρ4][[4,2]]

PulsePol

scaling factors κ

new derivation for PulsePol basic element 90y-τ-180x-τ90y

rf Propagator for basic element of duration τR

Urf0[t_]:=Piecewise[
{{opR[{π/2,"y"}],0≤t≤τR/2},
{opR[{π,"x"}].opR[{π/2,"y"}],τR/2≤t≤τR}}
]

Urf0[t]

derive rf Euler

U=opR[{π/2,"y"}];

RoundToSubmultipleOfπ[ang_]:=RoundToSubmultipleOfπ[ang,36]
RoundToSubmultipleOfπ[ang_,n_]:=Round[ang/(π/n)]*(π/n)

RoundToSubmultipleOfπ[3.1]

OperatorToEuler[U_]:=Map[RoundToSubmultipleOfπ,
-Reverse@AxesToEuler@Outer[OperatorAmplitude[U.#1.Adjoint[U]->#2]&,{opI["x"],opI["y"],opI["z"]},{opI["x"],opI["y"],opI["z"]}]]

OperatorToEuler[opR[{π/2,"y"}]]
N[%/°]

OperatorToEuler[opR[{π/2,"y"}]]
N[%/°]

OperatorToEuler[opR[{π,"x"}].opR[{π/2,"y"}]]
N[%/°]

rf Euler

Ωrf0[t_]:=
Piecewise[
{{OperatorToEuler[opR[{π/2,"y"}]],0≤t≤τR/2},
{OperatorToEuler[opR[{π,"x"}].opR[{π/2,"y"}]],τR/2≤t≤τR}}
]

αrf0[t_]:=MapAt[Part[#,1]&,Ωrf0[t],{{1,1,1},{1,2,1}}]

βrf0[t_]:=MapAt[Part[#,2]&,Ωrf0[t],{{1,1,1},{1,2,1}}]

γrf0[t_]:=MapAt[Part[#,3]&,Ωrf0[t],{{1,1,1},{1,2,1}}]

Ωrf0[t]

αrf0[t]

βrf0[t]

γrf0[t]

K factors (eq.37 of Sabba 2022)

K[m_,λ_,μ_]:=τR^(-1)Integrate[
Wignerd[λ,{μ,0}][-βrf0[t]]
Exp[I(μ γrf0[t]+m ωJ t)],
{t,0,τR},Assumptions->{τR∈Reals,τR>0}
]

K[1,1,1]

K[1,1,1]//.{τR->n τJ/bigN,τJ->2π/ωJ}

Simplify[K[1,1,1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ}]

Arg[%]//Simplify

Exp[-I μ (π ν/bigN)]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ,μ->1}

Arg[%]//Simplify

κ factors (eq.36 of Sabba 2022)

Clear[κ]

κ[l_,m_,λ_,μ_]:=Exp[-I μ (π ν/bigN)]K[m,λ,μ]

Exp[-I μ (π ν/bigN)]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ,μ->1}

Arg[%]//Simplify

Simplify[K[1,1,1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ}]
Arg[%]

Simplify[(Exp[-I μ (π ν/bigN)]K[1,1,1])//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ,μ->1}]
Arg[%]

Simplify[κ[1,1,1,1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ,μ->1}]

κ[1,1,1,1]

κ["R431"]=Simplify[κ[1,1,1,1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ}]

Abs[κ["R431"]]//N

Arg[κ["R431"]]//Simplify

Simplify[κ[1,-1,1,-1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ}]

propagation with R431

ρ1

ρ1=opI["z"]

MatrixRepresentation[ρ1]==
MatrixRepresentation[2 SingleTransitionOperator[{2,4},"z"]]

ρ2

ρ2=opI["x"]

MatrixRepresentation[ρ2]==
MatrixRepresentation[
2^(1/2)(
SingleTransitionOperator[{2,3},"x"]
+SingleTransitionOperator[{3,4},"x"]
)
]

MatrixRepresentation[ρ2]==
MatrixRepresentation[
2^(-1/2)(
Ket[2].Bra[3]+Ket[3].Bra[2]
+Ket[3].Bra[4]+Ket[4].Bra[3]
)
]

ρ3

Rx12[β_]:=Exp[-I β SingleTransitionOperator[{1,2},"x"]]

ρ3=Rx12[β].ρ2.Adjoint[Rx12[β]]

ExpressOperator[ρ3,CartesianProductOperatorBasis[]]

Normal@MatrixRepresentation[ρ3]==
Normal@MatrixRepresentation[
2^(1/2)(
(SingleTransitionOperator[{2,3},"x"]Cos[β/2]+SingleTransitionOperator[{1,3},"y"]Sin[β/2])
+SingleTransitionOperator[{3,4},"x"]
)
]

ρ4

ρ4=ExpressOperator[opR[{π/2,"x"}].ρ3.Adjoint[opR[{π/2,"x"}]]]

OperatorAmplitude[ρ4->SingleTransitionOperator[{2,4},"y"]]

aDQ[“PulsePol”]

aDQ["PulsePol"]=Simplify[MatrixRepresentation[ρ4][[2,4]]]

aDQF[“PulsePol”]

aDQF["PulsePol"]=ComplexExpand[Abs[aDQ["PulsePol"]]^2]

SLIC

terms

W (eq.76 [251029])

W[t_]:= opR[{π/2,"-y"}].Exp[-I HJ t].opR[{-ωJ t,"z"}]

W[t]

HΔtilde (eq.83 [251029])

HΔtilde[t_]:=Adjoint[W[t]].HΔ.W[t]

HΔtilde[t]

OperatorAmplitude[HΔtilde[t]->#]&/@
{SingleTransitionOperator[{1,2},"x"],SingleTransitionOperator[{1,4},"x"],SingleTransitionOperator[{1,4},"y"]}

ExpToTrig[%]//Simplify

HΣtilde (eq.83 [251029])

HΣtilde[t_]:=Adjoint[W[t]].HΣ.W[t]

HΣtilde[t]

OperatorAmplitude[HΣtilde[t]->#]&/@
{opI["x"],opI["y"],SingleTransitionOperator[{1,4},"y"]}

ExpToTrig[%]//Simplify

Hbar (eq.85 [251029])

Hbar=-κ ωΔ SingleTransitionOperator[{1,2},"x"]-ϵrf ωJ opI["z"]

ExpressOperator@Hbar

SLIC(a)

eq.87 [251029]

Adjoint@W[0].opI["z"].W[0]

ExpressOperator[%]

MatrixRepresentation[%]==MatrixRepresentation[2^(1/2)(SingleTransitionOperator[{2,3},"x"]+SingleTransitionOperator[{3,4},"x"])]

eqs.90 to 92 [251029]

R12=Exp[-I Hbar T]/.{ϵrf->0,T->-β12/(ωΔ κ)}

OperatorAmplitude[
R12.opI["x"].Adjoint[R12]->#
]&/@{SingleTransitionOperator[{2,3},"x"],SingleTransitionOperator[{1,3},"y"],SingleTransitionOperator[{3,4},"x"]}

OperatorAmplitude[
(R12.opI["x"].Adjoint[R12]/.β12->-2π)->#
]&/@{SingleTransitionOperator[{2,3},"x"],SingleTransitionOperator[{1,3},"y"],SingleTransitionOperator[{3,4},"x"]}

eq.93 [251029]

ρf["SLIC(a)"]=W[T].R12.opI["x"].Adjoint[R12].Adjoint@W[T]

OperatorAmplitude[ρf["SLIC(a)"]->SingleTransitionOperator[{2,4},"y"]]//ExpToTrig//Simplify

eq.94 [251029]

aDQ["SLIC(a)"]=MatrixRepresentation[ρf["SLIC(a)"]][[2,4]]//ExpToTrig//Simplify

SLIC(b)

add a (π/2)x pulse on the  end

ρf["SLIC(b)"]=opR[{π/2,"x"}].ρf["SLIC(a)"].Adjoint[opR[{π/2,"x"}]]

OperatorAmplitude[ρf["SLIC(b)"]->SingleTransitionOperator[{2,4},"y"]]//ExpToTrig//Simplify

aDQ["SLIC(b)"]=MatrixRepresentation[ρf["SLIC(b)"]][[2,4]]//ExpToTrig//Simplify

more

SLIC

general

W transformation

W[t_]:= opR[{π/2,"-y"}].Exp[-I HJ t].opR[{-ωJ t,"z"}]

W[t]

HΔtilde

HΔtilde[t_]:=Adjoint[W[t]].HΔ.W[t]

HΔtilde[t]

OperatorAmplitude[HΔtilde[t]->#]&/@
{SingleTransitionOperator[{1,2},"x"],SingleTransitionOperator[{1,4},"x"],SingleTransitionOperator[{1,4},"y"]}

ExpToTrig[%]//Simplify

HΣtilde

HΣtilde[t_]:=Adjoint[W[t]].HΣ.W[t]

HΣtilde[t]

OperatorAmplitude[HΣtilde[t]->#]&/@
{opI["x"],opI["y"],SingleTransitionOperator[{1,4},"y"]}

ExpToTrig[%]//Simplify

SLIC(a) [with exact rf field]

step 1 (eq.84 [v251028])

ExpressOperator[Adjoint[W[0]].opI["z"].W[0]]

OperatorAmplitude[%->#]&/@
{SingleTransitionOperator[{2,3},"x"],SingleTransitionOperator[{3,4},"x"]}

step 2 (eq.89 [v251028])

β12=-2π

ExpressOperator[Exp[-I β12 SingleTransitionOperator[{1,2},"x"]].opI["x"].Exp[+I β12 SingleTransitionOperator[{1,2},"x"]]]

OperatorAmplitude[%->#]&/@
{SingleTransitionOperator[{2,3},"x"],SingleTransitionOperator[{3,4},"x"]}

step 3

Clear[β12]

ExpressOperator[W[T].Exp[-I β12 SingleTransitionOperator[{1,2},"x"]].opI["x"].Exp[+I β12 SingleTransitionOperator[{1,2},"x"]].Adjoint[W[T]]]

aDQ: eq.90 [v251028]

OperatorAmplitude[%->#]&/@
{SingleTransitionOperator[{2,4},"y"]}

ExpToTrig[%]

Simplify[%]

SLIC(b) [with rf field error]

H1bar eq.82

H1bar=-κ ωΔ SingleTransitionOperator[{1,2},"x"]-ϵrf ωJ opI["z"]

step 1 (eq.84 [v251028])

ExpressOperator[Adjoint[W[0]].opI["z"].W[0]]

OperatorAmplitude[%->#]&/@
{SingleTransitionOperator[{2,3},"x"],SingleTransitionOperator[{3,4},"x"]}

step 2 (eq.89 [v251028]) but with finite ϵ

ExpressOperator[Exp[-I H1bar T].opI["x"].Exp[+I  H1bar T]]

temp=%

temp[[1]]

temp[[1,3,2,3,2]]

temp[[1,3,2,3,2]][[1]]

temp[[1]]//.temp[[1,3,2,3,2]][[1]]->- κp^2 ωΔ^2
PowerExpand[%]

temp//.temp[[1,3,2,3,2]][[1]]->- κp^2 ωΔ^2
temp2=PowerExpand[%]

temp2[[1]]

temp2[[2]]

Series[temp2,{ωΔ,0,1}]

temp3=Normal[%]

temp4=ExpressOperator[W[T].temp3.Adjoint[W[T]]]

MatrixRepresentation[temp4][[2,4]]

%//ExpToTrig//Simplify

%/.{ϵrf->(ωnut-ωJ)/ωJ}

Simplify[%]

more

general

Symbols

Format[Ω1]=(Subscript[Ω, 1])
Format[Ω2]=(Subscript[Ω, 2])
Format[ωΔ]=(Subscript[ω, "Δ"])
Format[ωΣ]=(Subscript[ω, "Σ"])
Format[ωJ]=(Subscript[ω, "J"])
Format[θST]=(Subscript[θ, "ST"])

Format[ωnut]=(Subscript[ω, "nut"])

Format[Δωnut]=(Subscript[Δω, "nut"])

Format[θe]=(Subscript[θ, "e"])
Format[ωe]=(Subscript[ω, "e"])

Format[ϕB]=(Subscript[ϕ, "B"])
Format[ϕC]=(Subscript[ϕ, "C"])

SubstitutionRules

SubstitutionRules1={ωJ^2+ ωΔ^2->Ω^2}
SubstitutionRules2={Ω^2->ωΔ^2+ωJ^2}
SubstitutionRules3={ωJ->Ω Cos[θST],ωΔ-> Ω Sin[θST]}
SubstitutionRules4={Ω->Sqrt[ωJ^2+ωΔ^2],θST->ArcTan[ωJ,ωΔ]}

INADEQUATE for general 2-spin-1/2 system

Hamiltonian

factor out overall offset

H=(1/2)ωΔ (opI[1,"z"]- opI[2,"z"])+ ωJ opI[1].opI[2]

ρini

ρini=opI["z"]

propagation through 3-pulse sequence

calculate and save ρfinal (execution disabled)

ρfinal=Fold[
PowerExpand[ComplexExpand[ExpressOperator[#2[#1],CartesianProductOperatorBasis[]]]/.SubstitutionRules1]&,
ρini,
{RotationSuperoperator[{π/2,"y"}],
Exp[-I τ CommutationSuperoperator[H]],
RotationSuperoperator[{π,"x"}],
Exp[-I τ CommutationSuperoperator[H]],
RotationSuperoperator[{π/2,"y"}]
}
]

SetDirectory[NotebookDirectory[]];
Save["rhofinal-INADQ.m",ρfinal]

read ρfinal

SetDirectory[NotebookDirectory[]];
Read["rhofinal-INADQ.m"]

ρfinal

DQ amplitude

ρfinal

iDQamp1=Simplify[I OperatorAmplitude[ρfinal->opI[1,"+"].opI[2,"+"]]/.SubstitutionRules2]

iDQamp2=Simplify@PowerExpand[iDQamp1/.SubstitutionRules3]

iDQamp3=FullSimplify[iDQamp2/.{ τ Ω Cos[θST]->τ ωJ}]

iDQamp=iDQamp3

iDQamp/.{θST->π/2,Ω->ωΔ}

iDQamp/.{θST->0,Ω->ωJ}

ρfinal

FullSimplify[
Simplify[PowerExpand[
MatrixRepresentation[ρfinal][[2,4]]/.SubstitutionRules2
]/.SubstitutionRules3
]/.{ τ Ω Cos[θST]->τ ωJ}
]

back substitute DQ amplitude and plot

iDQamp4=Simplify[iDQamp/.SubstitutionRules4]

parameters={J->50,ωJ->2π J,τ->(1/(4J)),ωΔ->2π Δ}

iDQamp4//.parameters

Plot[
Evaluate[iDQamp4//.parameters],
{Δ,0,1000},
PlotRange->All,
FrameLabel->{"Δ","iDQ"}
]

compare with simulation (* good *)

H//.parameters

ListPlot[
TransformationAmplitudeTable[opI["z"]->-I opI[1,"+"].opI[2,"+"],
{Pulse[{π/2,"y"}],Delay[τ],Pulse[{π,"x"}],Delay[τ],Pulse[{π/2,"y"}]}//.parameters,
{Δ,0,1000,{71}},
BackgroundGenerator->(H//.parameters)
]
]

contour plot vs Δ and τ

iDQamp4//.{J->50,ωJ->2π J,τ->Jτ /J,ωΔ->2π Δ}

ContourPlot[
Evaluate[
iDQamp4//.{J->50,ωJ->2π J,τ->Jτ /J,ωΔ->2π Δ}
],
{Δ,0,500},{Jτ,0,1},
FrameLabel->{"Δ","J τ"},
LabelStyle->16
]

ContourPlot[
Evaluate[
iDQamp4//.{J->50,ωJ->2π J,τ->Jτ /J,ωΔ->2π Δ}
],
{Δ,0,100},{Jτ,0,10},
FrameLabel->{"Δ","J τ"},
LabelStyle->Medium
]

contour plot in relative units

FullSimplify[PowerExpand[iDQamp4//.{ωJ->2π J,τ->Jτ /J,ωΔ->2π Δ,Δ->ΔbyJ J}]]

θSTtick[θST_]:={Tan[θST],N[θST/°],0.02};
θSTticks=SortBy[Join[{Tan[# °],None}&/@{5,10,30,50},θSTtick/@{20°,40°,60°,70°,80°}],N[First[#]]&];
graphics["INADQ-contour-plot"]=
ContourPlot[
Evaluate[
FullSimplify[PowerExpand[iDQamp4//.{ωJ->2π J,τ->Jτ /J,ωΔ->2π Δ,Δ->ΔbyJ J}]]
],
{ΔbyJ,0,6},{Jτ,0,2},
FrameLabel->{{"J τ","J τ"},{"Δ / J","(Subscript[θ, ST]) / °"}},
GridLines->{Tan[# °]&/@Range[10,85,10],Range[0.25,5,0.25]},
GridLinesStyle->White,
PlotLegends->BarLegend[Automatic,LegendLabel->"DQ amplitude"],
LabelStyle->16,
Contours->Chop@Range[-1,1,0.1],
FrameTicks->{{Automatic,{{0.25,"1/4"},{0.75,"3/4"},{1.25,"5/4"},{1.75,"7/4"}}},{Automatic,θSTticks}},
FrameTicksStyle->{{Darker@Gray,Darker@Gray},{Darker@Gray,14}}
]

SetDirectory[NotebookDirectory[]]
SetDirectory[ParentDirectory[]]
SetDirectory[ParentDirectory[]]
SetDirectory["manuscript/graphics"]

Export["INADQ-contour-plot.pdf",graphics["INADQ-contour-plot"]]

H0 in single-transition operators

Hamiltonian

factor out overall offset

H0=(1/2)ωΔ (opI[1,"z"]- opI[2,"z"])+(1/2)ωΣ (opI[1,"z"]+ opI[2,"z"])+ ωJ opI[1].opI[2]

SetBasis[SingletTripletBasis[]]

MatrixRepresentation[H0]//MatrixForm

MatrixRepresentation[
ωΔ SingleTransitionOperator[{1,3},"x"]
-ωJ SingleTransitionOperator[{1,3},"z"]
-(ωJ/4)SingleTransitionOperator[{1,3},"E"]
+ωΣ SingleTransitionOperator[{2,4},"z"]
+(ωJ/4)SingleTransitionOperator[{2,4},"E"]
]//MatrixForm

H0STops=Exp[I θST CommutationSuperoperator[SingleTransitionOperator[{1,3},"y"]]][
-Ω SingleTransitionOperator[{1,3},"z"]
]+ωΣ SingleTransitionOperator[{2,4},"z"]+(ωJ/4)(SingleTransitionOperator[{2,4},"E"]-SingleTransitionOperator[{1,3},"E"]);

MatrixRepresentation[
H0STops
]//Normal;
%/.SubstitutionRules4//Simplify//MatrixForm

Simplify[Normal@MatrixRepresentation[H0]]==Simplify[Normal@MatrixRepresentation[H0STops]/.SubstitutionRules4]

eigenstates

ket1p=Exp[I θST SingleTransitionOperator[{1,3},"y"]][Ket[1]]
VectorRepresentation[ket1p]//Normal//Simplify//MatrixForm

ket2p=Exp[I θST SingleTransitionOperator[{1,3},"y"]][Ket[2]]
VectorRepresentation[ket2p]//Normal//Simplify//MatrixForm

ket3p=Exp[I θST SingleTransitionOperator[{1,3},"y"]][Ket[3]]
VectorRepresentation[ket3p]//Normal//Simplify//MatrixForm

ket4p=Exp[I θST SingleTransitionOperator[{1,3},"y"]][Ket[4]]
VectorRepresentation[ket4p]//Normal//Simplify//MatrixForm

eigenvalues

eval2

eval2 = ωΣ/2  + ωJ/4

VectorRepresentation[
H0STops [ket2p]-eval2*ket2p
]//Normal//Simplify

eval3

eval3 =Ω/2  - ωJ/4

VectorRepresentation[
H0STops [ket3p]-eval3*ket3p
]//Normal//Simplify

eval4

eval4 = -ωΣ/2  + ωJ/4

VectorRepresentation[
H0STops [ket4p]-eval4*ket4p
]//Normal//Simplify

approx eigenvalues

eval1

eval1

Series[%/.SubstitutionRules4,{ωΔ,0,2}]//PowerExpand//Normal

eval2

eval2

Series[%,{ωΔ,0,2}]//PowerExpand//Normal

eval3

eval3

Series[%/.SubstitutionRules4,{ωΔ,0,2}]//PowerExpand//Normal

eval4

eval4

Series[%/.SubstitutionRules4,{ωΔ,0,2}]//PowerExpand//Normal

density operator transformations

geoDQ

ρ1

ρ1=opI["z"]

MatrixRepresentation[ρ1]==
MatrixRepresentation[2 SingleTransitionOperator[{2,4},"z"]]

ρ2

ρ2=opI["x"]

MatrixRepresentation[ρ2]==
MatrixRepresentation[
2^(1/2)(
SingleTransitionOperator[{2,3},"x"]
+SingleTransitionOperator[{3,4},"x"]
)
]

MatrixRepresentation[ρ2]==
MatrixRepresentation[
2^(-1/2)(
Ket[2].Bra[3]+Ket[3].Bra[2]
+Ket[3].Bra[4]+Ket[4].Bra[3]
)
]

ρ3

Πz13=Exp[-I π SingleTransitionOperator[{1,3},"z"]]

ρ3=Πz13.ρ2.Adjoint[Πz13]

MatrixRepresentation[ρ3]==
MatrixRepresentation[
2^(1/2)(
SingleTransitionOperator[{2,3},"y"]
-SingleTransitionOperator[{3,4},"y"]
)
]

MatrixRepresentation[ρ3]==
MatrixRepresentation[
2^(-1/2)I(
-Ket[2].Bra[3]+Ket[3].Bra[2]
+Ket[3].Bra[4]-Ket[4].Bra[3]
)
]

Commutator[SingleTransitionOperator[{1,3},"z"],SingleTransitionOperator[{2,3},"x"]]
%//MatrixRepresentation//MatrixForm

opA=SingleTransitionOperator[{1,3},"z"];
opB=SingleTransitionOperator[{2,3},"x"];
opC=SingleTransitionOperator[{2,3},"y"];
MatrixRepresentation[Commutator[opA,opB]]==
MatrixRepresentation[(I/2) opC]
MatrixRepresentation[Commutator[opA,opC]]==
MatrixRepresentation[-(I/2) opB]

opA=SingleTransitionOperator[{1,3},"z"];
opB=SingleTransitionOperator[{3,4},"x"];
opC=SingleTransitionOperator[{3,4},"y"];
MatrixRepresentation[Commutator[opA,opB]]==
MatrixRepresentation[-(I/2) opC]
MatrixRepresentation[Commutator[opA,opC]]==
MatrixRepresentation[+(I/2) opB]

ExpressOperator[ρ3, CartesianProductOperatorBasis[]]

ρ4

ρ4=ExpressOperator[opR[{π/2,"y"}].ρ3.opR[{-π/2,"y"}]]

ExpressOperator[ρ4,ShiftAndZOperatorBasis[]]

MatrixRepresentation[ρ4]==
MatrixRepresentation[
I(
Ket[4].Bra[2]-Ket[2].Bra[4]
)
]

MatrixRepresentation[ρ4][[2,4]]

spinorDQ

ρ1

ρ1=opI["z"]

MatrixRepresentation[ρ1]==
MatrixRepresentation[2 SingleTransitionOperator[{2,4},"z"]]

ρ2

ρ2=opI["x"]

MatrixRepresentation[ρ2]==
MatrixRepresentation[
2^(1/2)(
SingleTransitionOperator[{2,3},"x"]
+SingleTransitionOperator[{3,4},"x"]
)
]

MatrixRepresentation[ρ2]==
MatrixRepresentation[
2^(-1/2)(
Ket[2].Bra[3]+Ket[3].Bra[2]
+Ket[3].Bra[4]+Ket[4].Bra[3]
)
]

ρ3

C12=Exp[-I 2 π SingleTransitionOperator[{1,2},"x"]]

ρ3=C12.ρ2.Adjoint[C12]

MatrixRepresentation[ρ3]==
MatrixRepresentation[
2^(1/2)(
-SingleTransitionOperator[{2,3},"x"]
+SingleTransitionOperator[{3,4},"x"]
)
]

MatrixRepresentation[ρ3]==
MatrixRepresentation[
2^(-1/2)(
-Ket[2].Bra[3]-Ket[3].Bra[2]
+Ket[3].Bra[4]+Ket[4].Bra[3]
)
]

ExpressOperator[ρ3, CartesianProductOperatorBasis[]]

ρ4

ρ4=ExpressOperator[opR[{π/2,"x"}].ρ3.opR[{-π/2,"x"}]]

ExpressOperator[ρ4,ShiftAndZOperatorBasis[]]

MatrixRepresentation[ρ4]==
MatrixRepresentation[
I(
Ket[4].Bra[2]-Ket[2].Bra[4]
)
]

PulsePol theory

scaling factors κ

calculations using eq.49 of Sabba 2022

eq.49

κ[1,1,1,1]=2^(1/2)(bigN/(n π))(-1)^((bigN+(n-ν))/(2bigN))Sin[n π/(2bigN)]^2

κ[1,-1,1,-1]=2^(1/2)(bigN/(n π))(-1)^((bigN-(n-ν))/(2bigN))Sin[n π/(2bigN)]^2

κ[1,1,1,1]

κ[1,1,1,1]/.{bigN->4,n->3,ν->1}
%//N

Abs[%]

Arg[κ[1,1,1,1]/.{bigN->4,n->3,ν->1}]

κ[1,-1,1,-1]

κ[1,-1,1,-1]/.{bigN->4,n->3,ν->1}
%//N

Abs[%]

Arg[κ[1,-1,1,-1]/.{bigN->4,n->3,ν->1}]

conclusion

This is surprising: If the equations are correct, the κ scaling factors have phases of π/2 ± π/4. This is not what I expected!

This result also contradicts the table, where the scaling factor is shown as being negative. Maybe something is wrong with eq.49.

new derivation for PulsePol basic element 90y-τ-180x-τ90y

rf Propagator for basic element of duration τR

Urf0[t_]:=Piecewise[
{{opR[{π/2,"y"}],0≤t≤τR/2},
{opR[{π,"x"}].opR[{π/2,"y"}],τR/2≤t≤τR}}
]

Urf0[t]

derive rf Euler

U=opR[{π/2,"y"}];

RoundToSubmultipleOfπ[ang_]:=RoundToSubmultipleOfπ[ang,36]
RoundToSubmultipleOfπ[ang_,n_]:=Round[ang/(π/n)]*(π/n)

RoundToSubmultipleOfπ[3.1]

OperatorToEuler[U_]:=Map[RoundToSubmultipleOfπ,
-Reverse@AxesToEuler@Outer[OperatorAmplitude[U.#1.Adjoint[U]->#2]&,{opI["x"],opI["y"],opI["z"]},{opI["x"],opI["y"],opI["z"]}]]

OperatorToEuler[opR[{π/2,"y"}]]
N[%/°]

OperatorToEuler[opR[{π/2,"y"}]]
N[%/°]

OperatorToEuler[opR[{π,"x"}].opR[{π/2,"y"}]]
N[%/°]

rf Euler

Ωrf0[t_]:=
Piecewise[
{{OperatorToEuler[opR[{π/2,"y"}]],0≤t≤τR/2},
{OperatorToEuler[opR[{π,"x"}].opR[{π/2,"y"}]],τR/2≤t≤τR}}
]

αrf0[t_]:=MapAt[Part[#,1]&,Ωrf0[t],{{1,1,1},{1,2,1}}]

βrf0[t_]:=MapAt[Part[#,2]&,Ωrf0[t],{{1,1,1},{1,2,1}}]

γrf0[t_]:=MapAt[Part[#,3]&,Ωrf0[t],{{1,1,1},{1,2,1}}]

Ωrf0[t]

αrf0[t]

βrf0[t]

γrf0[t]

K factors (eq.37)

K[m_,λ_,μ_]:=τR^(-1)Integrate[
Wignerd[λ,{μ,0}][-βrf0[t]]
Exp[I(μ γrf0[t]+m ωJ t)],
{t,0,τR},Assumptions->{τR∈Reals,τR>0}
]

K[1,1,1]

K[1,1,1]//.{τR->n τJ/bigN,τJ->2π/ωJ}

Simplify[K[1,1,1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ}]

Arg[%]//Simplify

Exp[-I μ (π ν/bigN)]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ,μ->1}

Arg[%]//Simplify

κ factors (eq.36)

Clear[κ]

κ[l_,m_,λ_,μ_]:=Exp[-I μ (π ν/bigN)]K[m,λ,μ]

Exp[-I μ (π ν/bigN)]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ,μ->1}

Arg[%]//Simplify

Simplify[K[1,1,1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ}]
Arg[%]

Simplify[(Exp[-I μ (π ν/bigN)]K[1,1,1])//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ,μ->1}]
Arg[%]

Simplify[κ[1,1,1,1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ,μ->1}]

κ[1,1,1,1]

κ["R431"]=Simplify[κ[1,1,1,1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ}]

Abs[κ["R431"]]//N

Arg[κ["R431"]]//Simplify

Simplify[κ[1,-1,1,-1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ}]

K factor for general symmetry numbers

K[m_,λ_,μ_]:=τR^(-1)Integrate[
Wignerd[λ,{μ,0}][-βrf0[t]]
Exp[I(μ γrf0[t]+m ωJ t)],
{t,0,τR},Assumptions->{τR∈Reals,τR>0}
]

K[1,1,1]

Simplify[ExpToTrig[Exp[I n π/bigN]],Assumptions->{n∈Integers,bigN∈Integers}]

FullSimplify[ExpToTrig[K[1,1,1]//.{τR->n τJ/bigN,τJ->2π/ωJ}],Assumptions->{n∈Integers,bigN∈Integers}]

κ factors for general symmetry numbers

Clear[κ]

κ[l_,m_,λ_,μ_]:=Exp[-I μ (π ν/bigN)]K[m,λ,μ]

Exp[-I μ (π ν/bigN)]//.{τR->n τJ/bigN,τJ->2π/ωJ,μ->1}

Simplify[
ExpToTrig@Expand[
((Exp[-I μ (π ν/bigN)]K[1,1,1])//.{τR->n τJ/bigN,τJ->2π/ωJ,μ->1})
],
Assumptions->{n∈Integers,bigN∈Integers,ν∈Integers}
]

%/.{bigN->4,n->3,ν->1}
Abs[%]
Simplify[%]
%//N

Abs[κ["R431"]]//N

Arg[κ["R431"]]//Simplify

Simplify[κ[1,-1,1,-1]//.{bigN->4,n->3,ν->1,τR->n τJ/bigN,τJ->2π/ωJ}]

Simplify[
ExpToTrig@Expand[
((Exp[-I μ (π ν/bigN)]K[-1,1,-1])//.{τR->n τJ/bigN,τJ->2π/ωJ,μ->-1})
],
Assumptions->{n∈Integers,bigN∈Integers,ν∈Integers}
]

propagation with R431

ρ1

ρ1=opI["z"]

MatrixRepresentation[ρ1]==
MatrixRepresentation[2 SingleTransitionOperator[{2,4},"z"]]

ρ2

ρ2=opI["x"]

MatrixRepresentation[ρ2]==
MatrixRepresentation[
2^(1/2)(
SingleTransitionOperator[{2,3},"x"]
+SingleTransitionOperator[{3,4},"x"]
)
]

MatrixRepresentation[ρ2]==
MatrixRepresentation[
2^(-1/2)(
Ket[2].Bra[3]+Ket[3].Bra[2]
+Ket[3].Bra[4]+Ket[4].Bra[3]
)
]

ρ3

Rx12[β_]:=Exp[-I β SingleTransitionOperator[{1,2},"x"]]

ρ3=Rx12[β].ρ2.Adjoint[Rx12[β]]

ExpressOperator[ρ3,CartesianProductOperatorBasis[]]

Normal@MatrixRepresentation[ρ3]==
Normal@MatrixRepresentation[
2^(1/2)(
(SingleTransitionOperator[{2,3},"x"]Cos[β/2]+SingleTransitionOperator[{1,3},"y"]Sin[β/2])
+SingleTransitionOperator[{3,4},"x"]
)
]

2/(3 0.512)

propagator of one-half of R432

opR[{π,π/4}]
MatrixRepresentation[%]//MatrixForm

opR[{π,-π/4}].opR[{π,π/4}]
MatrixRepresentation[%]//MatrixForm

opR[{π,"z"}]
MatrixRepresentation[%]//MatrixForm

SLIC theory

H0

H0

HSLIC

HSLIC=H0 + ωnut opI["x"]

HTSLIC

T-transform

TransformOp[op_?OperatorQ]:=opR[{-π/2,"y"}].op.opR[{+π/2,"y"}]

HTSLIC

HTSLIC=ExpressOperator[TransformOp[HSLIC]]

HTSLIC//MatrixRepresentation//MatrixForm

decomposition

HA=-(ωΔ/Sqrt[2])SingleTransitionOperator[{1,4},"x"];
HB=ωnut  SingleTransitionOperator[{2,3},"z"];
HC=-(1/4)(ωJ+2ωnut)(SingleTransitionOperator[{1,4},"1"]-SingleTransitionOperator[{2,3},"1"]);

HD=-(ωΣ/Sqrt[2])(SingleTransitionOperator[{2,3},"x"]+SingleTransitionOperator[{3,4},"x"]);
HE=(ωnut-ωJ)SingleTransitionOperator[{1,4},"z"];
HF=(ωΔ/Sqrt[2])SingleTransitionOperator[{1,2},"x"];

(HA+HB+HC+HD+HE+HF)//MatrixRepresentation//MatrixForm

((HA+HB+HC+HD+HE)-HTSLIC)//MatrixRepresentation//MatrixForm

(((HA+HB+HC+HD+HE+HF)-HTSLIC)/.{ωnut->ωJ(1+ϵ)})//MatrixRepresentation//Normal//Expand//MatrixForm

check terms commute

Commutator[HA,HB]//MatrixRepresentation//MatrixForm

Commutator[HB,HC]//MatrixRepresentation//MatrixForm

Commutator[HA,HC]//MatrixRepresentation//MatrixForm

ρ1T

ρ1T=ExpressOperator[TransformOp[ρ1]]

ρ1T//MatrixRepresentation//MatrixForm

UTSLIC

propagators UA, UB, UC

UA=Exp[-I HA τ];
MatrixRepresentation[UA]//MatrixForm

UB=Exp[-I HB τ];
MatrixRepresentation[UB]//MatrixForm

UC=Exp[-I HC τ];
MatrixRepresentation[UC]//MatrixForm

apply UA to ρ1T

UA.ρ1T.Adjoint[UA]
%//MatrixRepresentation//MatrixForm

apply UB to ρ1T

UB.ρ1T.Adjoint[UB]
%//MatrixRepresentation//MatrixForm

apply UC to ρ1T

UC.ρ1T.Adjoint[UC]
%//MatrixRepresentation//MatrixForm

derive ρ2T

introduces κ = 1/sqrt{2}

ρ2T=(UA.UB.UC).ρ1T.Adjoint[(UA.UB.UC)]/.{ωΔ->κ ωΔ Sqrt[2]};
ρ2T//MatrixRepresentation//MatrixForm

OperatorAmplitude[ρ2T->Ket[1].Bra[3]]//Simplify

OperatorAmplitude[ρ2T->Ket[3].Bra[1]]//Simplify

OperatorAmplitude[ρ2T->Ket[2].Bra[3]]//Simplify

OperatorAmplitude[ρ2T->Ket[3].Bra[2]]//Simplify

OperatorAmplitude[ρ2T->Ket[3].Bra[4]]//Simplify

OperatorAmplitude[ρ2T->Ket[4].Bra[3]]//Simplify

OperatorAmplitude[ρ2T->SingleTransitionOperator[{1,3},"x"]]//ExpToTrig//Simplify

OperatorAmplitude[ρ2T->SingleTransitionOperator[{1,3},"y"]]//ExpToTrig//Simplify

OperatorAmplitude[ρ2T->SingleTransitionOperator[{3,4},"x"]]//ExpToTrig//Simplify

OperatorAmplitude[ρ2T->SingleTransitionOperator[{3,4},"y"]]//ExpToTrig//Simplify

OperatorAmplitude[ρ2T->SingleTransitionOperator[{2,3},"x"]]//ExpToTrig//Simplify

OperatorAmplitude[ρ2T->SingleTransitionOperator[{2,3},"y"]]//ExpToTrig//Simplify

MatrixRepresentation[ρ2T]//Normal//FullSimplify//MatrixForm

derive coefficient of important DQ term

ρ2=ExpressOperator[opR[{π/2,"y"}].ρ2T.opR[{-π/2,"y"}]];

OperatorAmplitude[ρ2->opI[1,"+"].opI[2,"+"]]

OperatorAmplitude[ρ2->opI[1,"-"].opI[2,"-"]]

ExpressOperator[opI[1,"+"].opI[2,"+"]-opI[1,"-"].opI[2,"-"]]

OperatorAmplitude[ρ2->opI[1,"x"].opI[2,"y"]+opI[1,"y"].opI[2,"x"]]//FullSimplify

ExpressOperator[TransformOp[opI[1,"x"].opI[2,"y"]+opI[1,"y"].opI[2,"x"]]]

OperatorAmplitude[ρ2T->2opI[1,"x"].opI[2,"z"]+2opI[1,"z"].opI[2,"x"]]//FullSimplify

transformation back to get ρ2

ρ2=ExpressOperator[opR[{π/2,"y"}].ρ2T.opR[{-π/2,"y"}]];

MatrixRepresentation[ρ2]//MatrixForm

MatrixRepresentation[ρ2][[2,4]]

MatrixRepresentation[ρ2][[2,4]]//FullSimplify

MapAt[TrigReduce,,3]

ApplyTo

%//TrigReduce

TrigReduce[Sin[θ]^2]

SLIC theory (v2)

test cases

TestCase[type:_String|{_String..},sa_SparseArray]:=TestCase[type,Normal[sa]]

TestCase["Magnetic Equivalence",expr_]:=FullSimplify[expr/.{ωΔ->0,θe->π/2,ωe->ϵ ωJ}]

TestCase["Exact SLIC",expr_]:=FullSimplify[expr/.{ϵ->0}]

TestCase[{"Magnetic Equivalence","Exact SLIC"},expr_]:=FullSimplify[expr//.{ωΔ->0,θe->π/2,ωe->ϵ ωJ,ϵ->0}]

H0

H0

HSLIC

HSLIC=H0 + ωnut opI["x"]

HTSLIC

T-transform

The sense is reversed compared to the previous version

TransformOp[op_?OperatorQ]:=ExpressOperator[opR[{π/2,"y"}].op.opR[{-π/2,"y"}]]

HTSLIC

HTSLIC=TransformOp[HSLIC]

HTSLIC//MatrixRepresentation//MatrixForm

SLIC match condition -> 12 resonance

HTSLIC/.{ωnut->ωJ}//MatrixRepresentation//MatrixForm

test decomposition

testop=
(
((
-(ωΔ/Sqrt[2])SingleTransitionOperator[{1,2},"x"])+
(+ϵ ωJ SingleTransitionOperator[{1,2},"z"])
)+
(-ωnut SingleTransitionOperator[{3,4},"z"])+
(-(1/4)(ωJ+2ωnut)(SingleTransitionOperator[{1,2},"1"]-SingleTransitionOperator[{3,4},"1"]))+
((ωΣ/Sqrt[2])(SingleTransitionOperator[{2,3},"x"]+SingleTransitionOperator[{3,4},"x"]))+
+
((ωΔ/Sqrt[2])SingleTransitionOperator[{1,4},"x"])
);

testop//MatrixRepresentation//MatrixForm

(HTSLIC-testop)//MatrixRepresentation//MatrixForm

((HTSLIC-testop)/.{ωnut->ωJ(1+ϵ)})//MatrixRepresentation//Normal//Expand//MatrixForm

HA

HA=((
-(ωΔ/Sqrt[2])SingleTransitionOperator[{1,2},"x"])+
(+ϵ ωJ SingleTransitionOperator[{1,2},"z"])
);

HB

HB=(-ωJ(1+ϵ) SingleTransitionOperator[{3,4},"z"]);

HC

HC=(-(1/4)(3ωJ+2ϵ ωJ)(SingleTransitionOperator[{1,2},"1"]-SingleTransitionOperator[{3,4},"1"]))

HD

HD=((ωΣ/Sqrt[2])(SingleTransitionOperator[{2,3},"x"]+SingleTransitionOperator[{3,4},"x"]));

HE

HE=((ωΔ/Sqrt[2])SingleTransitionOperator[{1,4},"x"]);

tests

sum terms

((HA+HB+HC+HD+HE)/.{ωnut->ωJ(1+ϵ)})//MatrixRepresentation//Normal//Expand//MatrixForm

(HA+HB+HC+HD+HE-HTSLIC)/.ωnut->ωJ(1+ϵ)//MatrixRepresentation//Normal//Expand//MatrixForm

check HA, HB, HC commute

Commutator[HA,HB]//MatrixRepresentation//MatrixForm

Commutator[HB,HC]//MatrixRepresentation//MatrixForm

Commutator[HA,HC]//MatrixRepresentation//MatrixForm

check reformulation of HA term

MatrixRepresentation[HA]//MatrixForm

(MatrixRepresentation[
Exp[-I θe SingleTransitionOperator[{1,2},"y"]].
(-ωe SingleTransitionOperator[{1,2},"x"]).
Exp[+I θe SingleTransitionOperator[{1,2},"y"]]
/.{ωe->Sqrt[(ωΔ/Sqrt[2])^2 + (ϵ ωJ)^2],θe->ArcTan[(ωΔ/Sqrt[2]),(ϵ ωJ)]}
]
)//Normal//FullSimplify//MatrixForm

UTSLIC

UA

UA=Exp[-I θe SingleTransitionOperator[{1,2},"y"]].Exp[-I(- ωe τ) SingleTransitionOperator[{1,2},"x"]].Exp[+I θe SingleTransitionOperator[{1,2},"y"]];
UA = Operator[MatrixRepresentation[UA]//Normal//FullSimplify];
MatrixRepresentation[UA]//Normal//FullSimplify//MatrixForm

test cases

TestCase["Magnetic Equivalence",MatrixRepresentation[UA]]//MatrixForm

TestCase["Exact SLIC",MatrixRepresentation[UA]]//MatrixForm

TestCase[{"Magnetic Equivalence","Exact SLIC"},MatrixRepresentation[UA]]//MatrixForm

UB

UB=Exp[-I HB τ]
MatrixRepresentation[UB]//MatrixForm

test cases

TestCase["Magnetic Equivalence",MatrixRepresentation[UB]]//MatrixForm

Normal@MatrixRepresentation[UB]

TestCase["Exact SLIC",Normal@MatrixRepresentation[UB]]//MatrixForm

TestCase[{"Magnetic Equivalence","Exact SLIC"},MatrixRepresentation[UB]]//MatrixForm

UC

UC=Exp[-I HC τ]
MatrixRepresentation[UC]//MatrixForm

test cases

TestCase["Magnetic Equivalence",MatrixRepresentation[UC]]//MatrixForm

TestCase["Exact SLIC",MatrixRepresentation[UC]]//MatrixForm

TestCase[{"Magnetic Equivalence","Exact SLIC"},MatrixRepresentation[UC]]//MatrixForm

UABC

UABC=Operator[Simplify@Normal@MatrixRepresentation[UA.UB.UC]]

UABC//MatrixRepresentation//MatrixForm

test cases

TestCase["Magnetic Equivalence",MatrixRepresentation[UABC]]//MatrixForm

TestCase["Exact SLIC",MatrixRepresentation[UABC]]//MatrixForm

TestCase[{"Magnetic Equivalence","Exact SLIC"},MatrixRepresentation[UABC]]//MatrixForm

ρ1T

ρ1T=ExpressOperator[TransformOp[ρ1]]

ρ1T//MatrixRepresentation//MatrixForm

ρ2T

derive ρ2T

ρ2T=Operator[
FullSimplify@ComplexExpand[
Normal[
MatrixRepresentation[UABC.ρ1T.Adjoint[UABC]]
]
]
]

MatrixRepresentation[ρ2T]//MatrixForm

explore ρ2T cases

TestCase["Magnetic Equivalence",MatrixRepresentation[ρ2T]]//MatrixForm

TestCase["Exact SLIC",MatrixRepresentation[ρ2T]]//MatrixForm

TestCase[{"Magnetic Equivalence","Exact SLIC"},MatrixRepresentation[ρ2T]]//MatrixForm

projections

projections onto pre-DQ

preDQ=2opI[1,"y"].opI[2,"z"]+2opI[1,"z"].opI[2,"y"]

TestCase["Magnetic Equivalence",OperatorAmplitude[ρ2T->preDQ]]//MatrixForm

TestCase["Exact SLIC",OperatorAmplitude[ρ2T->preDQ]]//MatrixForm

TestCase[{"Magnetic Equivalence","Exact SLIC"},OperatorAmplitude[ρ2T->preDQ]]//MatrixForm

project onto ket.bra

OperatorAmplitude[ρ2T->Ket[1].Bra[3]]//Simplify

OperatorAmplitude[ρ2T->Ket[3].Bra[1]]//Simplify

OperatorAmplitude[ρ2T->Ket[2].Bra[3]]//Simplify

OperatorAmplitude[ρ2T->Ket[3].Bra[2]]//Simplify

OperatorAmplitude[ρ2T->Ket[3].Bra[4]]//Simplify

OperatorAmplitude[ρ2T->Ket[4].Bra[3]]//Simplify

project onto Ix13

OperatorAmplitude[ρ2T->SingleTransitionOperator[{1,3},"x"]]//ExpToTrig//Simplify

Simplify[%/.{ϵ->(ωnut-ωJ)/ωJ}]

project onto Iy13

OperatorAmplitude[ρ2T->SingleTransitionOperator[{1,3},"y"]]//ExpToTrig//Simplify

Simplify[%/.{ϵ->(ωnut-ωJ)/ωJ}]

project onto Ix34

OperatorAmplitude[ρ2T->SingleTransitionOperator[{3,4},"x"]]//ExpToTrig//Simplify

Simplify[%/.{ϵ->(ωnut-ωJ)/ωJ}]

project onto Iy34

OperatorAmplitude[ρ2T->SingleTransitionOperator[{3,4},"y"]]//ExpToTrig//Simplify

Simplify[%/.{ϵ->(ωnut-ωJ)/ωJ}]

project onto Ix23

OperatorAmplitude[ρ2T->SingleTransitionOperator[{2,3},"x"]]//ExpToTrig//Simplify

FullSimplify[%/.{ϵ->(ωnut-ωJ)/ωJ}]

FullSimplify[%]

%/.Cos[(1/2) τ (ωe-ωJ-ωnut)]+Cos[(1/2) τ (ωe+ωJ+ωnut)]->2Cos[(1/2) τ ωe]Cos[(1/2) τ(ωJ+ωnut)]

%//Simplify

project onto Iy23

OperatorAmplitude[ρ2T->SingleTransitionOperator[{2,3},"y"]]//ExpToTrig//Simplify

FullSimplify[%/.{ϵ->(ωnut-ωJ)/ωJ}]

transformation back to get ρ2

ρ2=ExpressOperator[opR[{π/2,"y"}].ρ2T.opR[{-π/2,"y"}]];

MatrixRepresentation[ρ2]//MatrixForm

MatrixRepresentation[ρ2][[2,4]]

aDQ

aDQ=FullSimplify@ExpandAll[MatrixRepresentation[ρ2][[2,4]]/.ϵ->(ωnut-ωJ)/ωJ]

test cases

TestCase["Magnetic Equivalence",aDQ]/.ϵ->(ωnut-ωJ)/ωJ

FullSimplify[
TestCase["Exact SLIC",aDQ]/.{ϵ->(ωnut-ωJ)/ωJ,ωnut->ωJ,θe->0,ωe->κ ωΔ}
]

FullSimplify[
TestCase[{"Magnetic Equivalence","Exact SLIC"},aDQ]/.{ϵ->(ωnut-ωJ)/ωJ,ωnut->ωJ,θe->0,ωe->κ ωΔ}
]

SLIC theory (v3)

H0, HJ, Hrf

H0

HJ=ωJ opI[1].opI[2]

Hrf=ωnut opI["x"]

UJ

UJ[t_]:=Exp[-I HJ t]

W0

W0[t_]:=UJ[t].opR[{-ωJ t,"z"}]

MatrixRepresentation[W0[t]]//MatrixForm

Adjoint[W0[t]]

MatrixRepresentation[Adjoint[W0[t]]]//MatrixForm

MatrixRepresentation[W0[t].Adjoint[W0[t]]]//MatrixForm

W (* bugs fixed *)

W[t_]:=opR[{π/2,"y"}].W0[t]

W[t]

MatrixRepresentation[W[t]]//MatrixForm

Adjoint[W[t]]

MatrixRepresentation[Adjoint[W[t]]]//MatrixForm

MatrixRepresentation[W[t]].MatrixRepresentation[Adjoint[W[t]]]
%//MatrixForm

interaction frame transformations

transformation of rf field term

Adjoint[W[0]].(ωnut opI["x"]).W[0]

MatrixRepresentation[%]
MatrixRepresentation[ωnut opI["z"]]
%==%%

transformation of I1z term

Adjoint[W[0]].opI[1,"z"].W[0]

MatrixRepresentation[%]
MatrixRepresentation[-opI[1,"x"]]
%==%%

Adjoint[W[t]].opI[1,"z"].W[t]

MatrixRepresentation[%]
MatrixForm[%]

Dot@@(MatrixRepresentation/@{Adjoint[W[t]],opI[1,"z"],W[t]})
MatrixForm[%]

modified interaction frame transformations to target |1><2| transition

W0a

W0a[t_]:=W0[t]

MatrixRepresentation[W0a[t]]//MatrixForm

Adjoint[W0a[t]]

MatrixRepresentation[Adjoint[W0a[t]]]//MatrixForm

MatrixRepresentation[W0a[t].Adjoint[W0a[t]]]//MatrixForm

SLIC theory (v4)

new round of expressing interaction frame terms:

W0[t_]:=UJ[t].opR[{-ωJ t,"z"}]

W0[t]

HΔtilde

HΔtilde[t_]:=(ωΔ/2)Adjoint[W0[t]].(opI[1,"x"]-opI[2,"x"]).W0[t]

HΔtilde[t]

MatrixRepresentation[HΔtilde[t]]//MatrixForm

OperatorAmplitude[HΔtilde[t]->SingleTransitionOperator[{1,2},"x"]]

OperatorAmplitude[HΔtilde[t]->SingleTransitionOperator[{1,3},"x"]]

OperatorAmplitude[HΔtilde[t]->SingleTransitionOperator[{1,4},"x"]]//ExpToTrig//Simplify

OperatorAmplitude[HΔtilde[t]->SingleTransitionOperator[{1,4},"y"]]//ExpToTrig//Simplify

HΣtilde

HΣtilde[t_]:=(ωΣ/2)Adjoint[W0[t]].(opI[1,"x"]+opI[2,"x"]).W0[t]

HΣtilde[t]

MatrixRepresentation[HΣtilde[t]]//MatrixForm

OperatorAmplitude[HΣtilde[t]->SingleTransitionOperator[{2,3},"x"]]//ExpToTrig//Simplify

OperatorAmplitude[HΣtilde[t]->SingleTransitionOperator[{2,3},"y"]]//ExpToTrig//Simplify

OperatorAmplitude[HΣtilde[t]->SingleTransitionOperator[{3,4},"x"]]//ExpToTrig//Simplify

OperatorAmplitude[HΣtilde[t]->SingleTransitionOperator[{3,4},"y"]]//ExpToTrig//Simplify

ExpressOperator[SingleTransitionOperator[{2,3},"x"]+SingleTransitionOperator[{3,4},"x"]]

W

W[t_]:=opR[{-π/2,"y"}].W0[t]

W[t]

ρSLIC

ρfSLIC[τ_]=W[τ].Exp[-I β SingleTransitionOperator[{1,2},"x"]].opI["x"].Exp[+I β SingleTransitionOperator[{1,2},"x"]].Adjoint[W[τ]]

ExpressOperator[ρfSLIC[τ]]

OperatorAmplitude[ρfSLIC[τ]->SingleTransitionOperator[{2,4},"x"]]

OperatorAmplitude[ρfSLIC[τ]->SingleTransitionOperator[{2,4},"y"]]//ExpToTrig//Simplify

OperatorAmplitude[ρfSLIC[τ]->Ket[2].Bra[4]]//ExpToTrig//Simplify

MatrixRepresentation[ρfSLIC[τ]]//Normal//ExpToTrig//FullSimplify//MatrixForm

MatrixForm[
(MatrixRepresentation[ρfSLIC[τ]]//Normal//ExpToTrig//FullSimplify)/.{τ->(π/2)/ωJ}
]

MatrixForm[
(MatrixRepresentation[ρfSLIC[τ]]//Normal//ExpToTrig//FullSimplify)/.{τ->(3π/2)/ωJ}
]

MatrixForm[
(MatrixRepresentation[ρfSLIC[τ]]//Normal//ExpToTrig//FullSimplify)/.{τ->(5π/2)/ωJ}
]

OperatorAmplitude[ρfstrob->SingleTransitionOperator[{2,4},"y"]]

Density plots

INADEQUATE density plot vs θ and τ, as in the paper

aDQ["INADQ"][θST_,Jτ_]:=
Simplify@PowerExpand@FullSimplify[
(I Cos[θST]Sin[Ω τ]Cos[ωJ τ]-I(Sin[θST]^2 + Cos[θST]^2 Cos[Ω τ])Sin[ωJ τ])//.{
Ω->Sqrt[ωJ^2 + ωΔ^2],
ωJ->2π J,
τ->Jτ/ J,
ωΔ->ωJ Tan[θST]
}
]

aDQ["INADQ"][θST,Jτ]

aDQF["INADQ"][θST_,Jτ_]:=Abs[aDQ["INADQ"][θST,Jτ]]^2

aDQF["INADQ"][θST,Jτ]

DensityPlot[
Evaluate[PowerExpand@FullSimplify[aDQF["INADQ"][θST,Jτ]/.θST->N[θSTdeg °]]],
{θSTdeg,0,90},{Jτ,0,2},
FrameLabel->{{"J τ",None},{"(Subscript[θ, ST]) / °",None}},
LabelStyle->16,
PlotLegends->Automatic,
PlotPoints->100,
PlotLabel->"INADQ"
]

DensityPlot[
Evaluate[PowerExpand@FullSimplify[aDQF["INADQ"][θST,Jτ]/.θST->N[θSTdeg °]]],
{θSTdeg,0,90},{Jτ,0,40},
FrameLabel->{{"J τ",None},{"(Subscript[θ, ST]) / °",None}},
LabelStyle->16,
PlotLegends->Automatic,
PlotPoints->100,
PlotLabel->"INADQ"
]

DensityPlot[
Evaluate[PowerExpand@FullSimplify[aDQF["INADQ"][θST,Jτ]/.θST->N[θSTdeg °]]],
{θSTdeg,0,90},{Jτ,0,40},
FrameLabel->{{"J τ",None},{"(Subscript[θ, ST]) / °",None}},
LabelStyle->16,
PlotLegends->Automatic,
PlotPoints->100,
PlotLabel->"INADQ",
ScalingFunctions->{None,"Log"},
FrameTicks->{{{0.25,1,5,10,20},None},Automatic}
]

SLIC density plot vs θ and τ (approximate version)

aDQ["SLIC"][θST_,Jτ_]:=I Sin[ωJ τ]Sin[β12/4]//.{
β12->-κ ωΔ τ, κ->2^(-1/2),
ωJ->2π J,
τ->Jτ/ J,
ωΔ->ωJ Tan[θST]
};

aDQ["SLIC"][θST,Jτ]

aDQF["SLIC"][θST_,Jτ_]:=Chop[Abs[aDQ["SLIC"][θST,Jτ]]^2]

PowerExpand@FullSimplify[aDQF["SLIC"][θST,Jτ]/.θST->N[θSTdeg °]]

DensityPlot[
Evaluate[PowerExpand@FullSimplify[aDQF["SLIC"][θST,Jτ]/.θST->N[θSTdeg °]]],
{θSTdeg,0,90},{Jτ,0,40},
FrameLabel->{{"J τ",None},{"(Subscript[θ, ST]) / °",None}},
LabelStyle->16,
PlotLegends->Automatic,
PlotPoints->500,
PlotLabel->"SLIC (approximate)",
ScalingFunctions->{None,"Log"},
FrameTicks->{{{0.25,1,5,10,20},None},Automatic}
]

PulsePol density plot vs θ and τ (av. Ham. approx)

aDQ["PulsePol"][θST_,Jτ_]:=-I Sin[ωnut12 τ/4]^2//.{
ωnut12->κ ωΔ ,
 κ->2(1+Sqrt[2])/(3π),
ωJ->2π J,
τ->Jτ/ J,
ωΔ->ωJ Tan[θST]
};

aDQ["PulsePol"][θST,Jτ]

aDQF["PulsePol"][θST_,Jτ_]:=Chop[Abs[aDQ["PulsePol"][θST,Jτ]]^2]

PowerExpand@FullSimplify[aDQF["PulsePol"][θST,Jτ]/.θST->N[θSTdeg °]]

DensityPlot[
Evaluate[PowerExpand@FullSimplify[aDQF["PulsePol"][θST,Jτ]/.θST->N[θSTdeg °]]],
{θSTdeg,0,90},{Jτ,0,40},
FrameLabel->{{"J τ",None},{"(Subscript[θ, ST]) / °",None}},
LabelStyle->16,
PlotLegends->Automatic,
PlotPoints->100,
PlotLabel->"PulsePol (using av Ham)",
ScalingFunctions->{None,"Log"},
FrameTicks->{{{0.25,1,5,10,20,40},None},Automatic}
]

comparison with numerical calculations

SetBasis[SingletTripletBasis[]]

aDQ["INADQ","simulated"][θST_,Jτ_]:=
Module[{τ,ωΔ,J},(
(* take J=1 *)
J=1;
τ=Jτ/J;
ωΔ=2π J Tan[θST];
TransformationAmplitude[
opI["z"]->Ket[2].Bra[4],
{Pulse[{π/2,"x"}],
Delay[τ],
Pulse[{π,"y"}],
Delay[τ],
Pulse[{π/2,"x"}]
},
BackgroundGenerator->2π J opI[1].opI[2]+ωΔ (1/2)(opI[1,"z"]-opI[2,"z"])
]
)];

aDQ["INADQ","simulated"][π/2-(π/1000),1/4]

ListDensityPlot[
Evaluate[
Flatten[#,1]&@
Table[{θSTdeg,Jτ,Abs[aDQ["INADQ","simulated"][N[θSTdeg °],Jτ]]^2},
{θSTdeg,0,90,2},{Jτ,0,2,0.02}
]
],
FrameLabel->{{"J τ",None},{"(Subscript[θ, ST]) / °",None}},
LabelStyle->16,
PlotLegends->Automatic,
PlotLabel->"INADQ"
]

INADQ/PulsePol comparison

ResourceFunction["PolygonMarker"][All]

star=ResourceFunction["PolygonMarker"]["FivePointedStar"];

exptcoords["INADQ"]={4.3,2 325.0 (10^(-3))255.0}

ColorData[{"GrayTones","Reverse"}]

plot["INADQ"]=
Show[{
DensityPlot[
Evaluate[PowerExpand@FullSimplify[aDQF["INADQ"][θST,JT/2]/.θST->N[θSTdeg °]]],
{θSTdeg,0,30},{JT,1,200},
FrameLabel->{{"J T",None},{"(Subscript[θ, ST]) / °",None}},
LabelStyle->16,
PlotLegends->Automatic,
PlotPoints->100,
PlotLabel->"INADQ",
ScalingFunctions->{None,"Log"},
FrameTicks->{{{0.25,1,2,5,10,20,40,100,200},None},Automatic},
ColorFunction->ColorData[{"GrayTones","Reverse"}]
],
ListPlot[{exptcoords["INADQ"]},
PlotMarkers->{Graphics[{Yellow,star}],Large},
ScalingFunctions->{None,"Log"}]
}]

star=ResourceFunction["PolygonMarker"]["FivePointedStar"];

exptcoords["PulsePol"]={4.3,9 8 1.405 (10^(-3))255.0}

plot["PulsePol"]=
Show[{
DensityPlot[
Evaluate[PowerExpand@FullSimplify[aDQF["PulsePol"][θST,JT]/.θST->N[θSTdeg °]]],
{θSTdeg,0,30},{JT,1,200},
FrameLabel->{{"J T",None},{"(Subscript[θ, ST]) / °",None}},
LabelStyle->16,
PlotLegends->Automatic,
PlotPoints->100,
PlotLabel->"PulsePol",
ScalingFunctions->{None,"Log"},
FrameTicks->{{{0.25,1,2,5,10,20,40,100,200},None},Automatic},
ColorFunction->ColorData[{"GrayTones","Reverse"}]
],
ListPlot[{exptcoords["PulsePol"]},
PlotMarkers->{Graphics[{Yellow,star}],Large},
ScalingFunctions->{None,"Log"}]
}]

PulsePol transformation

Needs["SpinDynamica`"]
SetSpinSystem[2];
SetOperatorBasis[CartesianProductOperatorBasis[]]

ρ3 = 2^(1/2)(
SingleTransitionOperator[{2,3},"x"]Cos[β/2]+
SingleTransitionOperator[{1,3},"y"]Sin[β/2]+
SingleTransitionOperator[{3,4},"x"]
)

MatrixRepresentation[ρ3]//MatrixForm

ρ4=RotationSuperoperator[{π/2,"x"}][ρ3]

MatrixRepresentation[ρ4]//MatrixForm

OperatorAmplitude[ρ4->SingleTransitionOperator[{2,4},"y"]]

OperatorAmplitude[ρ4->Ket[2].Bra[4]]

cSLIC avHam theory (from Mohamed)

setup

Needs["SpinDynamica`"]

SetSpinSystem[2]

SetBasis[SingletTripletBasis[]]

SetOperatorBasis[CartesianProductOperatorBasis[]]

Hamiltonian terms

HJ = ωJ opI[1].opI[2]

HΔ = ωΔ SingleTransitionOperator[{1,3},"x"]

piecewise definition of rf field for cSLIC

Solve[ωstrong/(ωstrong+ωJ) ==α,ωstrong]

ϕ[t_]=Piecewise[{{ωa t,0<=t<=τa},{-ωb (t-τa)+ωa τa,τa<t<=τa+τb},{ωa (t-τa-τb)-ωb τb +ωa τa,τa+τb<t<=2τa+τb}}]

cSLIC interaction frame Hamiltonian

assumptions={τb<τa,ωJ∈Reals,τa∈PositiveReals,τb∈PositiveReals,t∈PositiveReals,ωa∈Reals,ωb∈Reals,ωstrong∈Reals,t∈PositiveReals,ϵrf∈Reals,α∈PositiveReals&&q∈PositiveReals}

useful pre-simplification

Wϵ[t_]=FullSimplify[((MatrixRepresentation@RotationOperator[{-π/2,π/2}].MatrixRepresentation@Exp[-I HJ t].MatrixRepresentation@RotationOperator[{-ϕ[t],"z"}])//Normal),assumptions];

interaction frame Hamiltonian

integrand=FullSimplify[Adjoint[Wϵ[t]].MatrixRepresentation[HΔ].Wϵ[t],assumptions];

cSLIC average Hamiltonian

(splitting up the integral makes it faster.)

integral=(1/(2τa+τb))(Integrate[integrand,{t,0,τa},Assumptions->{τa>0&&τb>0}]+Integrate[integrand,{t,τa,τa+τb},Assumptions->{τa>0&&τb>0}]+Integrate[integrand,{t,τa+τb,2τa+τb},Assumptions->{τa>0&&τb>0}]);

The following is essentially the exact average Hamiltonian.

this substitution makes life easier. we want to do away with ωa and ωb and define everything relative to ωJ.

Solve[ωstrong/(ωstrong+ωJ)==α,ωstrong]

HavExact=integral/.τa->α π/ωJ/.τb->2π(1- α)/ωJ/.ωb->(α/(1-α))ωa/.ωa->(1+ϵrf)ωJ//FullSimplify

the average Hamiltonian strictly contains (Subscript[I, x]^12) and (Subscript[I, x]^14) terms.

let’s rewrite the Hamiltonian as a sum of (Subscript[I, x]^12) and (Subscript[I, x]^14) operators, each multiplied by some coefficient. These two coefficients can then be called c12 and c14

OperatorAmplitude[Operator@simplerintegral->{SingleTransitionOperator[{1,2},"x"],SingleTransitionOperator[{1,4},"x"]}]

c12 and c14 are then associated with the resonant components f12 and f14:

f12==π ϵrf;f14==π (2+ϵrf);

such that c12 and c14 may be written as:

{c12,c14}=((-ωΔ)/Sqrt[2]){(α (1+f12/π)  Sinc[ α f12])/(1+α f12/π),-((α(1-f14/π))/(1-α f14/π))Sinc[ α f14]}//Simplify

sanity check that this result is true:

EqualTo[(({c12,c14}{SingleTransitionOperator[{1,2},"x"],SingleTransitionOperator[{1,4},"x"]}//Total//MatrixRepresentation//Normal)/.Sinc[x_]->Sin[x]/x/.f12->π ϵrf/.f14->π(2+ϵrf))][HavExact]//FullSimplify

in the limit of an infinitely short compensating pulse, the coefficients are incredibly simple:

Limit[{c12,c14},α->1]//FullSimplify

expression for DQ excitation amplitude under a general effective Hamiltonian of the form “c12 Ix12 - c14 Ix14”

we have the following basic rotation operator identity, which makes life easier:

Needs["MaTeX`"]

MaTeX["R_y^{24} (2\\theta_{\\epsilon})^{\\dagger} R_{x}^{24} ( \\omega_{\\epsilon} t) R_y^{24} (2\\theta_{\\epsilon}) = \\exp\\left(-i \\omega_{\\epsilon} \\left\\{\\cos(\\theta_{\\epsilon}) I_{x}^{12} - \\sin(\\theta_{\\epsilon}) I_{x}^{14}\\right\\} t\\right)",FontSize->20]

sanity check:

EqualTo[Exp[+I 2θϵ SingleTransitionOperator[{2,4},"y"]].Exp[-I ωϵ SingleTransitionOperator[{1,2},"x"] t].Exp[-I 2θϵ SingleTransitionOperator[{2,4},"y"]]//MatrixRepresentation//Normal][Exp[-I ωϵ (Cos[θϵ]SingleTransitionOperator[{1,2},"x"]-Sin[θϵ]SingleTransitionOperator[{1,4},"x"])t]//MatrixRepresentation//Normal]//FullSimplify

this is essentially the DQ excitation amplitude.

Clear[θϵ,ωϵ]

aDQ=OperatorAmplitude[-2opI[1,"x"].opI[2,"z"]-2opI[1,"z"].opI[2,"x"]->(Exp[+I 2θϵ SingleTransitionOperator[{2,4},"y"]].Exp[-I ωϵ SingleTransitionOperator[{1,2},"x"] t].Exp[-I 2θϵ SingleTransitionOperator[{2,4},"y"]]).opI["x"].Adjoint@(Exp[+I 2θϵ SingleTransitionOperator[{2,4},"y"]].Exp[-I ωϵ SingleTransitionOperator[{1,2},"x"] t].Exp[-I 2θϵ SingleTransitionOperator[{2,4},"y"]])]

(or, we can add imaginary factors just so Malcolm is happy:)

Bra[2].(RotationSuperoperator[{π/2,0}][(Exp[+I 2θϵ SingleTransitionOperator[{2,4},"y"]].Exp[-I ωϵ SingleTransitionOperator[{1,2},"x"] t].Exp[-I 2θϵ SingleTransitionOperator[{2,4},"y"]]).opI["x"].Adjoint@(Exp[+I 2θϵ SingleTransitionOperator[{2,4},"y"]].Exp[-I ωϵ SingleTransitionOperator[{1,2},"x"] t].Exp[-I 2θϵ SingleTransitionOperator[{2,4},"y"]])]).Ket[4]//FullSimplify

analytical forms for the case of cSLIC

let’s apply this expression to our case. we simply need to extract the parameters θϵ and ωϵ  and we’re good.

this is the detuning angle in general:

θϵexact=ArcTan[-c14/c12/.f14->f12+2π/.f12->π ϵrf]//FullSimplify

and in the delta pulse limit:

θϵ=Limit[θϵexact,α->1]/.Sinc[x_]->Sin[x]/x//FullSimplify//FullSimplify

this is the effective field magnitude in general:

ωϵexact=FullSimplify[Sqrt[(c14^2)+(c12^2)]/.f12->π ϵrf/.f14->π(2+ϵrf),ωΔ∈PositiveReals&&α∈PositiveReals]//FullSimplify

and in the delta pulse limit:

ωϵ=Limit[ωϵexact,α->1]//FullSimplify

basic trig identities show that these forms are equivalent:

Plot[{Abs@Sinc[π ϵrf]Sec[θϵ],√(((Sinc[π ϵrf])^2)+((Sinc[π (2+ϵrf)])^2))},{ϵrf,-8,8}]

Plot[{Sinc[π ϵrf]Sec[ArcTan[2+ϵrf,ϵrf]]Cos[ArcTan[2+ϵrf,ϵrf]],Sinc[π ϵrf]Sec[ArcTan[2+ϵrf,ϵrf]]Sin[ArcTan[2+ϵrf,ϵrf]],Cos[θϵ]Sinc[π ϵrf]Sec[θϵ],Sin[θϵ]Sinc[π ϵrf]Sec[θϵ]},{ϵrf,-8,8},PlotRange->All]

this is our expression in the delta pulse limit:

Cos[2 θϵ] ((Sin[(t ωϵ)/4])^2)

Plot[aDQ/.t->Sqrt[2] 2 π/ωΔ,{ϵrf,-3,1},PlotRange->All,FrameStyle->Directive[32,Thickness[0.0015],Black,FontFamily->"Helvetica Light"],ImageSize->600,FrameLabel->{"(Subscript[ϵ, rf])","i×(Subscript[a, DQ])"}]

this is what happens outside the delta pulse limit:

(((Cos[2 θϵexact] ((Sin[(t ωϵexact)/4])^2)))^2)

DensityPlot[(((Cos[2 θϵexact] ((Sin[(t ωϵexact)/4])^2)))^2)/.t->Sqrt[2] 2 π/ωΔ,{α,0.5,1},{ϵrf,-2-1,0+1},PlotPoints->{100,100},FrameLabel->{"α","(Subscript[ϵ, rf])"},PlotRange->All,FrameStyle->Directive[32,Thickness[0.0015],Black,FontFamily->"Helvetica Light"],ColorFunctionScaling->False,ColorFunction->(Blend[{Black,Purple,Darker@Darker@Blue,Darker@Blue,Blue,Lighter@Blue,Lighter@Lighter@Blue,Cyan,White},#]&),PlotLegends->True,ImageSize->500]

DensityPlot[(((Cos[2 θϵexact] ((Sin[(t ωϵexact)/4])^2)))^2)/.t->Sqrt[2] 2 π/ωΔ,{α,0,1},{ϵrf,-2-10,0+10},PlotPoints->{100,100},FrameLabel->{"α","(Subscript[ϵ, rf])"},PlotRange->All,ColorFunctionScaling->False,ColorFunction->(Blend[{Black,Purple,Darker@Darker@Blue,Darker@Blue,Blue,Lighter@Blue,Cyan,White},#]&)]