Second Order in T Phase Shift Calculation

Lagrangian and Equations of Motion

r={x[t],y[t],z[t]};

L=m*(((D[r,t].D[r,t])/2)-g*z[t])

Euler-Lagrange Equations

eqnsEL={D[D[L,x'[t]],t]==D[L,x[t]],D[D[L,y'[t]],t]==D[L,y[t]],D[D[L,z'[t]],t]==D[L,z[t]]}//FullSimplify

Time series solution to equation of motion

TOrder=2;
xTSeries=Series[Sum[xC[j]*(t^j),{j,0,TOrder}],{t,0,TOrder}]
yTSeries=Series[Sum[yC[j]*(t^j),{j,0,TOrder}],{t,0,TOrder}]
zTSeries=Series[Sum[zC[j]*(t^j),{j,0,TOrder}],{t,0,TOrder}]

eqnsELTSeries=eqnsEL/.{x->xTSeries,y->yTSeries,z->zTSeries}/.{Derivative[q_][f_][t]:>D[f,{t,q}],f_[t]->f}

initial conditions specify zeroth and first order coefficients, use Euler-Lagrange equations to solve for remaining coefficients to fourth order in time t.  Initial value of x at time 0 denoted x1, initial value of velocity in the x direction at time 0 denotes vx1, etc...

coefs=(Solve[eqnsELTSeries,Join[xC/@Range[2,TOrder],yC/@Range[2,TOrder],zC/@Range[2,TOrder]]]/.{xC[0]->x1,xC[1]->vx1,yC[0]->y1,yC[1]->vy1,zC[0]->z1,zC[1]->vz1})[[1]]

xSol=xTSeries/.coefs/.{xC[0]->x1,xC[1]->vx1,yC[0]->y1,yC[1]->vy1,zC[0]->z1,zC[1]->vz1}//Normal//Simplify
ySol=yTSeries/.coefs/.{xC[0]->x1,xC[1]->vx1,yC[0]->y1,yC[1]->vy1,zC[0]->z1,zC[1]->vz1}//Normal//Simplify
zSol=zTSeries/.coefs/.{xC[0]->x1,xC[1]->vx1,yC[0]->y1,yC[1]->vy1,zC[0]->z1,zC[1]->vz1}//Normal//Simplify

Folded Trajectories

specify times that kicks occur and what the momentum kicks received by the upper and lower arms are for each pulse.  Assume n*ℏ*k momentum transfer.  Keep terms only up to 2nd order in time, using a small parameter ϵ to keep track of how many powers of time.

timeList={0,T,2*T};

kickListUp={{0,0,n*ℏ*k},{0,0,-n*ℏ*k},{0,0,n*ℏ*k}};

kickListDown={{0,0,0},{0,0,n*ℏ*k},{0,0,0}};

trajUp=FoldList[{(xSol/.{x1->(#1[[1]]/.t->time),vx1->((D[#1[[1]],t]/.t->time)+(kickX/m)),y1->(#1[[2]]/.t->time),vy1->((D[#1[[2]],t]/.t->time)+(kickY/m)),z1->(#1[[3]]/.t->time),vz1->((D[#1[[3]],t]/.t->time)+(kickZ/m))}),(ySol/.{x1->(#1[[1]]/.t->time),vx1->((D[#1[[1]],t]/.t->time)+(kickX/m)),y1->(#1[[2]]/.t->time),vy1->((D[#1[[2]],t]/.t->time)+(kickY/m)),z1->(#1[[3]]/.t->time),vz1->((D[#1[[3]],t]/.t->time)+(kickZ/m))}),(zSol/.{x1->(#1[[1]]/.t->time),vx1->((D[#1[[1]],t]/.t->time)+(kickX/m)),y1->(#1[[2]]/.t->time),vy1->((D[#1[[2]],t]/.t->time)+(kickY/m)),z1->(#1[[3]]/.t->time),vz1->((D[#1[[3]],t]/.t->time)+(kickZ/m))})}/.{t->(t-time)}/.#2&,{x0+v0x*t,y0+v0y*t,z0+v0z*t},{{kickX->kickListUp[[1,1]],kickY->kickListUp[[1,2]],kickZ->kickListUp[[1,3]],time->timeList[[1]]},{kickX->kickListUp[[2,1]],kickY->kickListUp[[2,2]],kickZ->kickListUp[[2,3]],time->timeList[[2]]},{kickX->kickListUp[[3,1]],kickY->kickListUp[[3,2]],kickZ->kickListUp[[3,3]],time->timeList[[3]]}}];

trajUpSeries=Series[trajUp/.{t->ϵ*t,T->ϵ*T},{ϵ,0,TOrder}]//Simplify

trajDown=FoldList[{(xSol/.{x1->(#1[[1]]/.t->time),vx1->((D[#1[[1]],t]/.t->time)+(kickX/m)),y1->(#1[[2]]/.t->time),vy1->((D[#1[[2]],t]/.t->time)+(kickY/m)),z1->(#1[[3]]/.t->time),vz1->((D[#1[[3]],t]/.t->time)+(kickZ/m))}),(ySol/.{x1->(#1[[1]]/.t->time),vx1->((D[#1[[1]],t]/.t->time)+(kickX/m)),y1->(#1[[2]]/.t->time),vy1->((D[#1[[2]],t]/.t->time)+(kickY/m)),z1->(#1[[3]]/.t->time),vz1->((D[#1[[3]],t]/.t->time)+(kickZ/m))}),(zSol/.{x1->(#1[[1]]/.t->time),vx1->((D[#1[[1]],t]/.t->time)+(kickX/m)),y1->(#1[[2]]/.t->time),vy1->((D[#1[[2]],t]/.t->time)+(kickY/m)),z1->(#1[[3]]/.t->time),vz1->((D[#1[[3]],t]/.t->time)+(kickZ/m))})}/.{t->(t-time)}/.#2&,{x0+v0x*t,y0+v0y*t,z0+v0z*t},{{kickX->kickListDown[[1,1]],kickY->kickListDown[[1,2]],kickZ->kickListDown[[1,3]],time->timeList[[1]]},{kickX->kickListDown[[2,1]],kickY->kickListDown[[2,2]],kickZ->kickListDown[[2,3]],time->timeList[[2]]},{kickX->kickListDown[[3,1]],kickY->kickListDown[[3,2]],kickZ->kickListDown[[3,3]],time->timeList[[3]]}}];

trajDownSeries=Series[trajDown/.{t->ϵ*t,T->ϵ*T},{ϵ,0,TOrder}]//Simplify

Prop Phase

evaluate Lagrangian and carry out action integrals.  Note that a time derivative reduces the number of powers of t, so each time derivative should should reduce the power of our small parameter ϵ by 1.

LUpList=Table[(L/.{x->trajUpSeries[[i+1,1]],y->trajUpSeries[[i+1,2]],z->trajUpSeries[[i+1,3]]}/.{Derivative[q_][f_][t]:>D[(1/(ϵ^q))f,{t,q}],f_[t]->f})//Normal//Expand,{i,1,Length[timeList]-1}];

integration over time adds a power of time, so in performing action integrals add a power of ϵ.

SUpList=(1/ℏ)*Table[Series[ϵ*Integrate[LUpList[[i]],{t,timeList[[i]],timeList[[i+1]]}],{ϵ,0,TOrder}],{i,1,Length[timeList]-1}];

SUp=SUpList//Total//Simplify;

LDownList=Table[(L/.{x->trajDownSeries[[i+1,1]],y->trajDownSeries[[i+1,2]],z->trajDownSeries[[i+1,3]]}/.{Derivative[q_][f_][t]:>D[(1/(ϵ^q))f,{t,q}],f_[t]->f})//Normal//Expand,{i,1,Length[timeList]-1}];

SDownList=(1/ℏ)*Table[Series[ϵ*Integrate[LDownList[[i]],{t,timeList[[i]],timeList[[i+1]]}],{ϵ,0,TOrder}],{i,1,Length[timeList]-1}];

SDown=SDownList//Total//Simplify;

propagation phase is difference between action integrals for upper and lower paths

ϕProp=SUp-SDown//Simplify

Laser Phase

for each momentum transferring interaction, evaluate imprinted laser phase at position and time of interaction

ϕLUp=Sum[(1/ℏ)*kickListUp[[i]].trajUpSeries[[i]]/.t->timeList[[i]],{i,1,Length[timeList]}];

ϕLDown=Sum[(1/ℏ)*kickListDown[[i]].trajDownSeries[[i]]/.t->timeList[[i]],{i,1,Length[timeList]}];

ϕLaser=ϕLUp-ϕLDown//Simplify

Separation Phase

calculate separation phase, making sure to use canonical momentum

pCanonical={D[L,x'[t]],D[L,y'[t]],D[L,z'[t]]}

pCanonicalFinalUp=Series[(pCanonical/.{x->trajUpSeries[[Length[timeList]+1,1]],y->trajUpSeries[[Length[timeList]+1,2]],z->trajUpSeries[[Length[timeList]+1,3]]}/.{Derivative[q_][f_][t]:>D[(1/(ϵ^q))f,{t,q}],f_[t]->t})//Normal//Expand,{ϵ,0,2}]/.t->Last[timeList];

pCanonicalFinalDown=Series[(pCanonical/.{x->trajDownSeries[[Length[timeList]+1,1]],y->trajDownSeries[[Length[timeList]+1,2]],z->trajDownSeries[[Length[timeList]+1,3]]}/.{Derivative[q_][f_][t]:>D[(1/(ϵ^q))f,{t,q}],f_[t]->f})//Normal//Expand,{ϵ,0,2}]/.t->Last[timeList];

pCanonicalFinalAvg=((pCanonicalFinalUp+pCanonicalFinalDown)/2)//Simplify

Δx=(trajUpSeries[[Length[timeList]]]-trajDownSeries[[Length[timeList]]])/.t->Last[timeList]//Simplify

ϕSep=-(pCanonicalFinalAvg/ℏ).Δx//Simplify

Total phase shift

sum of propagation and separation phase

ϕProp+ϕSep//Simplify

laser phase

ϕLaser//Simplify

total phase shift

ϕTot=(ϕLaser+ϕProp+ϕSep//Simplify//Normal)/.ϵ->1//Expand

ϕTotSeries=Series[ϕTot,{T,0,2}]//ExpandAll//Normal//Expand

numbers={k->((2*π)/(780.*(10^(-9)))),n->2,g->9.8,T->1}

ϕTotSeries/.numbers