Magnus expansion and new analytic approximations applied to the NMR Example in the Schrödinger picture.

Álgebra of the Pauli sigma matrices

Cmt[x_,x_]:=0;
Cmt[a___, 0,     b___] := 0;
Cmt[a___, x_+y_, b___] :=    Cmt[a,x,b]+Cmt[a,y,b];
Cmt[a___, t^n_ x_,  b___] := t^n Cmt[a,x,b];
Cmt[a___, n_ x_Cmt, b___] := n   Cmt[a,x,b];
Cmt[a___, x_  y_,  b___] :=  Cmt[a, Expand[x y], b];
Cmt[a___, n_ sig[j_Integer], b___] := n *
     Cmt[a,sig[j],b];
     
 Cmt[sig[i_],sig[j_]]:= - Cmt[sig[j],sig[i]] /; i > j   
 
 Cmt[z_Cmt,sig[i_]]:= -Cmt[sig[i],z];
Cmt[Cmt[sig[i_],sig[j_]],z_Cmt]:=
  -Cmt[Cmt[sig[j],z],sig[i]] - 
   Cmt[Cmt[z,sig[i]],sig[j]];
Cmt[z_Cmt,Cmt[sig[i_],sig[j_]]]:=
  -Cmt[sig[i],Cmt[sig[j],z]] -
   Cmt[sig[j],Cmt[z,sig[i]]]; 
   
Cmt[sig[1],sig[2]]:= 2 I sig[3];
Cmt[sig[1],sig[3]]:= - 2 I sig[2];
Cmt[sig[2],sig[3]]:= 2 I sig[1];

Coefficient matrix (Hamiltonian computed in the interaction picture)

A[t_]:= -(I om0/2)sig[3] -I b (Cos[om t] sig[1] + Sin[om t] sig[2]);

Magnus expansion

First term

omg[1][t_] := Simplify[Integrate[A[t], t]];
ome[1][t_]:= Collect[omg[1][t] - (omg[1][t] /. t -> 0),{sig[1],sig[2],sig[3]}];

ome[1][t]

tper=2 Pi/om;

Recurrence

ss[n_Integer,1][t_]:= ss[n,1][t] = Simplify[
   Cmt[ome[n-1][t],A[t]]];
ss[n_Integer,j_Integer][t_]:= ss[n,j][t]=Sum[
   Simplify[Cmt[ome[m][t],ss[n-m,j-1][t]]], 
   {m,1,n-j}];
omg[n_Integer][t_]:= omg[n][t]=Simplify[Sum[ 
   (BernoulliB[j]/j!) * (
  Integrate[ss[n,j][t],t]), {j,1,n-1}]];
ome[n_Integer][t_]:= ome[n][t] = Collect[Simplify[omg[n][t] - 
    (omg[n][t] /. t->0)],{sig[1],sig[2],sig[3]}];

Magnus approximants 2 to 8

ome[2][t]

ome[3][t]

Timing[ome[4][t];]

Timing[ome[5][t];]

Timing[ome[6][t];]

Timing[ome[7][t];]

Timing[ome[8][t];]

Truncated Magnus series in terms of Pauli matrices

omex[j_]:=Collect[Simplify[Sum[ome[n][t],{n,1,j}]],{sig[1],sig[2],sig[3]}];

omex[1]

Do[omex[j]=Collect[Simplify[omex[j-1]+ome[j][t]],{sig[1],sig[2],sig[3]}],{j=2,9,1}]

ww[i_]:=ww[i]=Table[Simplify[Coefficient[(omex[i]/I),sig[j]]],{j,1,3}];

Table[ww[i],{i,1,8}];

Order 4 (in t)

ww[2]

Order 6 (in t)

ww[4]

Order 8 (in t) (large output)

ww[6];

Order 10 (in t) (large output)

ww[8];

Exact transition probability

ompri=Sqrt[4 b^2 + (om-om0)^2];

probex=(2 b Sin[ompri*t/2]/ompri)^2;

fundib=probex/.{b->0.4,om->4,om0->1};

Plot[fundib,{t,0,2 π}]

Transition probability computed with Magnus and difference with the exact result

diff[j_]:=Module[{vn,promag,vp1},vp1=ww[j]^2;vn=Sqrt[Total[vp1]];promag=(Sin[vn]/vn)^2 * Total[Take[vp1,2]];
Abs[(probex-promag)]/.{b->0.4,om->4,om0->1}
];

diferencias=Table[diff[j],{j,2,8,2}];

LogPlot[diferencias,{t,0,2Pi}]

New analytic approximations in terms of univariate integrals

Univariate integrals

B[i_]:=Expand[(1/t^(i)) Integrate[(u-t/2)^i  A[u],{u,0,t}]];

Order 4

co1=Collect[Simplify[Expand[Cmt[B[0],B[1]]]],{sig[1],sig[2],sig[3]}];

m4apro=Collect[Simplify[B[0]-co1],{sig[1],sig[2],sig[3]}];

m4apro

vm4=Table[Simplify[Coefficient[(m4apro/I),sig[j]]],{j,1,3}]

We compare with the Magnus approximant

Normal[Series[Expand[ww[2]],{t,0,4}]]==Normal[Series[Expand[vm4],{t,0,4}]]

dmagun4:=Module[{vn,prom,vp1},vp1=vm4^2;vn=Sqrt[Total[vp1]];prom=(Sin[vn]/vn)^2 * Total[Take[vp1,2]];
Abs[(probex-prom)]/.{b->0.4,om->4,om0->1}
];

LogPlot[{diff[2],dmagun4},{t,0,2Pi}]

Order 6

alfa1=(9/4)B[0]-15B[2]; alfa2=12 B[1]; alfa3=-15 B[0]+180 B[2];

c1=Collect[Simplify[Expand[Cmt[alfa1,alfa2]]],{sig[1],sig[2],sig[3]}];

c2=-(1/60)Collect[Simplify[Expand[Cmt[alfa1,2alfa3+c1]]],{sig[1],sig[2],sig[3]}];

c3=Collect[Simplify[Expand[Cmt[-20alfa1-alfa3+c1,alfa2+c2]]],{sig[1],sig[2],sig[3]}];

om6apro=Collect[Simplify[alfa1+(1/12)alfa3+(1/240)c3],{sig[1],sig[2],sig[3]}];

vm6=Table[Simplify[Coefficient[(om6apro/I),sig[j]]],{j,1,3}];

Normal[Series[Expand[ww[4]],{t,0,6}]]==Normal[Series[Expand[vm6],{t,0,6}]]

dmagun6:=Module[{vn,prom,vp1},vp1=vm6^2;vn=Sqrt[Total[vp1]];prom=(Sin[vn]/vn)^2 * Total[Take[vp1,2]];
Abs[(probex-prom)]/.{b->0.4,om->4,om0->1}
];

LogPlot[{diff[4],dmagun6},{t,0,2Pi}]

LogPlot[{diferencias[[1]],diferencias[[2]],diferencias[[4]],dmagun4,dmagun6},{t,0,2Pi},PlotStyle->{Thick,Thick,Thick,Dashed,Dashed},Frame->True,FrameLabel->{t,error},PlotRange->{{0,6.5},{10^(-10),0}}]