Two Level System and Rabi Flopping

sol=DSolve[{c1'[t]==-I*(Ω/2)*Exp[-I*ϕ]*Exp[I*δ*t]*c2[t],c2'[t]==-I*(Ω/2)*Exp[I*ϕ]*Exp[-I*δ*t]*c1[t],c1[0]==1,c2[0]==0},{c1[t],c2[t]},t][[1]];

δ=0

Plot[{((Abs[c1[t]])^2),((Abs[c2[t]])^2)}/.sol/.{Ω->1,δ->0,ϕ->0}//Evaluate,{t,0,4*π},PlotStyle->{Red,Blue}]

δ=Ω

Plot[{((Abs[c1[t]])^2),((Abs[c2[t]])^2)}/.sol/.{Ω->1,δ->1,ϕ->0}//Evaluate,{t,0,4*π},PlotStyle->{Red,Blue}]

δ=2*Ω

Plot[{((Abs[c1[t]])^2),((Abs[c2[t]])^2)}/.sol/.{Ω->1,δ->2,ϕ->0}//Evaluate,{t,0,4*π},PlotStyle->{Red,Blue}]

AC Stark Shifts

H=ℏ*{{0,(Ω/2)*Exp[-I*ϕ]},{(Ω/2)*Exp[I*ϕ],-δ}}

H//MatrixForm

FullSimplify[Eigensystem[H]/.{(Sqrt[(E^(2 I ϕ)) ((δ^2)+(Ω^2))])->(E^(I*ϕ))*(Sqrt[(δ^2)+(Ω^2)])}]

Bragg Diffraction Simulations

time in units of (1/(4 (Subscript[ω, r]))), frequencies in units of 4 (Subscript[ω, r])

Two-Photon Bragg

very efficient π pulse for low Rabi frequency

Ω=(0.3/4);
PosMomentumStates=15;
NegMomentumStates=5;
T=100;
coefs=Table[(Subscript[g, n])[τ],{n,-(NegMomentumStates+1),PosMomentumStates+1}];
dq=Table[(Subscript[g, n])'[τ]==-I*(Ω/2)*(Exp[I*1*τ])*(Subscript[g, n+1])[τ]-I*(Ω/2)*(Exp[-I*1*τ])*(Subscript[g, n-1])[τ]-I*((n^2))*(Subscript[g, n])[τ],{n,-NegMomentumStates,PosMomentumStates}];
initial=Join[{(Subscript[g, -(NegMomentumStates+1)])'[τ]==0,(Subscript[g, -(NegMomentumStates+1)])[0]==0, (Subscript[g, (PosMomentumStates+1)])'[τ]==0, (Subscript[g, (PosMomentumStates+1)])[0]==0}, Table[(Subscript[g, n])[0]==If[n==0,1,0],{n,-NegMomentumStates,PosMomentumStates}]];
eqs=Join[dq,initial];
dqsol=NDSolve[eqs,coefs,{τ,0,100},MaxSteps->22000000];

plot1=Plot[((((coefs/.dqsol[[1]])//Abs))^2)//Evaluate,{τ,0,(π/Ω)},PlotRange->{0,1.025},Frame->True,PlotStyle->Table[{Hue[((17*n)/(Length[coefs]))],Thickness[0.004]},{n,1,Length[coefs]}]]

((Abs[(coefs/.dqsol[[1]])[[NegMomentumStates+1+2]]/.τ->(π/Ω)])^2)

other levels start to get populated as Rabi frequency increases

Ω=(1/4);
PosMomentumStates=15;
NegMomentumStates=5;
T=100;
coefs=Table[(Subscript[g, n])[τ],{n,-(NegMomentumStates+1),PosMomentumStates+1}];
dq=Table[(Subscript[g, n])'[τ]==-I*(Ω/2)*(Exp[I*1*τ])*(Subscript[g, n+1])[τ]-I*(Ω/2)*(Exp[-I*1*τ])*(Subscript[g, n-1])[τ]-I*((n^2))*(Subscript[g, n])[τ],{n,-NegMomentumStates,PosMomentumStates}];
initial=Join[{(Subscript[g, -(NegMomentumStates+1)])'[τ]==0,(Subscript[g, -(NegMomentumStates+1)])[0]==0, (Subscript[g, (PosMomentumStates+1)])'[τ]==0, (Subscript[g, (PosMomentumStates+1)])[0]==0}, Table[(Subscript[g, n])[0]==If[n==0,1,0],{n,-NegMomentumStates,PosMomentumStates}]];
eqs=Join[dq,initial];
dqsol=NDSolve[eqs,coefs,{τ,0,100},MaxSteps->22000000];

plot1=Plot[((((coefs/.dqsol[[1]])//Abs))^2)//Evaluate,{τ,0,(π/Ω)},PlotRange->{0,1.025},Frame->True,PlotStyle->Table[{Hue[((17*n)/(Length[coefs]))],Thickness[0.004]},{n,1,Length[coefs]}]]

((Abs[(coefs/.dqsol[[1]])[[NegMomentumStates+1+2]]/.τ->(π/Ω)])^2)

Ω=(4/4);
PosMomentumStates=15;
NegMomentumStates=5;
T=100;
coefs=Table[(Subscript[g, n])[τ],{n,-(NegMomentumStates+1),PosMomentumStates+1}];
dq=Table[(Subscript[g, n])'[τ]==-I*(Ω/2)*(Exp[I*1*τ])*(Subscript[g, n+1])[τ]-I*(Ω/2)*(Exp[-I*1*τ])*(Subscript[g, n-1])[τ]-I*((n^2))*(Subscript[g, n])[τ],{n,-NegMomentumStates,PosMomentumStates}];
initial=Join[{(Subscript[g, -(NegMomentumStates+1)])'[τ]==0,(Subscript[g, -(NegMomentumStates+1)])[0]==0, (Subscript[g, (PosMomentumStates+1)])'[τ]==0, (Subscript[g, (PosMomentumStates+1)])[0]==0}, Table[(Subscript[g, n])[0]==If[n==0,1,0],{n,-NegMomentumStates,PosMomentumStates}]];
eqs=Join[dq,initial];
dqsol=NDSolve[eqs,coefs,{τ,0,100},MaxSteps->22000000];

plot1=Plot[((((coefs/.dqsol[[1]])//Abs))^2)//Evaluate,{τ,0,(π/Ω)},PlotRange->{0,1.025},Frame->True,PlotStyle->Table[{Hue[((17*n)/(Length[coefs]))],Thickness[0.004]},{n,1,Length[coefs]}]]

((Abs[(coefs/.dqsol[[1]])[[NegMomentumStates+1+2]]/.τ->(π/Ω)])^2)

Ω=(8/4);
PosMomentumStates=15;
NegMomentumStates=5;
T=100;
coefs=Table[(Subscript[g, n])[τ],{n,-(NegMomentumStates+1),PosMomentumStates+1}];
dq=Table[(Subscript[g, n])'[τ]==-I*(Ω/2)*(Exp[I*1*τ])*(Subscript[g, n+1])[τ]-I*(Ω/2)*(Exp[-I*1*τ])*(Subscript[g, n-1])[τ]-I*((n^2))*(Subscript[g, n])[τ],{n,-NegMomentumStates,PosMomentumStates}];
initial=Join[{(Subscript[g, -(NegMomentumStates+1)])'[τ]==0,(Subscript[g, -(NegMomentumStates+1)])[0]==0, (Subscript[g, (PosMomentumStates+1)])'[τ]==0, (Subscript[g, (PosMomentumStates+1)])[0]==0}, Table[(Subscript[g, n])[0]==If[n==0,1,0],{n,-NegMomentumStates,PosMomentumStates}]];
eqs=Join[dq,initial];
dqsol=NDSolve[eqs,coefs,{τ,0,100},MaxSteps->22000000];

plot1=Plot[((((coefs/.dqsol[[1]])//Abs))^2)//Evaluate,{τ,0,(π/Ω)},PlotRange->{0,1.025},Frame->True,PlotStyle->Table[{Hue[((17*n)/(Length[coefs]))],Thickness[0.004]},{n,1,Length[coefs]}]]

((Abs[(coefs/.dqsol[[1]])[[NegMomentumStates+1+2]]/.τ->(π/Ω)])^2)

It can help to use Gaussian shaped pulses instead of square pulses, though still do not want to make Rabi frequency too high

Ω=(8/4)*(E^(-((((τ-50))^2)/(1.2^2))));
PosMomentumStates=15;
NegMomentumStates=5;
T=100;
coefs=Table[(Subscript[g, n])[τ],{n,-(NegMomentumStates+1),PosMomentumStates+1}];
dq=Table[(Subscript[g, n])'[τ]==-I*(Ω/2)*(Exp[I*1*τ])*(Subscript[g, n+1])[τ]-I*(Ω/2)*(Exp[-I*1*τ])*(Subscript[g, n-1])[τ]-I*((n^2))*(Subscript[g, n])[τ],{n,-NegMomentumStates,PosMomentumStates}];
initial=Join[{(Subscript[g, -(NegMomentumStates+1)])'[τ]==0,(Subscript[g, -(NegMomentumStates+1)])[0]==0, (Subscript[g, (PosMomentumStates+1)])'[τ]==0, (Subscript[g, (PosMomentumStates+1)])[0]==0}, Table[(Subscript[g, n])[0]==If[n==0,1,0],{n,-NegMomentumStates,PosMomentumStates}]];
eqs=Join[dq,initial];
dqsol=NDSolve[eqs,coefs,{τ,0,100},MaxSteps->22000000];

plot1=Plot[((((coefs/.dqsol[[1]])//Abs))^2)//Evaluate,{τ,0,100},PlotRange->{0,1.025},Frame->True,PlotStyle->Table[{Hue[((17*n)/(Length[coefs]))],Thickness[0.004]},{n,1,Length[coefs]}]]

((Abs[(coefs/.dqsol[[1]])[[NegMomentumStates+1+2]]/.τ->100])^2)

Ω=(4/4)*(E^(-((((τ-50))^2)/(1.8^2))));
PosMomentumStates=15;
NegMomentumStates=5;
T=100;
coefs=Table[(Subscript[g, n])[τ],{n,-(NegMomentumStates+1),PosMomentumStates+1}];
dq=Table[(Subscript[g, n])'[τ]==-I*(Ω/2)*(Exp[I*1*τ])*(Subscript[g, n+1])[τ]-I*(Ω/2)*(Exp[-I*1*τ])*(Subscript[g, n-1])[τ]-I*((n^2))*(Subscript[g, n])[τ],{n,-NegMomentumStates,PosMomentumStates}];
initial=Join[{(Subscript[g, -(NegMomentumStates+1)])'[τ]==0,(Subscript[g, -(NegMomentumStates+1)])[0]==0, (Subscript[g, (PosMomentumStates+1)])'[τ]==0, (Subscript[g, (PosMomentumStates+1)])[0]==0}, Table[(Subscript[g, n])[0]==If[n==0,1,0],{n,-NegMomentumStates,PosMomentumStates}]];
eqs=Join[dq,initial];
dqsol=NDSolve[eqs,coefs,{τ,0,100},MaxSteps->22000000];

plot1=Plot[((((coefs/.dqsol[[1]])//Abs))^2)//Evaluate,{τ,0,100},PlotRange->{0,1.025},Frame->True,PlotStyle->Table[{Hue[((17*n)/(Length[coefs]))],Thickness[0.004]},{n,1,Length[coefs]}]]

((Abs[(coefs/.dqsol[[1]])[[NegMomentumStates+1+2]]/.τ->100])^2)

Six-Photon Bragg

Ω=(7.3/4)*(E^(-((((τ-50))^2)/(10^2))));
PosMomentumStates=15;
NegMomentumStates=5;
T=100;
coefs=Table[(Subscript[g, n])[τ],{n,-(NegMomentumStates+1),PosMomentumStates+1}];
dq=Table[(Subscript[g, n])'[τ]==-I*(Ω/2)*(Exp[I*3*τ])*(Subscript[g, n+1])[τ]-I*(Ω/2)*(Exp[-I*3*τ])*(Subscript[g, n-1])[τ]-I*((n^2))*(Subscript[g, n])[τ],{n,-NegMomentumStates,PosMomentumStates}];
initial=Join[{(Subscript[g, -(NegMomentumStates+1)])'[τ]==0,(Subscript[g, -(NegMomentumStates+1)])[0]==0, (Subscript[g, (PosMomentumStates+1)])'[τ]==0, (Subscript[g, (PosMomentumStates+1)])[0]==0}, Table[(Subscript[g, n])[0]==If[n==0,1,0],{n,-NegMomentumStates,PosMomentumStates}]];
eqs=Join[dq,initial];
dqsol=NDSolve[eqs,coefs,{τ,0,100},MaxSteps->22000000];

plot1=Plot[((((coefs/.dqsol[[1]])//Abs))^2)//Evaluate,{τ,0,100},PlotRange->{0,1.025},Frame->True,PlotStyle->Table[{Hue[((17*n)/(Length[coefs]))],Thickness[0.004]},{n,1,Length[coefs]}]]

((Abs[(coefs/.dqsol[[1]])[[NegMomentumStates+3+2]]/.τ->100])^2)

Twelve-Photon Bragg

Ω=(33/4)*(E^(-((((τ-50))^2)/(10^2))));
PosMomentumStates=15;
NegMomentumStates=5;
T=100;
coefs=Table[(Subscript[g, n])[τ],{n,-(NegMomentumStates+1),PosMomentumStates+1}];
dq=Table[(Subscript[g, n])'[τ]==-I*(Ω/2)*(Exp[I*6*τ])*(Subscript[g, n+1])[τ]-I*(Ω/2)*(Exp[-I*6*τ])*(Subscript[g, n-1])[τ]-I*((n^2))*(Subscript[g, n])[τ],{n,-NegMomentumStates,PosMomentumStates}];
initial=Join[{(Subscript[g, -(NegMomentumStates+1)])'[τ]==0,(Subscript[g, -(NegMomentumStates+1)])[0]==0, (Subscript[g, (PosMomentumStates+1)])'[τ]==0, (Subscript[g, (PosMomentumStates+1)])[0]==0}, Table[(Subscript[g, n])[0]==If[n==0,1,0],{n,-NegMomentumStates,PosMomentumStates}]];
eqs=Join[dq,initial];
dqsol=NDSolve[eqs,coefs,{τ,0,100},MaxSteps->22000000];

plot1=Plot[((((coefs/.dqsol[[1]])//Abs))^2)//Evaluate,{τ,0,100},PlotRange->{0,1.025},Frame->True,PlotStyle->Table[{Hue[((17*n)/(Length[coefs]))],Thickness[0.004]},{n,1,Length[coefs]}]]

((Abs[(coefs/.dqsol[[1]])[[NegMomentumStates+6+2]]/.τ->100])^2)