Teleportation, Harmonic Oscillator

Some Functions  and Definitions:

(* 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π];
{θ,ϕ}]

(* define a bloch sphere graph *)
With[{r=1},xy=ParametricPlot3D[{r Cos[θ],r Sin[θ],0},{θ,0,2 π},PlotStyle->Dashed,Boxed->False,Axes->False];
yz=ParametricPlot3D[{0,r Cos[θ],r Sin[θ]},{θ,0,2 π},PlotStyle->Dashed,Boxed->False,Axes->False];
xz=ParametricPlot3D[{r Cos[θ],0,r Sin[θ]},{θ,0,2 π},PlotStyle->Dashed,Boxed->False,Axes->False];
bloch=Show[Graphics3D[{Green,Opacity[.25],Sphere[{0,0,0},r],Opacity[1],Black,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[θ]}}]}]

(* draw vector for bloch sphere taking an arbitrary pure input state *)
blochvecstate[state_]:=blochvec@blochcoord[state]
(*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,θ,ϕ]
(* 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]
(* 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]}

id=(({{1, 0}, {0, 1}}));σx=(({{0, 1}, {1, 0}}));σy=(({{0, -I}, {I, 0}}));σz=(({{1, 0}, {0, -1}}));
tps[x_]:=Outer[Times,Conjugate[x],x]
tps[x_,y_]:=Outer[Times,Conjugate[x],y]
(* groundstate *) 
gs={1,0};
(* excited state *)
es={0,1};
(* superpositions state *)
xps={1,1}/(Sqrt[2]);
xms={1,-1}/(Sqrt[2]);
yps={1,I}/(Sqrt[2]);
yms={1,-I}/(Sqrt[2]);

TensorProduct[a_,b_]:=Table[(*a,b are nxn and mxm-martices*)a[[Ceiling[s/Length[b]],Ceiling[t/Length[b]]]]*b[[s-Floor[(s-1)/Length[b]]*Length[b],t-Floor[(t-1)/Length[b]]*Length[b]]],{s,1,Length[a]*Length[b]},{t,1,Length[a]*Length[b]}]

gg={1,0,0,0};ge={0,1,0,0};eg={0,0,1,0};ee={0,0,0,1};

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]

ZZ=TensorProduct[σz,σz];
RX[α_]:=Urot[α,1,π/2,0]
RY[α_]:=Urot[α,1,π/2,π/2]
RZ[α_]:=Urot[α,1,0,0]
RX1[α_]:=TensorProduct[Urot[α,1,π/2,0],(({{1, 0}, {0, 1}}))]
RY1[α_]:=TensorProduct[Urot[α,1,π/2,π/2],(({{1, 0}, {0, 1}}))]
RZ1[α_]:=TensorProduct[Urot[α,1,0,0],(({{1, 0}, {0, 1}}))]
RX2[α_]:=TensorProduct[(({{1, 0}, {0, 1}})),Urot[α,1,π/2,0]]
RY2[α_]:=TensorProduct[(({{1, 0}, {0, 1}})),Urot[α,1,π/2,π/2]]
RZ2[α_]:=TensorProduct[(({{1, 0}, {0, 1}})),Urot[α,1,0,0]]

Questions 3 :

CNOT=(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}}));

UZZ[α_]=MatrixExp[-I α ZZ];
UZZ[α]//MatrixForm

PG=-I (E^(I π/4))RZ2[-π/2].RZ1[-π/2].UZZ[π/4];
PG//MatrixForm

Its action on the states is given by

MatrixForm/@{PG.gg,PG.ge,PG.eg,PG.ee}

Instead of rotating  qubit 2 we rotate  qubit 1 :

RY1[-π/2].PG.RY1[π/2]//MatrixForm

A phase shift is missing :

MatrixForm/@(RY1[-π/2].PG.RY1[π/2].#&/@{gg,ge,eg,ee})

RZ2[π]//MatrixForm

I RZ2[π].RY1[-π/2].PG.RY1[π/2]//MatrixForm

CNOT21=I RZ2[π].RY1[-π/2].PG.RY1[π/2]

SWAP=CNOT.CNOT21.CNOT

initial=tps[gg]+tps[eg];
%//MatrixForm

SWAP.(gg+eg);
%
gg+ge

Teleportation

Antisymmetric Bell state |((Subscript[Ψ, 23])^(⟩-)) and state  |(Subscript[Ψ, 1])⟩  to be teleported :

Psi23m=(1/(Sqrt[2]))(eg-ge)
Psi1=a gs+b es

Tensor Product state |(SuperMinus[(Subscript[Ψ, 123])])⟩:

Psi123m=Flatten@Outer[Times,Psi1,(1/(Sqrt[2]))(eg-ge)]

Basis vectors :

(({{1}, {0}, {0}, {0}, {0}, {0}, {0}, {0}}))=|000>,(({{0}, {1}, {0}, {0}, {0}, {0}, {0}, {0}}))=|001>,(({{0}, {0}, {1}, {0}, {0}, {0}, {0}, {0}}))=|010>,(({{0}, {0}, {0}, {1}, {0}, {0}, {0}, {0}}))=|011>,(({{0}, {0}, {0}, {0}, {1}, {0}, {0}, {0}}))=|100>,(({{0}, {0}, {0}, {0}, {0}, {1}, {0}, {0}}))=|101>,(({{0}, {0}, {0}, {0}, {0}, {0}, {1}, {0}}))=|110>,(({{0}, {0}, {0}, {0}, {0}, {0}, {0}, {1}}))=|111>

basisstring={"ggg","gge","geg","gee","egg","ege","eeg","eee"}
basis=Table[RotateRight[{1,0,0,0,0,0,0,0},j],{j,0,7}]

MatrixForm/@basis

Projector onto |(SuperMinus[(Subscript[Ψ, 12])])⟩

tps[(1/(Sqrt[2]))(ge-eg)]

Include also Hilbert space of particle 3 :

Ppsi12m=TensorProduct[tps[(1/(Sqrt[2]))(ge-eg)],(({{1, 0}, {0, 1}}))]

Make a Bell-measurement of  |(SuperMinus[(Subscript[Ψ, 12])])⟩⟨ (SuperMinus[(Subscript[Ψ, 12])])| on the state|(SuperMinus[(Subscript[Ψ, 123])])⟩ :

psiproj=FullSimplify[((Ppsi12m.Psi123m)/(Sqrt[Psi123m\[Conjugate].Ppsi12m.Psi123m])),Assumptions->{a∈Reals,b∈Reals}]/.(b^2)->(1-(a^2))
basisstring

psiproj.basisstring//FullSimplify

Decomposing the state  again into a tensorproduct:

(psiproj.#)&/@basis[[{2,4,6,8}]]
(psiproj.#)&/@basis[[{1,3,5,7}]]

-> a ((| ge - eg ⟩)/(Sqrt[2])) | g⟩ +b ((| ge - eg⟩)/(Sqrt[2])) | e ⟩ = ((| ge - eg ⟩)/(Sqrt[2]))(a | g⟩ +b | e ⟩ )

i.e. the state can be written as (SuperMinus[(Subscript[Ψ, 12])])⊗(a |g⟩+b|e⟩)

Numerical Eigenfunctions and Eigenvalues of Harmonic Oscillator Hamiltonian

How to Solve Schrödinger's Equation Numerically

We can discretize the potential, and the wavefunction in space, over a suitable range.

ψ[x] becomes an array {ψ[(Subscript[x, 1])],ψ[(Subscript[x, 2])],...}
Likewise for V[x]

The step size is ((Subscript[x, i+1])-(Subscript[x, i]))=dx

Now we can make a first order approximation of the second derivative of the wavefunction:

ψ''[x] becomes ψ''[(Subscript[x, i])] = (1/(dx^2))(ψ[(Subscript[x, i+1])]-2ψ[(Subscript[x, i])]+ψ[(Subscript[x, i-1])])

Now we can write the Schrödinger equation as a matrix eigenvalue/eigenvector problem:

H ψ = E ψ

writing the hamiltonian operator -((ℏ^2)/(2m))((Subscript[∇, i])^2)+V[(Subscript[x, i])] now as a matrix.

This matrix now looks like this:

ExampleMatrix[size_]:=(DiagonalMatrix[Table[-((ℏ^2)/(2m))((-2)/(dx^2))+V["x"<>ToString[j]],{j,size}]]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{size-1}],1]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{size-1}],-1])

ExampleMatrix[5]//MatrixForm

DiffMatrix[m_,dx_,Vx_]:=(DiagonalMatrix[-((ℏ^2)/(2m))((-2)/(dx^2))+Vx]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{Length[Vx]-1}],1]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{Length[Vx]-1}],-1])

Example: A harmonic oscillator (by courtesy of P. Leek)

Vfunc[V_,x_]:=-V+30 (x^2)

Varray[V_,boundary_,dx_]:=Table[Vfunc[V,x],{x,-boundary,boundary,dx}];

For V = 100 and x from 0 to 5 with 500 divisions the potential looks like

harmpotential=Varray[100,5,2 5/500];

ListPlot[harmpotential,Filling->1.1Min[harmpotential],Frame->True,PlotRange->All]

Set up the Matrix for the Problem :

MProblem=DiffMatrix[1,2 5/500.,harmpotential];

Calulate the Eigenfunctions and Eigenvectors:

answers=Eigensystem[MProblem];

{evals,evects}=Sort[Transpose[answers],#1[[1]]<#2[[1]]&]//Transpose;

(* From now on we will just arbitrarily take the lowest 10 levels *)

stateenergies=Take[evals,10]

ListPlot[ Table[evects[[i]]+(i-1)/2,{i,1,10}],Joined->True,Filling->Table[n->(n-1)/2,{n,1,10}],PlotRange->All,ImageSize->600,Frame->True]

potplot=ListPlot[harmpotential,Filling->1.1Min[harmpotential],Joined->True,Frame->True,PlotRange->All,GridLines->{None,stateenergies},GridLinesStyle->Directive[Thickness[0.002],Red],Axes->False]

Plot the Eigenvectors on top of the Eigenenergies :

ScaleFactor=20;
stateplot=ListPlot[Table[ScaleFactor evects[[i]]+evals[[i]],{i,1,10}],Joined->True,Filling->Table[n->evals[[n]],{n,1,10}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->600,Frame->True,Axes->False]

Show[stateplot,potplot,PlotRange->{Automatic,{-100,0.1 100}}]

(* Very nearly perfectly equally spaced energy levels *)
ListPlot[tab=Table[stateenergies[[i+1]]-stateenergies[[i]],{i,1,10-1}],PlotRange->{Automatic,{7.5,8}},PlotStyle->PointSize[0.02],GridLines->{None,{Mean[tab]}},Frame->True,Axes->False]

(* Slight deviation due to numerical errors *)

Frequency response of a damped Harmonic Oscillator

Questions 4