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

Question 1 :

rho[r1_,r2_,r3_]:=1/2((({{1, 0}, {0, 1}}))+r1 σx+r2 σy+r3 σz)
MatrixToBloch[x_]:={Tr[σx.x],Tr[σy.x],Tr[σz.x]}

a)

rho1a=rho[0,0,0.3]

rho1a//MatrixForm

rho1a\[ConjugateTranspose]==rho1a

Tr[rho1a]

Eigenvalues[rho1a]

{rho1a\[ConjugateTranspose]==rho1a,Tr[rho1a],Eigenvalues[rho1a]}

vectoronsphere@blochvecmcart@MatrixToBloch[rho1a]

b)

rho1b=rho[.9,.6,0.3]

Norm[{.9,.6,0.3}]

rho1b//MatrixForm

{rho1b\[ConjugateTranspose]==rho1b,Tr[rho1b],Eigenvalues[rho1b]}//Chop

not physical!

vectoronsphere@blochvecmcart@MatrixToBloch[rho1b]

c)

rho1c=rho[0,1,0]

rho1c//MatrixForm

tps[{1,I},{1,I}]

{rho1c\[ConjugateTranspose]==rho1c,Tr[rho1c],Eigenvalues[rho1c]}//Chop

vectoronsphere@blochvecmcart@MatrixToBloch[rho1c]

Question 2: What' s the difference between mixed and pure states?

2 a)

ρemit=3/4ρgs + 1/4 ρes

mixed state generation and measurement

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

1/2((({{1, 0}, {0, 1}}))+Cos[ϕ]σz + Sin[ϕ]σx)

P[ϕ]//ComplexExpand//Simplify

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

ψemit=(Sqrt[3/4])gs+(Sqrt[1/4])es;

pure state generation and measurement

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

2 b)

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

2 c)

Measurement operators along x, y and z

Px=1/2((({{1, 0}, {0, 1}}))+σx)
Py=1/2((({{1, 0}, {0, 1}}))+σy)
Pz=1/2((({{1, 0}, {0, 1}}))+σz)

Measurement outcomes along x, y and z

Tr[Px.rho[r1,r2,r3]]//Simplify
Tr[Py.rho[r1,r2,r3]]//Simplify
Tr[Pz.rho[r1,r2,r3]]//Simplify

Solve[{mx==Tr[Px.rho[r1,r2,r3]],
my==Tr[Py.rho[r1,r2,r3]],
mz==Tr[Pz.rho[r1,r2,r3]]},{r1,r2,r3}]

Question 3: Density Matrix for thermal state

A harmonic oscillator is desribed by the Hamiltonian H=ℏ (Subscript[ω, r])((SuperDagger[a]).a +1/2) = ℏ (Subscript[ω, r]) (n+1/2), where (a^((†)))denotes the operator annihilating (creating)  an energy quantum, and n denotes the operator for the number of energy quanta, that is n|0⟩=0, n|1⟩=1,n|2⟩=2,... 
Find the density operator that describes the state of the harmonic oscillator, which  is  in thermal equilibrium with a heat bath at temperature T.

<<PhysicalConstants`

hbar=PlanckConstant/(2π)/Joule/Second;
kb=BoltzmannConstant/Joule Kelvin;

P[n_,T_,ω_]:=(1-Exp[- ((ℏ ω)/(kB T))])Exp[-((n ℏ ω)/(kB T))]

P[0,300,(10^9)2π]/.{ℏ->hbar,kB->kb}

System at 100 mK and frequency of 1 GHz

probability to find system in the ground state :

P[0,0.1,(10^9)2π]/.{ℏ->hbar,kB->kb}

probability to find system in the first excited state :

P[1,0.1,(10^9)2π]/.{ℏ->hbar,kB->kb}

Planck Law (average number of excitations at a given temperature <n> = Σ n p(n) ) :

Sum[P[n,T,ω]n,{n,0,∞}]//Simplify

again for System at 100 mK and 1 GHz

Sum[P[n,0.1,(10^9)2π]n,{n,0,∞}]/.{ℏ->hbar,kB->kb}

Plot for different temperatures

probs=Table[Tooltip[P[n,temp,(10^9)2π]/.{ℏ->hbar,kB->kb},ToString[temp]],{temp,1,301,50}];

Plot[probs,{n,0,5000},PlotRange->{0,0.0005},Frame->True,FrameLabel->{"n","(Subscript[p, ω])(n)"},BaseStyle->{FontSize->14}]

Manipulate[Plot[P[n,temp,(10^9)2π]/.{ℏ->hbar,kB->kb},{n,0,5000},PlotRange->{0,0.0005},Frame->True,FrameLabel->{"n","(Subscript[p, ω])(n)"},BaseStyle->{FontSize->14}],{temp,0.0001,300}]

Mixed state :

ρthermal[T_,nmax_,ω_]:=DiagonalMatrix[Table[P[n,T,ω],{n,1,nmax}]]

Wrong Normalization (not adapted to size of Hilbert space taken into account)

ρthermal[T,3,ω]//MatrixForm

Tr[ρthermal[300,400,(10^9)2π]/.{ℏ->hbar,kB->kb}]

Sum[Exp[ ((-n ℏω)/(kB T))],{n,0,∞}]

ρthermal2[T_,nmax_,ω_]:=DiagonalMatrix[Table[Exp[ ((-n ℏω)/(kB T))],{n,0,nmax}]/Sum[Exp[ ((-n ℏω)/(kB T))],{n,0,nmax}]]

ρthermal2[300,100,(10^9)2π]/.{ℏ->hbar,kB->kb}//Tr

<<BarCharts`

temperature = 0.1 K, ν = 1 GHz

BarChart3D[ρthermal2[0.1,20.,(10^9)2π]/.{ℏ->hbar,kB->kb},BarStyle->{Opacity[0.7],Red},PlotRange->All]

Typical temperatures and frequencies in experiments with superconducting quantum circuits coupled to a harmonic oscillator (microwave resonator) are about 20mK and the typical frequencies are ν=6GHz:

BarChart3D[ρthermal2[0.02,20,6 (10^9)2π]/.{ℏ->hbar,kB->kb},BarStyle->{Opacity[0.7],Red},PlotRange->All]

Unitary Evolution of a mixed state

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]