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

Schrödinger Cat : Unitary evolution

unit=(({{(1/(Sqrt[2])), 0, 0, -1/(Sqrt[2])}, {0, 1, 0, 0}, {0, 0, 1, 0}, {(1/(Sqrt[2])), 0, 0, 1/(Sqrt[2])}}))
unit.unit\[ConjugateTranspose]//MatrixForm

Simple Pauli Operators and their effect on states :

σz ... Phaseflip of basis state; more general: rotation by π about the z-axis

σz.{1,0}

{1,0} is an eigenvector to the eigenvalue +1

GraphicsRow[{vectoronsphere@blocharrowstate[{1,0}],vectoronsphere@blocharrowstate[σz.{1,0}]},ImageSize->700]

σz.{0,1}

{0, 1} is an eigenvector to the eigenvalue -1

GraphicsRow[{vectoronsphere@blocharrowstate[{0,1}],vectoronsphere@blocharrowstate[σz.{0,1}]},ImageSize->700]

factor -1 not visible on the Bloch sphere, since global phase is immaterial

superposition state on equatorial line is rotated about z-axis by an angle of π:

σz.{(1/(Sqrt[2])),(1/(Sqrt[2]))}

GraphicsRow[{vectoronsphere@blocharrowstate[{(1/(Sqrt[2])),(1/(Sqrt[2]))}],vectoronsphere@blocharrowstate[σz.{(1/(Sqrt[2])),(1/(Sqrt[2]))}]},ImageSize->700]

σx ... Spinflip of basis state; more general: rotation by π about the x-axis

σx.{0,1}

GraphicsRow[{vectoronsphere@blocharrowstate[{0,1}],vectoronsphere@blocharrowstate[σx.{0,1}]},ImageSize->700]

σx.{(1/(Sqrt[2])),(1/(Sqrt[2]))}

{(1/(Sqrt[2])),(1/(Sqrt[2]))} is eigenvector to eigenvalue +1

GraphicsRow[{vectoronsphere@blocharrowstate[{(1/(Sqrt[2])),(1/(Sqrt[2]))}],vectoronsphere@blocharrowstate[σx.{(1/(Sqrt[2])),(1/(Sqrt[2]))}]},ImageSize->700]

σx.{(1/(Sqrt[2])),I(1/(Sqrt[2]))}

GraphicsRow[{vectoronsphere@blocharrowstate[{(1/(Sqrt[2])),I(1/(Sqrt[2]))}],vectoronsphere@blocharrowstate[σx.{(1/(Sqrt[2])),I(1/(Sqrt[2]))}]},ImageSize->700]

σy ... Spin+Phaseflip

σy.{0,1}//TrigToExp

again, phaseflip is not visible on Blochsphere for ground/excited state, since global phase is immaterial

GraphicsRow[{vectoronsphere@blocharrowstate[{0,1}],vectoronsphere@blocharrowstate[σy.{0,1}]},ImageSize->700]

GraphicsRow[{vectoronsphere@blocharrowstate[{(1/(Sqrt[2])),(1/(Sqrt[2]))}],vectoronsphere@blocharrowstate[σy.{(1/(Sqrt[2])),(1/(Sqrt[2]))}]},ImageSize->700]

σy.{(1/(Sqrt[2])),(I/(Sqrt[2]))}

{(1/(Sqrt[2])),(i/(Sqrt[2]))} is eigenvector to eigenvalue +1

GraphicsRow[{vectoronsphere@blocharrowstate[{(1/(Sqrt[2])),(I/(Sqrt[2]))}],vectoronsphere@blocharrowstate[σy.{(1/(Sqrt[2])),(I/(Sqrt[2]))}]},ImageSize->700]

connection to matrix exponentiation :

MatrixExp[-I ((ω t)/2)σx]/.{t->π/ω}

is equivalent to the operation (e^(-i π/2))σx.

Projectors onto eigenstates:

eg. 
|±><±|=(Subscript[P, z±])=1/2(1±σz):

In matrix representation :

(Subscript[P, z+])=(({{1, 0}, {0, 0}}))
(Subscript[P, z-])=(({{0, 0}, {0, 1}}))

Describe the simplest possible observables corresponding to yes-no questions.
Example: The question "Is the system in the groundstate?" corresponds to applying the projector  | 0⟩⟨0|. After this application the system is in the state   | 0⟩.  The expectation value ⟨ψ | 0⟩⟨0|ψ⟩ tells us what the probability is to find the system in the ground state.

Properties :

*) (P^2)=P 
*) Eigenvalues are either 0 or 1

Pzp=(({{1, 0}, {0, 0}}));
Pzm=(({{0, 0}, {0, 1}}));

Pzp.(({{1}, {0}}))

has eigenvalue 1, i.e. the system is found with certainty in state (({{1}, {0}})).

GraphicsRow[{Show[bloch,blocharrowstate[{1,0}]],Show[bloch,blocharrowstate[{1,0}],blocharrowstate[{1,0},color->Green]]},ImageSize->700]

Pzp.(({{0}, {1}}))

has eigenvalue 0, i.e. the system is never found in state (({{0}, {1}})).

Show[bloch,blocharrowstate[{0,1}],blocharrowstate[{1,0},color->Green]]

xpstate=1/(Sqrt[2]){1,0}+1/(Sqrt[2]){0,1}

has both components in z + and z - direction

vectoronsphere@blocharrowstate[xpstate]

project onto eigenstates :

Pzp.xpstate
Pzm.xpstate

probability to find system in either of the basis states :

Conjugate[xpstate].Pzp.xpstate
Conjugate[xpstate].Pzm.xpstate

Problem Sheet 3 :

Problem Sheet 3 - Question 1

σv[{vx_,vy_,vz_}]:=(({{vz, vx-I vy}, {vx+I vy, -vz}}))

π-flip about axis defined by (vx, vy, vz) :

Module[{vector},
vector={0,0,1};
Show[bloch,blocharrowstate[{(1/(Sqrt[2])),(I/(Sqrt[2]))}],blocharrowstate[σv[vector].{(1/(Sqrt[2])),(I/(Sqrt[2]))},color->Green],blochvecmarrowcart[vector,color->Darker@Blue]]]

vector in polar coordinates :

Module[{vector},
θ=π/3.;
ϕ=0.3;
vector={Cos[ϕ]Sin[θ],Sin[ϕ]Sin[θ],Cos[θ]};
Show[bloch,blocharrowstate[{(1/(Sqrt[2])),(I/(Sqrt[2]))}],blocharrowstate[σv[vector].{(1/(Sqrt[2])),(I/(Sqrt[2]))},color->Green],blochvecmarrowcart[vector,color->Darker@Blue]]]

Eigenvalues :

Simplify[
Eigenvalues[(({{vz, vx+I vy}, {vx-I vy, -vz}}))],Assumptions->((vx^2)+(vy^2)+(vz^2))==1]

Eigenvectors

evs=FullSimplify[
Eigenvectors[(({{vz, vx+I vy}, {vx-I vy, -vz}}))/.{vx->Sin[θ]Cos[ϕ],vy->Sin[θ]Sin[ϕ],vz->Cos[θ]}]
,Assumptions->{Im[θ]==0,Im[ϕ]==0}]

Normalize the eigenvectors :

evs1=evs[[1]]
evs2=evs[[2]]

norm1=FullSimplify[(Sqrt[Conjugate[evs1].evs1]),Assumptions->{Im[θ]==0,Im[ϕ]==0,0<θ<π}]
norm2=FullSimplify[(Sqrt[Conjugate[evs2].evs2]),Assumptions->{Im[θ]==0,Im[ϕ]==0,0<θ<π}]

nevs1=evs1/norm1
nevs2=evs2/norm2

Visualize them on the Blochsphere

GraphicsRow[{vectoronsphere@blocharrowstate[nevs1/.{θ->π/4,ϕ->0}],vectoronsphere@blocharrowstate[nevs2/.{θ->π/4,ϕ->0}]},ImageSize->700]

Projector onto eigenstate in arbitrary direction :

(Subscript[P, n])=1/2(1+(OverVector[n]).(OverVector[σ]))

Pp[θ_,ϕ_]:=(1/2)((({{1, 0}, {0, 1}}))+{Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]}.{σx,σy,σz})
Pm[θ_,ϕ_]:=(1/2)((({{1, 0}, {0, 1}}))-{Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]}.{σx,σy,σz})

State[π/2,0]
Pp[π/2,0].State[π/2,0]
Pm[π/2,0].State[π/2,0]

Problem Sheet 3 - Question 2

B-field in z-direction

Hamilton Operator :

σz//MatrixForm

H[A_]:= A σz

The unitary matrix can be calculated  by using MatrixExp:

Rz[A_,t_]:=MatrixExp[-I A σz t]

Rz[(Subscript[μ, n])B0,t]//MatrixForm

groundstate={0,1};
excitedstate={1,0};
xpstate={(1/(Sqrt[2])),(1/(Sqrt[2]))};

Manipulate[
vectoronsphere@blochvecstate[Rz[1,t].{0,1.}],{t,0,2π}]

B-field in x-direction

Hamilton Operator :

Hx[B_]:= B σx

exponentiation:

Rx[B_,t_]:=MatrixExp[-I Hx[B]t]

Rx[(Subscript[μ, n])B0,t]//MatrixForm

rotstate=Rx[1,t].groundstate

Projection to excited state :

Pzp=1/2((({{1, 0}, {0, 1}}))+σz)

Pzp.rotstate

Expectation value :

rotstate\[Conjugate].Pzp.rotstate//ComplexExpand

Rabi Oscillations :

Manipulate[
GraphicsRow[{Show[bloch,
blocharrowstate[Rx[1,t].groundstate],blochvecmarrowcart[{1,0,0},color->Darker@Blue]],
Plot[((Abs[(Rx[1,t0].groundstate)[[1]]])^2),{t0,0,t},PlotRange->{{0,2π},{0,1}}]},ImageSize->700],{t,0.0000001,2π}]

B-field in y - direction :

-> Rotation about y axis instead of rotation about x - axis

Ry[A_,t_]:=MatrixExp[-I A σy t];

Manipulate[
GraphicsRow[{Show[bloch,
blocharrowstate[Ry[1,t].groundstate],blochvecmarrowcart[{0,1,0},color->Darker@Blue]],
Plot[((Abs[(Ry[1,t0].groundstate)[[1]]])^2),{t0,0,t},PlotRange->{{0,2π},{0,1}}]},ImageSize->700],{t,0.0000001,2π}]

Rotating Wave Approximation :

Larmor Precession of superposition state :

Rz[1,t]

Animate[Show[bloch,blocharrowstate[Rz[1,t].xpstate],SphericalRegion->True,ImageSize->300],{t,0,2π},AnimationRunning->False]

Time - dependent field :

Animate[Show[bloch,Arrow3D[{0,0,0},{Cos[t],0,0},color->Green],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->Green],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]

make a transformation into a frame rotating in the same direction as the Larmor precession (yellow arrow):

Trans=MatrixExp[I ((ωr t)/2)σz];
termplus=Trans.(Cos[ω t]σx + Sin[ω t] σy).Trans\[ConjugateTranspose]//Simplify[#,Assumptions->{t∈Reals,ωr ∈Reals,ωr==ω}]&
termminus=Trans.(Cos[ω t]σx - Sin[ω t] σy).Trans\[ConjugateTranspose]//Simplify[#,Assumptions->{t∈Reals,ωr ∈Reals,ωr==ω}]&;
1/2(Tr[termplus.σx]"σx"+Tr[termplus.σy]"σy")
1/2(Tr[termminus.σx]"σx"+Tr[termminus.σy]"σy")//ExpToTrig

Discarding the terms at 2ω amounts to making the Rotating Wave Approximation :

Animate[vectoronsphere@blocharrowstate@(RX[3t].{1,0}),{t,0,2π/3-0.03,0.03},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[{ϵ=3,Ω=70,ω=70,ϕ0=0},
solnrwa=u/.First@NDSolve[{Hnum[ϵ,Ω,ω,ϕ0,t].u[t]==I u'[t],u[0]=={1,0}},u,{t,0,2π/ϵ},StartingStepSize->0.01,MaxStepSize->(10^(-5)),MaxSteps->(10^6)]]

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[.3,t],t}],
vectoronsphere@blocharrow[{θsolnrwa[t],ϕsolnrwa[t]}]],
{t,100$MachineEpsilon,2π/3-0.005},AnimationRate->0.1,AnimationRunning->False]