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

Mixed States

(* groundstate *) 
gs={0,1};
(* excited state *)
es={1,0};
(* superpositions state *)
xps={1,1}/(Sqrt[2]);
xms={1,-1}/(Sqrt[2]);
(* generate a projector matrix from a pure state *)
tps[x_]:=Outer[Times,Conjugate[x],x]
tps[x_,y_]:=Outer[Times,Conjugate[x],y]

State[π/10,0]//MatrixForm

tps[State[π/10,0]]//N

ρarb=0.5tps[es]+0.5 tps[State[π/10,0]]
ρarb//MatrixForm

Properties :

*) Hermitian (all eigenvalues are real)
*) Tr[ρ]=1 (normalization)
*) positive (all eigenvalues are positive)

*) Tr[(ρ^2)]≤Tr[ρ], where equality holds for ρ representing a pure state.
*) State evolution by U ρ (SuperDagger[U]).

e.g.ρ=(({{0.5, 0}, {0, 0.5}}))

ρ1=0.5 tps[gs,gs]+0.5 tps[es,es]
ρ1//MatrixForm

Test properties

{ρ1\[ConjugateTranspose]==ρ1,Tr[ρ1],Eigenvalues[ρ1]}

{ρarb\[ConjugateTranspose]==ρarb,Tr[ρarb],Eigenvalues[ρarb]}

Measurement :

Tr[ρ.M] with ρ the density matrix and M the measurement operator (e.g. (Subscript[σ, x]),(Subscript[σ, y]),(Subscript[σ, z]),(Subscript[P, x+]),...)

e.g. probability to find the system in the mixed state ρarb in the ground state

Tr[ρarb.(({{0, 0}, {0, 1}}))]

Expectation value of σz :

Tr[ρarb.σz]

Density Matrix representation on the Bloch ball:

Mixed states can also be depicted on the bloch ball using the parameterization
ρ=((I+(OverVector[r]).(OverVector[σ]))/2)with ||(OverVector[r])||≤1. (OverVector[r]) gives the Bloch vector.

Determine (OverVector[r]) by measuring the expecation value of the Pauli matrices:

Since any Hermitian matrix can unambiguously be written as a linear combination of Pauli matrices we can calculate the coefficients by taking the trace:

(Subscript[r, x])=Tr[(Subscript[ρσ, x])]=Tr[((id (Subscript[σ, x]))/2)]+Tr[(1/2) (Subscript[r, 1]) (Subscript[σ, x]) (Subscript[σ, x])]+Tr[(1/2) (Subscript[r, 2]) (Subscript[σ, y]) (Subscript[σ, x])]+Tr[(1/2) (Subscript[r, 3]) (Subscript[σ, z]) (Subscript[σ, x])]=0+(Subscript[r, 1])+0+0=(Subscript[r, 1])

(* Bloch-sphere coordinates from arbitrary state *)
MatrixToBloch[x_]:={Tr[σx.x],Tr[σy.x],Tr[σz.x]}

pure state | 1⟩ ⟨1 |:

vectoronsphere@blochvecmarrowcart@MatrixToBloch[(({{1, 0}, {0, 0}}))]

incoherent superposition of 3/4 | 0⟩ ⟨0 | and 1/4 | 1⟩ ⟨1 |:

vectoronsphere@blochvecmarrowcart@MatrixToBloch[1/2(({{3/4, 0}, {0, 1/4}}))]

Totally mixed state - incoherent superposition of |0⟩⟨0| and | 1⟩ ⟨1 |:

vectoronsphere@blochvecmarrowcart@MatrixToBloch[1/2(({{1, 0}, {0, 1}}))]

incoherent superposition from above:

vectoronsphere@blochvecmarrowcart@MatrixToBloch[ρarb]

arbitrary mixed state :

Manipulate[vectoronsphere@blochvecmarrowcart@MatrixToBloch[1/2((({{1, 0}, {0, 1}}))+r {Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]}.{σx,σy,σz})],{r,0,1},{θ,0,π},{ϕ,0,2π}]

State Evolution :

U.ρ.U† with U the unitary evolution operator

unitary[θ_]:=MatrixExp[-I (θ/2)σx]
Manipulate[vectoronsphere@blochvecmarrowcart@MatrixToBloch[1/2unitary[θ].(({{7/8, 0}, {0, 1/8}})).unitary[θ]\[ConjugateTranspose]],{θ,0,2π}]

Bloch equations :

BE=MapThread[Equal,{{x'[τ],y'[τ],z'[τ]},Cross[{x[τ],y[τ],z[τ]},{ϵx[τ],ϵy[τ],ϵz[τ]}]}];
BE//MatrixForm

DSolve[{r'[t]==(r[t]+1)/T1,r[0]==1},r[t],t]

inital state = excited state; B - field in z - direction

Block[{ϵx,ϵy,ϵz,x,y,z,ini,tini,tfinal,rx,ry,rz},

ϵx[τ_]=0;ϵy[τ_]=0;ϵz[τ_]=1;

ini={x[0]==0,y[0]==0,z[0]==1};

tini=0; tfinal=2π;
v1[τ_]={x[τ],y[τ],z[τ]}/.First@NDSolve[{Sequence@@BE,Sequence@@ini},{x[τ],y[τ],z[τ]},{τ,tini,tfinal}];
Animate[vectoronsphere@blochvecmarrowcart[v1[t]],{t,tini,tfinal},AnimationRunning->False]
]

inital state = excited state; B - field in x - direction

Block[{ϵx,ϵy,ϵz,x,y,z,ini,tini,tfinal,rx,ry,rz},

ϵx[τ_]=1;ϵy[τ_]=0;ϵz[τ_]=0;

ini={x[0]==0,y[0]==0,z[0]==1};

tini=0; tfinal=2π;
v2[τ_]={x[τ],y[τ],z[τ]}/.First@NDSolve[{Sequence@@BE,Sequence@@ini},{x[τ],y[τ],z[τ]},{τ,tini,tfinal}];
Animate[vectoronsphere@blochvecmarrowcart[v2[t]],{t,tini,tfinal},AnimationRunning->False]
]

with decay :

equs=MapThread[Equal,{{x'[τ],y'[τ],z'[τ]},Cross[{x[τ],y[τ],z[τ]},{ϵx[τ],ϵy[τ],ϵz[τ]}]-{((x[τ])/T2),((y[τ])/T2),(((z[τ]-zs))/T1)}}];
equs//MatrixForm

For (OverVector[ϵ])=(0,0,1), decay of excited state

Block[{ϵx,ϵy,ϵz,x,y,z,T1,T2,zs,ini,tini,tfinal,rx,ry,rz},

ϵx[τ_]=0;ϵy[τ_]=0;ϵz[τ_]=1;

zs=-1; (* steady state *)
ini={x[0]==0,y[0]==0,z[0]==1};

T1=1;
T2=.0001;

tini=0; tfinal=2π;
vec1[τ_]={x[τ],y[τ],z[τ]}/.First@NDSolve[{Sequence@@equs,Sequence@@ini},{x[τ],y[τ],z[τ]},{τ,tini,tfinal}];
Animate[vectoronsphere@blochvecmarrowcart[vec1[t]],{t,tini,tfinal},AnimationRunning->False]
]

For (OverVector[ϵ])=(0,0,1), decay of ground state

Block[{ϵx,ϵy,ϵz,x,y,z,T1,T2,zs,ini,tini,tfinal},

ϵx[τ_]=0;ϵy[τ_]=0;ϵz[τ_]=1;

zs=-1; (* steady state *)
ini={x[0]==0,y[0]==0,z[0]==-1};

T1=1;
T2=0.001;

tini=0; tfinal=2π;
vec2[τ_]={x[τ],y[τ],z[τ]}/.First@NDSolve[{Sequence@@equs,Sequence@@ini},{x[τ],y[τ],z[τ]},{τ,tini,tfinal}];
Animate[vectoronsphere@blochvecmarrowcart[vec2[t]],{t,tini,tfinal},AnimationRunning->False]
]

For (OverVector[ϵ])=(0,0,1), dephasing of superposition state

Block[{ϵx,ϵy,ϵz,x,y,z,T1,T2,zs,ini,tini,tfinal},

ϵx[τ_]=0;ϵy[τ_]=0;ϵz[τ_]=3;

zs=-1; (* steady state *)
ini={x[0]==0,y[0]==1,z[0]==0};

T1=∞;
T2=2;

tini=0.; tfinal=2π;
vec3[τ_]={x[τ],y[τ],z[τ]}/.First@NDSolve[{Sequence@@equs,Sequence@@ini},{x[τ],y[τ],z[τ]},{τ,tini,tfinal}];
Animate[Show[ParametricPlot3D[vec3[θ],{θ,0,a}],vectoronsphere@blochvecmarrowcart[vec3[a]]],{{a,$MachineEpsilon},tini,tfinal},AnimationRunning->False]
]

For (OverVector[ϵ])=(0,0,1), dephasing and relaxation of superposition state

Block[{ϵx,ϵy,ϵz,x,y,z,T1,T2,zs,ini,tini,tfinal},

ϵx[τ_]=0;ϵy[τ_]=0;ϵz[τ_]=3;

zs=-1; (* steady state *)
ini={x[0]==0,y[0]==1,z[0]==0};

T1=3;
T2=1.5;

tini=0.; tfinal=4π;
vec4[τ_]={x[τ],y[τ],z[τ]}/.First@NDSolve[{Sequence@@equs,Sequence@@ini},{x[τ],y[τ],z[τ]},{τ,tini,tfinal}];
Animate[Show[ParametricPlot3D[vec4[θ],{θ,0,a}],vectoronsphere@blochvecmarrowcart[vec4[a]]],{{a,$MachineEpsilon},tini,tfinal},AnimationRunning->False]
]

Steady state solutions at resonance and in rotating frame at the driving frequency

Block[{ϵx,ϵy,ϵz,x,y,z,T1,T2,zs,ini,tini,tfinal},

ϵx[τ_]=1;ϵy[τ_]=0;ϵz[τ_]=.0;

zs=-1; (* steady state *)
ini={x[0]==0,y[0]==0,z[0]==-1};

T1=1;
T2=1;

tini=0.; tfinal=π;
(* For small ϵx: *)
vec5s[τ_]={x[τ],y[τ],z[τ]}/.First@NDSolve[{Sequence@@equs,Sequence@@ini},{x[τ],y[τ],z[τ]},{τ,tini,tfinal}];

(* For medium ϵx: *)
ϵx[τ_]=2;
vec5m[τ_]={x[τ],y[τ],z[τ]}/.First@NDSolve[{Sequence@@equs,Sequence@@ini},{x[τ],y[τ],z[τ]},{τ,tini,tfinal}];

(* For large ϵx: *)
ϵx[τ_]=20;
vec5l[τ_]={x[τ],y[τ],z[τ]}/.First@NDSolve[{Sequence@@equs,Sequence@@ini},{x[τ],y[τ],z[τ]},{τ,tini,tfinal}];

Animate[Show[ParametricPlot3D[vec5s[θ],{θ,0,a+$MachineEpsilon}],vectoronsphere@blochvecmarrowcart[vec5s[a]],
ParametricPlot3D[vec5m[θ],{θ,0,a+$MachineEpsilon},PlotStyle->Directive[Yellow]],
blochvecmarrowcart[vec5m[a],color->Yellow],
ParametricPlot3D[vec5l[θ],{θ,0,a+$MachineEpsilon},PlotStyle->Directive[Green]],
blochvecmarrowcart[vec5l[a],color->Green]],{{a,100$MachineEpsilon},tini+100$MachineEpsilon,tfinal},AnimationRunning->False]
]

© Stefan Filipp, 10/2009