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

Projection operators:

A test whether the system is in the ground state yields either yes (0) or no (1). The measurement corresponding to this test is represented by the operator | 0⟩ ⟨0 | . The complementary operator is  | 1⟩ ⟨1 |, which tests whether the qubit is in its excited state or not.

in another notation
 |±><±|=(Subscript[P, z±])=1/2(1±σz):

In matrix representation :

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

In Quantum mechanics the question "Is the system in the groundstate?" corresponds to applying the projector  | 0⟩ ⟨0 | . After this application the system is in the state   ((P|ψ⟩)/(p(0)))=((|0⟩)/(Sqrt[⟨ψ | P | ψ⟩]))=((|0⟩)/(Sqrt[⟨ψ | 0⟩ ⟨0 | ψ⟩])).  The expectation value ⟨ψ | 0⟩ ⟨0 | ψ⟩ tells us what the probability is to find the system in the ground state.

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

Is the system in the ground state?

Answer : "yes", if system is in the ground state

Pzm.(({{0}, {1}}))//MatrixForm

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

|0><0| applied onto the ground state |0> has expectation value 1, i.e. the system is found with certainty in state (({{0}, {1}})).

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

Answer "no", if system is in excited state

Pzm.(({{1}, {0}}))//MatrixForm
(({{0, 1}})).Pzm.(({{1}, {0}}))

|0><0| applied onto the excited state |1>  has expectation value 0, i.e. the system is never found in state (({{0}, {1}})).

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

Answer "maybe", if system is in superposition state

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

has both components in z + and z - direction

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

state ofter projective measurement onto the positive or negative z-axis :

Pzp.xpstate//MatrixForm
Pzm.xpstate//MatrixForm

probability to find system in either of the basis states :

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

Measurement along arbitrary axis

Is the system in the state a |0⟩+ b |1⟩? This test is implemented by the operator (a|0⟩+ b|1⟩)(a*⟨0|+ b*⟨1|) and corresponds to measure along the axis defined by the coefficients a and b.

In another notation this operator can be written as

(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})

Pp[θ,ϕ]//MatrixForm
Pm[θ,ϕ]//MatrixForm

Property of projection (measurement) operators: (P^2)=P

(asking a question a second time does not give you more information)

Pp[θ,ϕ].Pp[θ,ϕ]==Pp[θ,ϕ]//Simplify
Pm[θ,ϕ].Pm[θ,ϕ]==Pm[θ,ϕ]//Simplify

Probability to get outcome + 1 when measuring the state | 0⟩ along an axis in the z-y plane (θ,ϕ=0):

(({{0, 1}})).Pp[θ,0].(({{0}, {1}}))//FullSimplify

In optics this is called  Malus' law.

Plot[(({{0, 1}})).Pp[θ,0].(({{0}, {1}})),{θ,0,2π},Frame->True,FrameLabel->{θ,"<(SubPlus[P])>"},BaseStyle->{FontSize->14}]

Mixed States

How to describe statistical mixtures of pure states in Quantum Mechanics? 

Formally, this is not possible with state vectors, since a | ψ⟩+ b | ϕ⟩ is a superposition state and leads to quantum interference phenomena due to the well defined phase relation, i.e. coherence between states. Instead one introduces the notion of density matrices (operators). I. e. an incoherent mixture of pure states forms a mixed state.

ρ =(∑_i)(Subscript[r, i])|(Subscript[ψ, i])⟩⟨(Subscript[ψ, i])|

where the state have to be normalized (the probability to find the system in state | (Subscript[ψ, i])⟩ ⟨(Subscript[ψ, i]) | must not exceed 1), i.e. (∑_i)(Subscript[r, i])=1 and (Subscript[r, i])≥0.

(* groundstate *) 
gs={1,0};
(* excited state *)
es={0,1};
(* 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]

ρarb=0.5tps[gs,gs]+0.3 tps[State[π/10,0],State[π/10,0]]+0.2 tps[State[π/10,π],State[π/10,π]]
ρ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]).

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

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

Applications :

Describe the state of a quantum system.

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.

How do we determine (OverVector[r])?
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]}

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

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

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

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

another incoherent superposition from above:

vectoronsphere@blochvecmcart@MatrixToBloch[ρarb]

arbitrary mixed state :

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

Measurement :

Tr[ρ.M] with ρ the density matrix and M the Hermitian measurement operator

Tr[tps[gs].tps[State[π/3,π/4]]]

State Evolution :

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

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

What' s the difference between mixed and pure states?

ρemit=0.7ρgs + 0.3 ρes

mixed state generation and measurement

P[ϕ_]:=tps[State[ϕ,0]]

Plot[Tr[P[ϕ].ρemit],{ϕ,0,2π},Frame->True,FrameLabel->{"ϕ[rad]","⟨(Subscript[P, z↑])⟩"},BaseStyle->{FontSize->14}]

ψemit=(Sqrt[0.7])gs+(Sqrt[0.3])es;

pure state generation and measurement

Plot[Tr[P[ϕ].tps[ψemit]],{ϕ,0,2π},Frame->True,FrameLabel->{"ϕ[rad]","⟨(Subscript[P, z↑])⟩"},BaseStyle->{FontSize->14}]

Decomposition freedom :

There are many ways to obtain the same density matrix from different statistical ensembles.
E.g.

ρ1=3/4tps[gs]+1/4 tps[es]
ψa=(Sqrt[3/4])gs+(Sqrt[1/4])es;
ψb=(Sqrt[3/4])gs-(Sqrt[1/4])es;
ρ2=1/2 tps[ψa]+1/2 tps[ψb]