Lecture 5 - Coherent control and pulse design

by Leonid Pryadko

Goal: make given evolution matrix by the end of the interval τ

Initialization (numerical)

Definition of Pauli matrices

(Subscript[σ, x])={{0,1},{1,0}};(Subscript[σ, y])={{0,-I},{I,0}};(Subscript[σ, z])={{1,0},{0,-1}};(Subscript[σ, 0])={{1,0},{0,1}};

WF will be a two-component spinor

ψ={y1[t],y2[t]};

Equation generator from Lecture 3

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]

Simplest case: decoupling

Example : constant external Hamiltonian with single spin

H1= (Subscript[Δ, x])(Subscript[σ, x])+ (Subscript[Δ, y])(Subscript[σ, y])+ (Subscript[Δ, z])(Subscript[σ, z]);

U1=Simplify[MatrixExp[-I H1 t],Assumptions->(Subscript[Δ, x])>0 &&(Subscript[Δ, y])>0&& (Subscript[Δ, z])>0 && t>0];

U1

Upx=-I (Subscript[σ, x]);Umx=-Upx;Upy=-I (Subscript[σ, y]);Umy=-Upx;

Sequence {X, Y, -X, Y}

U=(U1.Upy.U1).(U1.Umx.U1).(U1.Upy.U1).(U1.Upx.U1);

U=Series[U,{t,0,2}];Simplify[U,Assumptions->(Subscript[Δ, x])>0 &&(Subscript[Δ, y])>0&& (Subscript[Δ, z])>0 && t>0]//MatrixForm

Sequence {X, Y, -X, Y, Y, -X, Y, X}

Ua=Simplify[U1.Upx.U1.U1.Upy.U1.U1.Umx.U1.U1.Upy.U1];

Ub=Simplify[U1.Upy.U1.U1.Umx.U1.U1.Upy.U1.U1.Upx.U1];

U=Ua.Ub;

U=Series[U,{t,0,4}];Simplify[U,Assumptions->(Subscript[Δ, x])>0 &&(Subscript[Δ, y])>0&& (Subscript[Δ, z])>0 && t>0]//MatrixForm

Sequence symmetry plays a special role.  Can we try doing the same with finite width pulses?

Time - dependent perturbation theory.

n=3;

Table[ToExpression["u" <> ToString[i] <> "[t]"], {i,1,n^2}]

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

Flatten[U]

n=.

WFPert[TT_,  (* duration *) 
       H0_,    (* main Hamiltonian *)       
       H1_,    (* pert Hamiltonian *)
       Nord_:10, (* max order *)
       eps_:10^-3 (* how to round off the errors *)] :=     
    Block[{n, U, Eq0, EqSol0, Eq, EqSol, AnsList, HH1, S, S0, Sf}, 
       n=Length[H0];(* dimension of the Hamiltonian *)
U = Partition[Table[ToExpression["u" <> ToString[i] <> "[t]"], {i,1,n^2}],n];
       Eq0 = DoEq[Flatten[D[U, t] + I H0 . U], Flatten[U] /. t -> 0,Flatten[IdentityMatrix[n]]]; 
       EqSol0 = NDSolve[Evaluate[Eq0], Evaluate[Flatten[U]], {t, 0.,TT},AccuracyGoal->10]; 
       HH1=(ConjugateTranspose[U]. H1. U)/.EqSol0[[1]];
       S0=IdentityMatrix[n];
       AnsList={Sf=Chop[(U/.EqSol0[[1]])/.t->TT,eps/100]}; (* final value of zeroth order S *)        S = Partition[Table[ToExpression["s" <> ToString[i] <> "[t]"], {i,1,n^2}],n];
       For[i=1,i<=Nord,i++,
Eq = DoEq[Flatten[D[S, t] + I HH1 . S0], Flatten[S] /. t -> 0,Flatten[ConstantArray[0,{n,n}]]]; 
           EqSol=NDSolve[Evaluate[Eq], Evaluate[Flatten[S]], {t, 0., TT},AccuracyGoal->10]; 
           S0=S/.EqSol[[1]];       
           AppendTo[AnsList,Sf=Chop[(S0/.t->TT),eps/100]];
];
       AnsList];

Here how it works for a Gaussian pulse:

V01[t_Real]:=Evaluate[Sqrt[Pi] 
   (Exp[-t^2/tau^2]/tau)/.{tau->0.05,t->t-0.5-Floor[t]}];

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

1-NIntegrate[V01[t],{t,0,1}]/Pi

H0=.

H0[t_]:=0.5 (Subscript[σ, x]) V01[t];H1=(* 0.45 (Subscript[σ, x])+0.3(Subscript[σ, y])+ *)0.7(Subscript[σ, z]);

ans=WFPert[1.0,H0[t],H1,4];For[i=1,i≤Length[ans],i++,Print[i,"  ",MatrixForm[ans[[i]]]]];

Error appears already in the first (linear) order.  Here is the numerical way of optimizing this -- we need to kill the integrals∫Sin[α[t]]d t, ∫Cos[α[t]]d t over the pulse duration.

Specify a symmetric pulse shape for the interval 0<t<2 π:

t=.;A0=.;A1=.;A2=.;

Omega[x_]:=A0+A1 Cos[x-Pi]+A2 Cos[2 x-2 Pi];

This is the corresponding integral

Alpha[x_]:=A0 x+A1 (Sin[x-Pi])+A2 Sin[2x-2 Pi]/2;

Plot[Evaluate[{Omega[x],Alpha[x]}/.{A0->0.5,A1->1.2, A2->0}],{x,0,2Pi},PlotRange->All]

Solve by direct integration

This is a slow way for computing the various integrals

Fun[xx_]:=Block[{A0,A1,A2},
{A0,A1,A2}=xx;
{NIntegrate[Evaluate[Cos[Alpha[x]]],{x,0, 2Pi},{AccuracyGoal->14,MaxRecursion->20}],NIntegrate[Evaluate[Sin[Alpha[x]]],{x,0, 2Pi},{AccuracyGoal->14,MaxRecursion->20}](*,NIntegrate[NIntegrate[Cos[Alpha[t]-Alpha[x]],{t,0,x},{AccuracyGoal->14,MaxRecursion->20}],{x,0, 2Pi},{AccuracyGoal->14,MaxRecursion->20}],NIntegrate[NIntegrate[Sin[Alpha[t]-Alpha[x]],{t,0,x},{AccuracyGoal->14,MaxRecursion->20}],{x,0, 2 Pi},{AccuracyGoal->14,MaxRecursion->20}]*)}
]

Fun[{0.5,1.2,0}]

This function computes the  rms deviation

Gr[c0_Real,c1_Real,c2_Real]:=Block[{res,val},
res=Fun[{c0,c1,c2}];
 val=Norm[res];
Print[{c0,c1,c2},"  -> ",res," -> ",val]; 
val]

Gr[0.5,1.17,0.]

QQ[a1_,a2_]:=FindMinimum[Gr[0.5,A1,A2],{{A1,a1},{A2,a2}},
MaxIterations->5,AccuracyGoal->5];

ans=QQ[1.2,0.]

Here is the obtained function

Plot[(Omega[x]/.ans[[2]])/.A0->0.5,{x,0,2 Pi}]

NIntegrate[(Omega[x]/.ans[[2]])/.A0->0.5,{x,0,2 Pi}]

(Omega[x]/.ans[[2]])/.A0->0.5

Try using this with perturbation series

H0[t_]:=0.5(Subscript[σ, x])((Omega[t]/.ans[[2]])/.A0->0.5)

WFPert[2. Pi,H0[t],H1,2]

Now, error is in the second order!

Try solving the differential equation instead (faster)

t=.;A0=.;A1=.;A2=.;

Fun1[A0_Real,A1_Real,A2_Real]:=Block[{t},
sol=NDSolve[{
D[A[t],t]==A0+A1 Cos[t-Pi]+A2 Cos[2 t-2 Pi],A[0]==0,
D[c1[t],t]==Cos[A[t]],c1[0]==0,
D[s1[t],t]==Sin[A[t]],s1[0]==0},
{A,c1,s1},{t,0,2 Pi}];
{c1[2 Pi],s1[2 Pi]}/.sol[[1]]]

Fun1[0.5,1.2,0.0]

QQ1[a1_Real,a2_Real]:=Norm[Fun1[0.5,a1,a2]]

QQ1[1.2,0.0]

ans1=NMinimize[{QQ1[A1,A2]+0.5 A2^2,Norm[{A1,A2}]<5},{A1,A2}]

Plot[0.5+A1 Cos[t-Pi]+A2 Cos[2 t-2 Pi]/.ans1[[2]],{t,0,2 Pi}]

H0[t_]:=0.5(Subscript[σ, x])(0.5+A1 Cos[t-Pi]+A2 Cos[2 t-2 Pi]/.ans1[[2]])

WFPert[2. Pi,H0[t],H1,2]

ans2=FindMinimum[Evaluate[QQ1[A1,A2]],{{A1,1.2},{A2,0}}]

Plot[0.5+A1 Cos[t-Pi]+A2 Cos[2 t-2 Pi]/.ans2[[2]],{t,0,2 Pi}]

H0[t_]:=0.5(Subscript[σ, x])(0.5+A1 Cos[t-Pi]+A2 Cos[2 t-2 Pi]/.ans2[[2]])

WFPert[2. Pi,H0[t],H1,2]

Clearly, the errors in second order differ - we can try to do better than that!

Same to 2nd order

Fun2[A1_Real,A2_Real, A3_Real]:=Block[{A0=0.5,t},
sol=NDSolve[{
D[A[t],t]==A0+A1 Cos[t]+A2 Cos[2 t]+A3 Cos[3 t],A[0]==0,
D[c1[t],t]==Cos[A[t]],c1[0]==0,
D[s1[t],t]==Sin[A[t]],s1[0]==0,
D[c2[t],t]==Sin[A[t]] s1[t]+Cos[A[t]] c1[t],c2[0]==0,
D[s2[t],t]==Sin[A[t]] c1[t]-Cos[A[t]] s1[t],s2[0]==0},
{A,c1,s1,c2,s2},{t,0,2 Pi}];
{c1[2 Pi],s1[2 Pi],c2[2 Pi],s2[2 Pi]}/.sol[[1]]]

Fun2[-1.137,1.57,0.68]

NMinimize[{Norm[Fun2[b1,b2,b3]],Norm[{b1,b2,b3}]<5},{b1,b2,b3}]

ans3=FindMinimum[Norm[Fun2[b1,b2,b3]],{{b1,-1.2},{b2,0},{b3,0}}]

Plot[0.5+b1 Cos[t]+b2 Cos[2 t]+b3 Cos[3 t]/.ans3[[2]],{t,0, 2 Pi}]

H0[t_]:=0.5(Subscript[σ, x])(0.5+b1 Cos[t]+b2 Cos[2 t]+b3 Cos[3 t]/.ans3[[2]])

WFPert[2. Pi,H0[t],0.5H1,3]

Additional lobes suppress errors to quadratic order!

Here are some better shapes which also vanish at the ends of the interval:

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]);

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}];
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}];

First-order self-refocusing pulse

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

Second-order self-refocusing pulse

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

These pulses are more carefully designed to have more of derivatives to vanish

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

How it works for a sequence

Here is the function to generate a sequence of pulses into a continuous function

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

seq4={{A,0,-A,0},{0,A,0,A}};
seq8={{A,0,-A,0,0,-A,0,A},{0,A,0,A,A,0,A,0}};

H1=0.3 (Subscript[σ, x])+.5 (Subscript[σ, y])+H1

Try 4 - pulse sequence with 1 st - order pulses

SEQ=seq4;
Vx[t_]:=DoField[t,SEQ[[1]]/.A->S1[t]];
Vy[t_]:=DoField[t,SEQ[[2]]/.A->S1[t]];
H0[t_]:=0.5 (Subscript[σ, x])Vx[t]+0.5 (Subscript[σ, y])Vy[t];

Plot[{Vx[t],Vy[t]},{t,0,4},PlotStyle->{Hue[0.3],Hue[0.8]}]

WFPert[4.,H0[t],H1,3]

Try 8 - pulse sequence with 1 st - order pulses

SEQ=seq8;
Vx[t_]:=DoField[t,SEQ[[1]]/.A->S1[t]];
Vy[t_]:=DoField[t,SEQ[[2]]/.A->S1[t]];
H0[t_]:=0.5 (Subscript[σ, x])Vx[t]+0.5 (Subscript[σ, y])Vy[t];

Plot[{Vx[t],Vy[t]},{t,0,8},PlotStyle->{Hue[0.3],Hue[0.8]}]

WFPert[8.,H0[t],H1,3]

Try 4 - pulse sequence with 2nd - order pulses

SEQ=seq4;
Vx[t_]:=DoField[t,SEQ[[1]]/.A->Q1[t]];
Vy[t_]:=DoField[t,SEQ[[2]]/.A->Q1[t]];
H0[t_]:=0.5 (Subscript[σ, x])Vx[t]+0.5 (Subscript[σ, y])Vy[t];

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

WFPert[4.,H0[t],H1,4]

Try 8 - pulse sequence with 2nd - order pulses

SEQ=seq8;
Vx[t_]:=DoField[t,SEQ[[1]]/.A->Q1[t]];
Vy[t_]:=DoField[t,SEQ[[2]]/.A->Q1[t]];
H0[t_]:=0.5 (Subscript[σ, x])Vx[t]+0.5 (Subscript[σ, y])Vy[t];

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

WFPert[8.,H0[t],H1,4]

H1

Initialization (analytical)

Define Dot and Outer products for multi-spin operators

(Subscript[σ, x])=.;(Subscript[σ, y])=.;(Subscript[σ, z])=.;(Subscript[σ, 0])=.;

Off[General::spell]

Unprotect[Dot,NonCommutativeMultiply];

SpinQ[x_]:=(Head[x]==S);

SpinQ[S[t]]

NumbQ[x_]:=FreeQ[x,S[_]]&&FreeQ[x,A[_]];

S[x]**S[0]**S[0]//FullForm

Dot[NonCommutativeMultiply[a__],NonCommutativeMultiply[b__]]:=
Block[{aa,bb,ans,i,n},
aa={a};
bb={b};
n=Length[aa];
If[Length[bb]≠n,Return[False]];
For[i=1;ans={},i≤n,i++,
AppendTo[ans,aa[[i]].bb[[i]]];
];
Apply[NonCommutativeMultiply,ans]];

(S[x]**S[0]**S[0]).(S[y]**S[z]**S[0])

S[x].S[y]:=I S[z];S[y].S[z]:=I S[x];S[z].S[x]:=I S[y];
S[y].S[x]:=-I S[z];S[z].S[y]:=-I S[x];S[x].S[z]:=-I S[y];
S[x_].S[x_]:=S[0];S[0].S[x_]:=S[x];S[x_].S[0]:=S[x];
S[0].A[x_]:=A[x];A[x_].S[0]:=A[x];

Dot[A_?NumbQ x_,y_]:=A (x.y);
Dot[x_,A_?NumbQ y_]:=A (x.y);
Dot[x_,z_+y_]:=(x.y)+(x.z);
Dot[x_+y_,z_]:=(x.z)+(y.z);
Dot[x_,A_?NumbQ]:=A x;
Dot[A_?NumbQ,x_]:=A x;

NonCommutativeMultiply[x_,z_+y_]:=(x**y)+(x**z);
NonCommutativeMultiply[x_+y_,z_]:=(x**z)+(y**z);

NonCommutativeMultiply[A_?NumbQ x_,y_]:=A (x**y);
NonCommutativeMultiply[x_,A_?NumbQ y_]:=A (x**y);

NonCommutativeMultiply[x_,A_?NumbQ]:=A x;
NonCommutativeMultiply[A_?NumbQ,x_]:=A x;

Format[S[x_]]:=(Subscript[σ, x]);
Format[S[x_],TeXForm]:="{\\sigma_{x}}";
Format[A[x_]]:=(Subscript[A, x]);
Format[A[x_],TeXForm]:="{{\\rm A}_{x}}";
Format[S[0]]:=1;Format[S[0],TeXForm]:="\\openone";

(S[x]**S[0]**A[0]).(S[y]**S[z]**A[z])

Examples

We can now do analytical calculations with multiple spins

(S[x]**S[z]).(S[y]**S[0])

Or even with unknown operators

H=S[x]**A[x]+S[y]**A[y]+S[z]**A[z]+S[0]**A[0]

Here is the square of this operator

H.H

Evolution operator to quadratic order

sub1x={S[y]->c1 S[y]+s1 S[z],S[z]->c1 S[z]-s1 S[y],S[x]->x1 S[x],A[0]->x1 A[0]};
sub1y=sub1x/.{x->y,y->z,z->x};
sub1z=sub1y/.{x->y,y->z,z->x};
sub2x=sub1x/.{c1->c2,s1->s2,x1->x2};
sub2y=sub1y/.{c1->c2,s1->s2,x1->x2};
sub2z=sub1z/.{c1->c2,s1->s2,x1->x2};

These substitutions account for pulse symmetry

sub2={c1 c2->0,s1 s2->s^2/2,c1 s2->-α/2,s1 c2->α/2,
         c1 x2->-ζ,x1 c2->ζ,s1 x2->s/2,x1 s2->s/2,x1 x2->1/2};

sub1={c1->0,s1->s,x1->1};

Ux=S[0]**S[0]-I T(H/.sub1x)-T^2(H/.sub1x).(H/.sub2x);
Ux=Expand[(-I S[x]**S[0].(Expand[Ux]/.sub2))/.sub1];

Collect[Ux,T]

Uy=S[0]**S[0]-I T(H/.sub1y)-T^2(H/.sub1y).(H/.sub2y);
Uy=Expand[(-I S[y]**S[0].(Expand[Uy]/.sub2))/.sub1];

Uz =S[0]**S[0]-I T(H/.sub1z)-T^2(H/.sub1z).(H/.sub2z);
Uz=Expand[(-I S[z]**S[0].(Expand[Uz]/.sub2))/.sub1];

Vx=Collect[-Ux/.{s->-s,α->-α},T];
Vy=Collect[-Uy/.{s->-s,α->-α},T];
Vz=Collect[-Uz/.{s->-s,α->-α},T];

DoProd[vars__]:=Block[{a,n,z1,z2,UU,UUz},
a={vars};
n=Length[a];
If[n==1,Return[a[[1]]]];
a=Reverse[a];
z1=CoefficientList[a[[1]]+T^3,T];
z2=CoefficientList[a[[2]]+T^3,T];
UU={(z2[[1]]).(z1[[1]]), (z2[[2]]).(z1[[1]])+(z2[[1]]).(z1[[2]]),
z2[[3]].z1[[1]]+z2[[2]].z1[[2]]+z2[[1]].z1[[3]],
z2[[4]].z1[[1]]+z2[[3]].z1[[2]]+z2[[2]].z1[[3]]+z2[[1]].z1[[4]]};
For[i=3,i≤n,i++,
z2=CoefficientList[a[[i]]+T^3,T];
UUz={z2[[1]].UU[[1]], z2[[2]].UU[[1]]+z2[[1]].UU[[2]], z2[[3]].UU[[1]]+z2[[2]].UU[[2]]+z2[[1]].UU[[3]],
z2[[4]].UU[[1]]+z2[[3]].UU[[2]]+z2[[2]].UU[[3]]+z2[[1]].UU[[4]]};
UU=UUz;
];
Collect[UU[[1]]+T UU[[2]]+T^2 UU[[3]],{T,NonCommutativeMultiply[__]}]
]

$RecursionLimit = Infinity;

Expand[Ux/.T^2->0]

Effective Hamiltonian

tt=Collect[I (S[x]**S[0]).Ux/.{A[x]->0,A[y]->0,A[0]->0},{T,s,α,ζ}]-S[0]**S[0]

Collect[ Expand[-I (tt-DoProd[tt,tt]/2+DoProd[tt,tt,tt]/3)],{T,NonCommutativeMultiply[__]}]

U4=Expand[DoProd[Ux,Uy,Vx,Uy]/.{s->0,α->0}]

Effective Hamiltonian

tt=Collect[U4,T]-S[0]**S[0];
Collect[ Expand[I (tt-DoProd[tt,tt]/2)/(4 T)],
{T,NonCommutativeMultiply[__]}]

U8=Expand[DoProd[Ux,Uy,Vx,Uy,Uy,Vx,Uy,Ux]/.{s->0,α->0}]

tt=Collect[U8,T]-S[0]**S[0];
Collect[ Expand[I (tt-DoProd[tt,tt]/2)],{T,NonCommutativeMultiply[__]}]/(8 T)

To quadratic order, this gives the effective Hamiltonian is the external Hamiltonian.