Coherent control and numerical time-dependent perturbation theory

Initialization and examples

Initialization of needed packages

Off[General::spell]

Introduce NCM

Unprotect[NonCommutativeMultiply];
NonCommutativeMultiply[A_?MatrixQ,B_?MatrixQ]:=KroneckerProduct[A,B];
$RecursionLimit=Infinity;

define Pauli matrices

Sx={{0,1},{1,0}};Sz={{1,0},{0,-1}};Sy={{0,-I},{I,0}};S0={{1,0},{0,1}};
sx=Sx/2;sy=Sy/2;sz=Sz/2;      (*spin-1/2 operators*)
sp=sx+I sy;sm=sx-I sy;

We can now easily construct direct product of matrices

sx**sz**S0//MatrixForm

define pulse shapes

R1=2 Pi (a0+a1 Cos[2 Pi t]+a2 Cos[4 Pi t]+a3 Cos[6 Pi t]+a4 Cos[8 Pi t]+a5 Cos[10 Pi t]); (* generic Fourier expansion *)

0th order shapes

G010[t_Real]:=Evaluate[Sqrt[Pi]/(1-a/2) (1-a t^2/tau^2)(Exp[-t^2/tau^2]/tau)/.{tau->0.10,a->0,t->t-0.5-Floor[t]}];
K1[t_Real]:=Evaluate[R1/.{a0->0.5,a1->-0.5,a2->0.,a3->0.,a4->0.,a5->0.}];

Plot[{G010[t],K1[t]},{t,0,1}]

NIntegrate[K1[t]/Pi,{t,0,1}]

1st order shapes

H010[t_Real]:=Evaluate[Sqrt[Pi]/(1-a/2) (1-a t^2/tau^2)(Exp[-t^2/tau^2]/tau)/.{tau->0.10,a->9609317217/10000000000,t->t-0.5-Floor[t]}];
S1[t_Real]:=Evaluate[R1/.{a0->0.5,a1->-1.2053194466,a2->0.4796460175,a3->0.225673429,a4->0,a5->0}];
S2[t_Real]:=Evaluate[R1/.{a0->0.5,a1->-1.1950755990,a2->0.7841246569,a3->0.0738054432,a4->-0.1628545011,a5->0}];

Plot[{H010[t],S1[t],S2[t]},{t,0,1}]

NIntegrate[S1[t],{t,0,1}]

eqs={D[phi[t],t]==S1[t],D[s[t],t]==Sin[phi[t]],s[0]==0,phi[0]==0};

sol=NDSolve[eqs,{s[t],phi[t]},{t,0,1}]

(s[t]/.sol[[1]])/.t->1

2nd order shapes

Q1[t_Real]:=Evaluate[R1/.{a0->0.5,a1->-1.1374003264,a2->1.5774784244,a3->-0.6825954606,a4->-0.2574826374,a5->0}];
Q2[t_Real]:=Evaluate[R1/.{a0->0.5,a1->-1.0965122417,a2->1.5309957409,a3->-1.1470791601,a4->0.0020722004,a5->0.2105234605}];

Plot[{Q1[t],Q2[t]},{t,0,1},PlotRange->All]

We can now plot

Plot[{H010[t],G010[t]},{t,0,1},PlotRange->All]

Define Single-qubit refocusing sequences

SEQ0s={{0},{0}};(*no refocusing*)
SEQ2a={{A,-A},{0,0}};
SEQ4p={{A,0,A,0},{0,-A,0,A}};

SEQ8p={{A,0,A,0,-A,0,-A,0},{0,-A,0,A,0,A,0,-A}};

SEQ16p={{A,0,A,0,-A,0,-A,0,-A,0,-A,0,A,0,A,0},{0,-A,0,A,0,A,0,-A,0,-A,0,A,0,A,0,-A}};

SEQ8s={{A,0,A,0,0,A,0,A},{0,-A,0,A,A,0,-A,0}};

SEQ8a={{A,0,A,0,0,-A,0,-A},{0,-A,0,A,-A,0,A,0}};

This function defines a field  as a sequence of pulses

DoField[tt_Real,X_]:=Block[{i,tmod,num},
num=Length[Flatten[X]];
tmod=tt-num Floor[tt/num];
(Take[Flatten[X][[1+IntegerPart[tmod]]]])];

We can now define

Vx[t_Real]=DoField[t,SEQ4p[[1]]/.A->G010[t]];

SEQ4p[[1]]

Plot[Vx[t],{t,0,8},PlotRange->All]

SEQ4p[[2]]

Vy[t_Real]=DoField[t,SEQ4p[[2]]/.A->G010[t]];

Plot[{Vx[t],Vy[t]},{t,0,4},PlotRange->All]

Solve evolution equations

this constructs a system of ODEs with RHS x,for variables Psi with initial conditions Psi0

DoEq[x_,Psi_,Psi0_]:=Block[{i,y={}},
For[i=1,i≤Length[x],i++,y=Append[y,(x[[i]]==0)]];
For[i=1,i≤Length[Psi],i++,y=Append[y,(Psi[[i]]==Psi0[[i]])]];
y];

DoEq[{D[phi[t],t]-G010[t],D[s[t],t]-Sin[phi[t]]},{s[0],phi[0]},{0,0}]

NDSolve[%,{s[t],phi[t]},{t,0,1}]

Plot[Evaluate[{s[t],phi[t]}/.%56],{t,0,1}]

ψ={y1[t],y2[t]};DoEq[D[ψ,t]-sz . ψ,ψ/.t->0,{1,0}]

plot evolution samples

DoPlots[R_, (* density matrix as a fun of t *)
observables_, (* list of pairs {operator, name} *)
tint_ (* pair t0,tmax *)
]:=Block[{obs,opt,name,val,tmin,tmax,val0,val1,plotlist={}},
{tmin, tmax}=tint;
For[obs=1,obs≤Length[observables],obs++,{opt,name}=observables[[obs]];val=((Tr[R.opt])/(Tr[R]));val0=ToString[Chop[val/.t->tmin]];val1=ToString[Chop[val/.t->tmax]];PlotFuncts=Plot[Evaluate[{Re[val],Im[val]}],{t,tmin,tmax},
PlotRange->Automatic,
(* PlotStyle->{RGBColor[0,1,0],RGBColor[0,0,1]}, *) 
AxesLabel->{"t",name},
PlotLabel->Style[name<>": "<>val0<>" to "<>val1,FontWeight->"Bold"]];
AppendTo[plotlist,PlotFuncts];];
(* Print[GraphicsRow[plotlist]]; *)
plotlist]

Full Unitary matrix

Table[ToExpression["u"<>ToString[i]<>"[t]"],{i,1,4}]

U=Partition[Table[ToExpression["u"<>ToString[i]<>"[t]"],{i,1,4}],2]//MatrixForm

U=.

WFEvol[Rinit_:{},(*initial DM for plots only*)
Time_,(*duration*)
Htot_,(*full Hamiltonian*)
observables_:{} (*list of pairs {obs,name} if plots are required*)]:=Block[{N,U,Eq,EqSol,R},
N=Length[Htot];(*size of the hilbert space*)
U=Partition[Table[ToExpression["u"<>ToString[i]<>"[t]"],{i,1,N^2}],N];
Eq=DoEq[Flatten[D[U,t]+I Htot.U],Flatten[U]/.t->0,Flatten[IdentityMatrix[N]]];
EqSol=NDSolve[Evaluate[Eq],Evaluate[Flatten[U]],{t,0,Time},MaxSteps->Infinity,AccuracyGoal->10];
If[Length[observables]>0,
R=(U.Rinit.ConjugateTranspose[U])/.EqSol[[1]];
plots1=DoPlots[R,observables,{0,Time}]];
U/.EqSol[[1]]/.t->Time];(*return final evolution matrix*)

sz

Plot[{Vx[t],Vy[t]},{t,0,4},PlotRange->All]

ans=WFEvol[(S0+Sy)/2,4,0.5sz+sx Vx[t]+sy Vy[t],{{sx,"(Subscript[s, x])"},{sy,"(Subscript[s, y])"},{sz,"(Subscript[s, z])"}}];
Print[GraphicsRow[plots1]]

SEQ2a

Vx[t_Real]=DoField[t,SEQ2a[[1]]/.A->K1[t]];
Vy[t_Real]=DoField[t,SEQ2a[[2]]/.A->K1[t]];

Plot[{Vx[t],Vy[t]},{t,0,2},PlotRange->All]

ans=WFEvol[(S0+Sx)/2,2,0.5sz+sx Vx[t]+sy Vy[t],{{sx,"(Subscript[s, x])"},{sy,"(Subscript[s, y])"}}];Print[GraphicsRow[plots1]]

Perturbative expansion for evolution

WFPert[TT_,H0_,H1_,PrintMatr_:1,Nord_:6,eps_:(1/(10^3))]:=Block[{N,U,Eq0,EqSol0,Eq,EqSol,AnsList,HH1,S,S0,Sf},N=Length[H0];U=Partition[Table[ToExpression["u"<>ToString[i]<>"[t]"],{i,1,(N^2)}],N];Eq0=DoEq[Flatten[(∂_t)U+I H0.U],Flatten[U]/.t->0,Flatten[IdentityMatrix[N]]];EqSol0=NDSolve[Evaluate[Eq0],Evaluate[Flatten[U]],{t,0,TT},MaxSteps->∞,AccuracyGoal->10];HH1=ConjugateTranspose[U].H1.U/.EqSol0[[1]];
S0=IdentityMatrix[N];
AnsList={Sf=Chop[U/.EqSol0[[1]]/.t->TT,(eps/100)]};If[PrintMatr==2,Print["U0[T]=",Sf,"  Norm[DU0]=",Norm[Flatten[Sf-IdentityMatrix[N]]]],If[PrintMatr==1,Print["Norm[DU0]=",Chop[Norm[Flatten[Sf-IdentityMatrix[N]]]]]]];S=Partition[Table[ToExpression["s"<>ToString[i]<>"[t]"],{i,1,(N^2)}],N];For[i=1,i≤Nord,i++,
Eq=DoEq[Flatten[(∂_t)S+I HH1.S0],Flatten[S]/.t->0,Flatten[ConstantArray[0,{N,N}]]];EqSol=NDSolve[Evaluate[Eq],Evaluate[Flatten[S]],{t,0,TT},MaxSteps->∞,AccuracyGoal->10];
S0=S/.EqSol[[1]];AppendTo[AnsList,Sf=Chop[S0/.t->TT,(eps/100)]];If[PrintMatr==2,Print["U",i,"[T]=",Sf,"  Norm[U",i,"[T]]=",Norm[Flatten[Sf]]],If[PrintMatr==1,Print["Norm[U",i,"[T]]=",Chop[Norm[Flatten[Sf]]]]]];];If[PrintMatr≠-1,Print["Order=",DoOrder[AnsList,eps]]];AnsList];

Compute the sequence order

DoOrder[AnsLis_,(*perturbative expansion of U*)
eps_:10^-3,(*accuracy*)
cut_:{1,1} (*cut this fraction of size-- use {nn,nn} when working with an oscillator*)]:=Block[{i,N,N1,N2,matr,order=-1},
N=Length[AnsLis[[1]]];
{N1,N2}=cut;
N1=N/N1;N2=N/N2;
matr=Take[AnsLis[[1]]-IdentityMatrix[N],N1,N2];
If[Norm[Flatten[matr]]<eps,order++];
For[i=2,i≤Length[AnsLis],i++,matr=Take[AnsLis[[i]],N1,N2];
If[order==i-2&&Norm[Flatten[matr]]<eps,order++];];;
If[Length[AnsLis]>order+1&&Norm[Flatten[matr=Take[AnsLis[[order+2]],N1,N2]]]<100 eps,
order=ToString[order]<>"*",order=ToString[order]];
order];

SEQ2a

Vx[t_Real]=DoField[t,SEQ2a[[1]]/.A->G010[t]];
Vy[t_Real]=DoField[t,SEQ2a[[2]]/.A->G010[t]];

ans=WFEvol[{},2,0.5sz+sx Vx[t]+sy Vy[t]]

ans2=WFPert[2,sx Vx[t]+sy Vy[t],0.5sz,1,10];

ans2

Chop[ans-Apply[Plus,ans2],10^-5]

Try same pulse sequence with different pulse shapes

shapelist={"G010","S1","Q1"};
seq=SEQ2a;
len=Length[seq[[1]]];
For[shapeno=1,shapeno≤Length[shapelist],shapeno++,Vx[t_Real]=DoField[t,seq[[1]]/.A->ToExpression[shapelist[[shapeno]]<>"[t]"]];Vy[t_Real]=DoField[t,seq[[2]]/.A->ToExpression[shapelist[[shapeno]]<>"[t]"]];Print[Style["Working with "<>shapelist[[shapeno]],FontSize->17,FontWeight->"Bold"]];
WFPert[len,sx Vx[t]+sy Vy[t],1.8` sz,1];
res=WFEvol[(1/2) {{1,1},{1,1}},len,1.8` sz+sx Vx[t]+sy Vy[t]];
Print[res];
NN=2;
Print["1-F=",(NN^2-Abs[((Tr[res])^2)])/(NN^2+NN)];];

seq={{A,0,-A,0,0,-A,0,A},{0,A,0,A,A,0,A,0}};

seq={{A,0,-A,0},{0,A,0,-A}};

shapelist={"G010","S1","Q1"};
len=Length[seq[[1]]];
For[shapeno=1,shapeno≤Length[shapelist],shapeno++,V1[t_Real]=DoField[t,seq[[1]]/.A->ToExpression[shapelist[[shapeno]]<>"[t]"]];V2[t_Real]=DoField[t,seq[[2]]/.A->ToExpression[shapelist[[shapeno]]<>"[t]"]];Print[Style["Working with "<>shapelist[[shapeno]],FontSize->17,FontWeight->"Bold"]];
Print[Plot[{V1[t],V2[t]},{t,0,len},PlotRange->All]];
H0=sx**S0 V1[t]+S0**sy V2[t];H1=0.5 sz**sz+0.4 sx**sx+0.4 sy**sy+0.2 sz**S0+0.1 S0**sz;
WFPert[len,H0,H1,1];
res=WFEvol[{},len,H0+ H1];
NN=Length[H0];
Print["1-F=",(NN^2-Abs[((Tr[res])^2)])/(NN^2+NN)];];

shapelist={"G010","S1","Q1"};
len=Length[seq[[1]]];
For[shapeno=1,shapeno≤Length[shapelist],shapeno++,V1[t_Real]=DoField[t,seq[[1]]/.A->ToExpression[shapelist[[shapeno]]<>"[t]"]];V2[t_Real]=DoField[t,seq[[2]]/.A->ToExpression[shapelist[[shapeno]]<>"[t]"]];Print[Style["Working with "<>shapelist[[shapeno]],FontSize->17,FontWeight->"Bold"]];
Print[Plot[{V1[t],V2[t]},{t,0,len},PlotRange->All]];
H0=(sx**S0**S0+S0**S0**sx) V1[t]+S0**sy**S0 V2[t];
H1=0.5 (sz**sz**S0+S0**sz**sz)+0.4( sx**sx**S0+S0**sx**sx)+0.4( sy**sy**S0+S0**sy**sy)+0.2 sz**S0**S0+0.1 S0**sz**S0-0.25 S0**S0**sz;
WFPert[len,H0,H1,1,3];
res=WFEvol[{},len,H0+ 2 H1];
NN=Length[H0];
Print["1-F=",(NN^2-Abs[((Tr[res])^2)])/(NN^2+NN)];];