Init

Some Functions

(* general pure state in polar coordinates *)
State[θ_,ϕ_]:={Cos[θ/2],Sin[θ/2]Exp[I ϕ]}
(*normalized state with real first argument *)
StateNorm[x_List]:=Module[{arg1=Arg[x[[1]]],arg2=Arg[x[[2]]],normx=x/Norm[x],θ,ϕ},
θ=2Chop@ArcCos[normx[[1]]Exp[-I arg1]];
ϕ=Mod[arg2-arg1,2π];
{Cos[θ/2],Sin[θ/2]Exp[I ϕ]}]
(* outputs the polar coordinate from a state vector *)
blochcoord[state_]:=Module[{arg1=Arg[state[[1]]],arg2=Arg[state[[2]]],normx=state/Norm[state],θ,ϕ},
θ=Abs[2Chop@ArcCos[normx[[1]]Exp[-I arg1]]];
ϕ=Abs@Mod[arg2-arg1,2π];
{θ,ϕ}]
blochcoordcart[state_]:={Sin[#[[1]]]Cos[#[[2]]],Sin[#[[1]]]Sin[#[[2]]],Cos[#[[1]]]}&[blochcoord[state]]

(* define a bloch sphere graph *)
With[{r=1},xy=ParametricPlot3D[{r Cos[θ],r Sin[θ],0},{θ,0,2 π},PlotStyle->{Dashed,Opacity[0.3]},Boxed->False,Axes->False];
yz=ParametricPlot3D[{0,r Cos[θ],r Sin[θ]},{θ,0,2 π},PlotStyle->{Dashed,Opacity[0.3]},Boxed->False,Axes->False];
xz=ParametricPlot3D[{r Cos[θ],0,r Sin[θ]},{θ,0,2 π},PlotStyle->{Dashed,Opacity[0.3]},Boxed->False,Axes->False];
bloch=Show[Graphics3D[{Blue,Opacity[.15],Sphere[{0,0,0},r],Black,Opacity[1],Thick,Line[{{0,r,0},{0,-r,0}}],Line[{{0,0,r},{0,0,-r}}],Line[{{r,0,0},{-r,0,0}}],Text[Style["z+",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,0,1.3 r}],Text[Style["z-",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,0,-1.3 r}],Text[Style["x+",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{1.3 r,0,0}],Text[Style["x-",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{-1.3 r,0,0}],Text[Style["y+",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,1.3 r,0}],Text[Style["y-",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,-1.3 r,0}]}],xy,xz,yz,Boxed->False]];
(* generate vector for bloch sphere representation *)
Options[blochvec]={color->Blue};
blochvec[{θ_,ϕ_},OptionsPattern[]]:=Graphics3D[{OptionValue[color],Thickness[0.02],Line[{{0,0,0},{Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]}}]}]
blocharrow[{θ_,ϕ_},opts:OptionsPattern[]]:=Arrow3D[{0,0,0},{Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]},opts]

σx=(({{0, 1}, {1, 0}}));σy=(({{0, -I}, {I, 0}}));σz=(({{1, 0}, {0, -1}}));
(* draw vector for bloch sphere taking an arbitrary pure input state *)
blochvecstate[state_]:=blochvec@blochcoord[state]
blocharrowstate[state_,opts:OptionsPattern[]]:=blocharrow[blochcoord[state],opts]
(*bloch vector for mixed states *)
blochvecmsphere[r_,θ_,ϕ_]:=Graphics3D[{Blue,Thickness[0.02],Line[{{0,0,0},r{Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]}}]}]
blochvecmsphere[{r_,θ_,ϕ_}]:=blochvec[r,θ,ϕ]
blochvecmarrowsphere[{r_,θ_,ϕ_},opts:OptionsPattern[]]:=Arrow3D[{0,0,0},r{Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]},opts]
(* cartesian coordinates *)
blochvecmcart[x_,y_,z_]:=Graphics3D[{Blue,Thickness[0.02],Line[{{0,0,0},{x,y,z}}]}]
blochvecmcart[{x_,y_,z_}]:=blochvecmcart[x,y,z]
blochvecmarrowcart[{x_,y_,z_},opts:OptionsPattern[]]:=Arrow3D[{0,0,0},{x,y,z},opts]
(* show vector with sphere *)
vectoronsphere[vec_]:=Show[bloch,vec,SphericalRegion->True]
(* Bloch-sphere coordinates from arbitrary state *)
MatrixToBloch[x_]:={Tr[σx.x],Tr[σy.x],Tr[σz.x]}

Urot[α_,t_,θ_,ϕ_]:=(({{Cos[α t/2] - I Cos[θ]Sin[α t /2], -I Exp[-I ϕ] Sin[θ] Sin[α t/2]}, {(-I Exp[I ϕ]) Sin[θ] Sin[α t/2], Cos[α t/2] + I Cos[θ]Sin[α t /2]}}))

RX[α_]:=Urot[α,1,π/2,0]
RY[α_]:=Urot[α,1,π/2,π/2]
RZ[α_]:=Urot[α,1,0,0]

SetOptions[ParametricPlot3D,PlotStyle->{Thick,Red},PlotRange->{{-1,1},{-1,1},{-1,1}},ImageSize->300,Axes->False,Boxed->False];
BlochSphere=Graphics3D[{Opacity[0.4],Sphere[{0,0,0},1]},Boxed->False];
BlochMesh=SphericalPlot3D[1,{θ,0,2 Pi},{ϕ,0,2Pi},MeshFunctions->{#4&},Mesh->{{Pi/2}},PlotStyle->Directive[Opacity[0]],Boxed->False,Axes->False];

Options[Arrow3D]={color->Red};
Arrow3D[a_,b_,OptionsPattern[]]:=Module[{l,n,rm,a3d,offset,col,arrowheadlength,thickness,arrowslope},
       col=OptionValue[color];
n=b-a;
l=Norm[n];
offset=l-1;
thickness=0.02;
arrowslope=0.3;
arrowheadlength=0.2;
rm=If[Abs[n]==l {0,0,1},Sign[n[[3]]]*IdentityMatrix[3],RotationMatrix[{{0,0,1},n}]];
If[l>0,
a3d={ParametricPlot3D[rm.{(arrowslope v) Cos[u],(arrowslope v) Sin[u],1-v+offset},{u,0,2 π},{v,0,arrowheadlength},Mesh->None,PlotStyle->Directive[col]],ParametricPlot3D[rm.{v Cos[u],v Sin[u],1-arrowheadlength+offset},{u,0,2 π},{v,0,arrowslope*arrowheadlength},PlotStyle->Directive[col],Mesh->None],Graphics3D[{col,Cylinder[{a,b-arrowheadlength/l (b-a)},thickness]}]},
a3d=Graphics3D[{col,Sphere[a,3thickness]}]];
a3d];

Question 1

Initial state

ψi={0,1}

vectoronsphere@blochvecstate[ψi]

1 a)

rotation by π about any axis on the equator :

MatrixExp[-I ω t σx/2].ψi

MatrixExp[-I ω t σx/2]//MatrixForm

Choose t = π/ω for a π - Pulse :

ψfa=MatrixExp[-I ω t σx/2].ψi/.t->π/ω

MatrixExp[-I ω t σx/2]/.t->π/ω//MatrixForm
MatrixExp[-I ω t σy/2]/.t->π/ω//MatrixForm

Show[vectoronsphere@blochvecstate[ψfa]]

1 b

rotation by π/2 about y-axis on the equator :

Choose t = π/(2ω) for a π/2 - Pulse :

MatrixExp[-I ω t (σx)/2]

ψfb=MatrixExp[-I ω t (σy)/2].ψi/.t->π/(2ω)

MatrixExp[I ω t σy/2]/.t->π/(2ω)//MatrixForm

Show[vectoronsphere@blochvecstate[ψfb]]

1 c

angle θ = 3 π/4 and the minus sign (phase factor of π) tells us that the state lies in the x-z plane

-> rotation by 3π/4 about y-axis on the equator :

Choose t = 3π/(4ω) for a 3π/4 - Pulse :

ψfc=MatrixExp[-I ω t σy/2].ψi/.t->3π/(4ω)

Cos[π/8]==Sin[3π/8]
Sin[π/8]==Cos[3π/8]

MatrixExp[-I ω t σy/2]/.t->3π/(4ω)//MatrixForm

Show[vectoronsphere@blochvecstate[ψfc]]

1 d

angle θ = 3 π/4 and a phase factor of π/4

Two options: 
either we rotate about the axis in the equatorial plane with an angle ϕ=π/4+π/2
or we rotate again about the y-axis and consequently apply a rotation about the z-axis

First option :

-> rotation by 3π/4 about ϕ-π/4-axis on the equator :

Choose t = 3π/(4ω) for a 3π/4 - Pulse :

σϕ[ϕ_]:=Cos[ϕ]σx+Sin[ϕ]σy//TrigToExp

σϕ[ϕ]
σϕ[0] (* rotation about x-axis *)
σϕ[π/2] (* rotation about y-axis *)

ψfd1=MatrixExp[-I ω t σϕ[ϕ]/2].ψi/.{t->3π/(4ω),ϕ->π/4}

(E^(I π/2))

MatrixExp[-I ω t σϕ[ϕ]/2].ψi/.{t->3π/(4ω),ϕ->π/4}//MatrixForm

Show[vectoronsphere@blochvecstate[ψfd1]]

Second option :

-> rotation by 3π/4 about y-axis followed by a z-Rotation for π/4:

U1=MatrixExp[-I ω t σy/2]/.t->3π/(4ω)
U2=MatrixExp[I ω t σz/2]/.t->π/(4ω)

Show[vectoronsphere@blochvecstate[U1.ψi]]

Show[vectoronsphere@blochvecstate[U2.U1.ψi]]

ψfd2= (E^((I π)/8)) U2.U1.ψi

U2.U1//MatrixForm

Question 2

Cos[θ]Cos[θ/2]+Sin[θ]Sin[θ/2]//FullSimplify

Sin[θ]Cos[θ/2]-Cos[θ]Sin[θ/2]//FullSimplify

Rotating Wave Approximation :

Hamiltonian (in lab frame)

H=((ℏ Ω)/2)σz+ℏ ϵ cos (ω t+ ϕ) σx

Larmor precession of superposition state |x> = (1/(Sqrt[2]))(({{1}, {1}})) for H=ℏΩ/2 σz

Manipulate[
vectoronsphere@blochvecstate[RZ[t].{(1/(Sqrt[2])),(1/(Sqrt[2]))}],{t,0,2π}]

Time - dependent field and its decomposition into rotating and counterrotating component:

Animate[Show[bloch,Arrow3D[{0,0,0},{Cos[t],0,0},color->Black],ImageSize->300],{t,0,2π},AnimationRunning->False]

Decompose state into field rotating clockwise and anti - clockwise:

Animate[Show[bloch,Arrow3D[{0,0,0},{Cos[t],0,0},color->Black],Arrow3D[{0,0,0},{1/2Cos[t],1/2Sin[t],0},color->Yellow],Arrow3D[{0,0,0},{1/2Cos[t],-1/2Sin[t],0},color->Blue],ImageSize->300],{t,0,2π-.1,0.1},AnimationRunning->False]

Transformation into rotating frame (rotating in the same direction as the Larmor precession (yellow arrow)):

Hamiltonian with  rotating and counterrotating field:

H=((ℏ Ω)/2)σz+((ℏ ϵ)/2) (cos (ω t+ ϕ) σx+sin(ω t + ϕ) σy)+((ℏ ϵ)/2) (cos (ω t+ ϕ) σx-sin(ω t + ϕ) σy)

Transformation matrix

U=(e^(i ωr t σz/2))

Hamiltonian in rotating frame

Hp = U H U† - i ℏ U (d U†/dt)

Urf=MatrixExp[I ((ωr t)/2)σz];

transformed σz

Urf.σz.Urf\[ConjugateTranspose]//ComplexExpand

transformed σx

Urf.σx.Urf\[ConjugateTranspose]//ComplexExpand//Simplify//MatrixForm
Cos[ωr t]σx-Sin[ωr t]σy//MatrixForm

transformed σy

Urf.σy.Urf\[ConjugateTranspose]//ComplexExpand//Simplify//MatrixForm
Cos[ωr t]σy+Sin[ωr t]σx//MatrixForm

Gauge Term

D[Urf\[ConjugateTranspose],t]

- I ℏ Urf D[Urf\[ConjugateTranspose],t]//ComplexExpand//MatrixForm
-(( ℏ ωr)/2)σz

first term in Hamiltonian (+) with ω = ωrf

termplus=Urf.(Cos[ωr t+ϕ]σx + Sin[ωr t+ϕ] σy).Urf\[ConjugateTranspose]//ComplexExpand//Simplify;
%//MatrixForm

Cos[ϕ]σx+Sin[ϕ]σy//MatrixForm

second term in Hamiltonian (-) with ω = ωrf

termminus=Urf.(Cos[ωr t+ϕ]σx - Sin[ωr t+ϕ] σy).Urf\[ConjugateTranspose]//ComplexExpand//Simplify;
%//MatrixForm

Discarding the terms at 2 ω amounts to making the Rotating Wave Approximation (RWA)

resulting Hamiltonian :

H=((ℏ (Ω-ω))/2)σz+(ϵ/2)Cos[ϕ]σx+(ϵ/2)Sin[ϕ]σy

Hrwa[ϵ_,Ω_,ω_,ϕ0_]:=(ϵ/2)(Cos[ϕ0])σx+(ϵ/2)( Sin[ϕ0])σy+(((Ω-ω))/2)σz

rotation  about effective magnetic field axis (ϵ/2 cos[ϕ],ϵ/2 sin[ϕ],(Ω-ω)/2):

Manipulate[Show[bloch,Arrow3D[{0,0,0},{ϵ Cos[ϕ],ϵ Sin[ϕ],(1-ω)},color->Darker@Orange]],
{ϕ,0,2π},
{ϵ,0,1},
{ω,0,2}]

Rotating about static axis in rotating frame

finalt=100;
Module[{ϵ=0.1,Ω=1,ω=1,ϕ0=2},
Animate[
vectoronsphere@blochvecstate@(MatrixExp[-I t Hrwa[ϵ,Ω,ω,ϕ0]].{0,1}),{t,0,finalt,finalt/100.},AnimationRunning->False]]

Numerical Solution without discarding counterrotating terms

Hcounter[ϵ_,Ω_,ω_,ϕ0_,t_]:=(ϵ/2)(Cos[ϕ0]+Cos[ϕ0+ 2 t ω])σx+(ϵ/2)(Sin[ϕ0]-Sin[ϕ0+ 2 t ω])σy+(((Ω-ω))/2)σz

finalt=100;
(* ϵ is typically small compared to Ω *)
Block[{ϵ=0.1,Ω=1,ω=1,ϕ0=0},
solncounter=u/.First@NDSolve[{Hcounter[ϵ,Ω,ω,ϕ0,t].u[t]==I u'[t],u[0]=={0,1}},u,{t,0,finalt},MaxSteps->(10^5)]]

(* complicated way in getting x,y,z coordinates and polar coordinates from the solution vector *)
Block[{interv},
interv=1000;
xaxiscounter=Table[t,{t,solncounter[[1,1,1]],solncounter[[1,1,2]],(solncounter[[1,1,2]]-solncounter[[1,1,1]])/interv}];
vecsolncounter=Table[{Cos[#2]Sin[#1],Sin[#2]Sin[#1],Cos[#1],#1,#2}&[Sequence@@blochcoord@solncounter[t]],{t,solncounter[[1,1,1]],solncounter[[1,1,2]],(solncounter[[1,1,2]]-solncounter[[1,1,1]])/interv}];
]
{xsolncounter,ysolncounter,zsolncounter,θsolncounter,ϕsolncounter}=Table[Interpolation@MapThread[{#1,#2}&,{xaxiscounter,vecsolncounter[[All,jj]]}],{jj,1,5}]

Animate[
Show[ParametricPlot3D[{xsolncounter[τ],ysolncounter[τ],zsolncounter[τ]},{τ,t-Min[10,t],t},PlotStyle->Directive[Blue,Thickness[0.01]]],
vectoronsphere@blochvec[{θsolncounter[t],ϕsolncounter[t]}],PlotRange->{{-1,1},{-1,1},{-1,1}}],
{t,100$MachineEpsilon,finalt},AnimationRate->5,AnimationRunning->False]

Without Rotating Wave Approximation :

No simple analytic solution, numerical simulation:

Hamiltonian :

Hnum[ϵ_,Ω_,ω_,ϕ0_,t_]=({ϵ Cos[ω t+ϕ0],0,Ω/2}.(({{σx}, {σy}, {σz}})))[[1]]

Solve the Schroedinger Equation numerically :

Block[{ϵ=0.1,Ω=1,ω=1,ϕ0=0},
solnrwa=u/.First@NDSolve[{Hnum[ϵ,Ω,ω,ϕ0,t].u[t]==I u'[t],u[0]=={0,1}},u,{t,0,finalt},MaxSteps->(10^5)]]

(* complicated way in getting x,y,z coordinates and polar coordinates from the solution vector *)
Block[{interv},
interv=1000;
xaxis=Table[t,{t,solnrwa[[1,1,1]],solnrwa[[1,1,2]],(solnrwa[[1,1,2]]-solnrwa[[1,1,1]])/interv}];
vecsolnrwa=Table[{Cos[#2]Sin[#1],Sin[#2]Sin[#1],Cos[#1],#1,#2}&[Sequence@@blochcoord@solnrwa[t]],{t,solnrwa[[1,1,1]],solnrwa[[1,1,2]],(solnrwa[[1,1,2]]-solnrwa[[1,1,1]])/interv}];
]
{xsolnrwa,ysolnrwa,zsolnrwa,θsolnrwa,ϕsolnrwa}=Table[Interpolation@MapThread[{#1,#2}&,{xaxis,vecsolnrwa[[All,jj]]}],{jj,1,5}]

Animate[
Show[ParametricPlot3D[{xsolnrwa[τ],ysolnrwa[τ],zsolnrwa[τ]},{τ,t-Min[10,t],t},PlotStyle->Directive[Blue,Thickness[0.01]]],
vectoronsphere@blochvec[{θsolnrwa[t],ϕsolnrwa[t]}],PlotRange->{{-1,1},{-1,1},{-1,1}}],
{t,100$MachineEpsilon,finalt},AnimationRate->3,AnimationRunning->False]