Basic Intro to Mathematica

2 parts :

frontend (Notebook) - Interface to the user
kernel - process doing the calculations (plattform independent)

additionally:

Mathlink - Interface between Frontend and Kernel
Packages - Extensions to the built-in functions

Tips on using Mathematica

Run pieces of code by pressing Shift+Enter with the cursor on the code.

Use curly brackets to define lists of anything:    {1,2,3,4}      {"One","Two","Three"}      {(x^2),12,"Quantum"}

Use square brackets when calling and defining functions

Use the esc key for special characters eg esc-p-esc gives π

Variables

myvar=5

Functions

f[x_]=x^2

f[3]

Attention with function definitions :

f[x]=x+3

f[3]

f[x]

Immediate/delayed assignment

(* Immediate assignment *)
myrand1[x_]=x Random[];

(* Delayed assignment *)
myrand2[x_]:=x Random[];

Table[myrand1[3],{x,1,10}]
Table[myrand2[3],{x,1,10}]

Replacement rules

x

x /. x -> y

x

<< PhysicalConstants`

SpeedOfLight

kb T

kb T /. {kb -> BoltzmannConstant, T -> 300 Kelvin}

kb T

Lists - Everything is a List!

primeList=Table[Prime[n],{n,1,200}]

mylist=Table[{x,(x^2)},{x,0,10}]

Plots

analytplot=Plot[(x^2),{x,0,10},AxesLabel->{"x","(x^2)"}]

numplot=ListPlot[mylist,AxesLabel->{"x","(x^2)"},PlotRange->All]

Show[numplot,analytplot]

plotParams={Frame->True,FrameStyle->Directive[16,Black],FrameLabel->{"Which prime","The prime","The first 200 prime numbers",""},ImageSize->450};

primePlot=ListPlot[primeList,plotParams]

(* Fill in your local file path here *)
fullPath="D:\\";
Export[fullPath<>"primePlot.pdf",primePlot]

Help system

Use the Help system!!
Just press F1 on a function name, and the correct Help page will appear.

E.g.

Integrate
Sum
Fourier
Expand

Some linear algebra for QIPC

State vectors

A (pure) quantum state is represented by a normalized state vector |ψ⟩ ∈ H, where H denotes the Hilbert space (complex vector space).
As long as there is a finite number of states, or as long as we choose to only consider a finite number of states, in the system we se a list for representing the quantum state.
E.g.:

(* groundstate *) 
gs={1,0};
(* excited state *)
es={0,1};
(* superpositions state *)
xps={1,1}/(Sqrt[2]);
xms={1,-1}/(Sqrt[2]);

some function 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],θ,ϕ},
θ=2Chop@ArcCos[normx[[1]]Exp[-I arg1]];
ϕ=Mod[arg2-arg1,2π];
{θ,ϕ}]

(* arbitrary state {Cos[θ/2], (e^(I ϕ))Sin[θ/2]}

State[π/10,π/4]
blochcoord[%]//FullSimplify

A basis is composed of a maximal set of linearly independent vectors, e.g. the set {gs,es}

Solve[{c1 gs+c2 es=={0,0}},{c1,c2}]

Solve[{c1 xps+c2 xms=={0,0}},{c1,c2}]

The norm of a vector is simply

Norm[gs]
(Sqrt[gs.gs])

Bloch Sphere Representation :

(* 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["|0⟩",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,0,1.3 r}],Text[Style["|1⟩",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,0,-1.3 r}],Text[Style["|+⟩",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{1.3 r,0,0}],Text[Style["|-⟩",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{-1.3 r,0,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]

ground + excited state :

GraphicsRow[{Show[bloch,blochvec[{0,0},color->Red],SphericalRegion->True],Show[bloch,blochvec[{π,0},color->Blue],SphericalRegion->True]},ImageSize->700]

Superposition States

GraphicsRow[{Show[bloch,blochvec[{π/2,0},color->Blue],blochvec[{π/2,π/2},color->Red],SphericalRegion->True,PlotLabel->Style["((|0⟩+ | 1⟩))/(Sqrt[2])) and ((|0⟩ + I | 1⟩))/(Sqrt[2]))",FontSize->14]],Show[bloch,blochvec[{π/3,π/4},color->Blue],SphericalRegion->True,PlotLabel->Style["Cos (π/3)|0⟩+Sin[(π/3)](e^((Iπ)/4))|1⟩",FontSize->14]]},ImageSize->700]

Operators and States as Matrices

Matrices form a representation of operators acting on the states, which map vectors into vectors.

E.g the map from the ground state {1,0} to the excited state {0, 1} and vice versa is given by the matrix

piflip=(({{0, 1}, {1, 0}}));
piflip.gs
piflip.es

There are some canonical types of matrices that occur frequently in Quantum Mechanics:

Hermitian Matrix (Self-adjoint operator)

Definition :

(SuperDagger[A])==A, (SuperDagger[where ])denotes the conjugate transpose.

herm=(({{1, -0.3I}, {0.3I, -1}}));

herm\[ConjugateTranspose]==herm

They can be formed by linear combinations of the Pauli Matrices + Identity(with real coefficients):

id=(({{1, 0}, {0, 1}}));σx=(({{0, 1}, {1, 0}}));σy=(({{0, -I}, {I, 0}}));σz=(({{1, 0}, {0, -1}}));

herm2=c2 σy+c3 σz
herm2==herm2\[ConjugateTranspose]
herm2==herm2\[ConjugateTranspose]//FullSimplify[#,Assumptions->{c2∈Reals,c3∈Reals}]&

Properties :

*) Real Eigenvalues

Eigenvalues[herm2]

*) Orthogonal Eigenstates

{e1,e2}=Eigenvectors@herm2;
{ne1,ne2}={e1/Norm[e1],e2/Norm[e2]}

ne1.ne2//FullSimplify
FullSimplify[
Conjugate[ne1].ne2,Assumptions->{c1 ∈Reals,c2∈Reals,c3∈Reals}]

*) Normalized eigenvectors form a complete orthonormal basis set

Solve[d1 ne1+d2 ne2==0,{d1,d2}]

Applications:

*) Measurements in Quantum Mechanics are represented as Hermitian operators - property of real eigenvalues corresponds to observable quantities.
*) Hamiltonian of a system is Hermitian, since only Hermitian matrices produce unitary evolutions (as generators of the group of unitary matrices - Stones theorem).

Unitary Matrix

Transform one state vector into another while preserving its norm.

Definition:

U (SuperDagger[U])==(SuperDagger[U])U==1.

Property :

Eigenvalues with modulus 1

unitary[θ_]=(({{Cos[θ], I Sin[θ]}, {I Sin[θ], Cos[θ]}}))

unitary[θ]\[ConjugateTranspose].unitary[θ]//FullSimplify[#,Assumptions->{θ∈Reals}]&

TrigToExp@Eigenvalues@unitary[θ]

Applications:

*) Describe evolution of a quantum system.

Animate[Show[bloch,blochvec[blochcoord[unitary[θ].gs],color->Blue],SphericalRegion->True,PlotLabel->Style["U[θ].",FontSize->14]],{θ,0,2π},AnimationRate->0.1]

Examples: Hadamard Gate

HG=(1/(Sqrt[2]))(({{1, 1}, {1, -1}}));

Is it unitary?

HG\[ConjugateTranspose].HG//MatrixForm

What does it do to the ground state?

Show[bloch,blochvec@blochcoord[HG.gs],blochvec[blochcoord[gs],color->Red]]

What does it do to the excited state?

Show[bloch,blochvec@blochcoord[HG.es],blochvec[blochcoord[es],color->Red]]

Examples: bit flip Gate (σx)

bf=(({{0, 1}, {1, 0}}));
bf\[ConjugateTranspose].bf//MatrixForm

GraphicsRow[{Show[bloch,blochvec@blochcoord[bf.gs],blochvec[blochcoord[gs],color->Red]],Show[bloch,blochvec@blochcoord[bf.es],blochvec[blochcoord[es],color->Red]]},ImageSize->600]

Examples: bit-phase flip Gate ((Subscript[σ, y]))

bpf=(({{0, -I}, {I, 0}}));
bpf\[ConjugateTranspose].bpf//MatrixForm

GraphicsRow[{Show[bloch,blochvec@blochcoord[bpf.gs],blochvec[blochcoord[gs],color->Red]],Show[bloch,blochvec@blochcoord[bpf.es],blochvec[blochcoord[es],color->Red]]},ImageSize->600]

GraphicsRow[{Show[bloch,blochvec@blochcoord[bpf.xps],blochvec[blochcoord[xps],color->Red]],
Show[bloch,blochvec@blochcoord[bpf.xms],blochvec[blochcoord[xms],color->Red]]},ImageSize->600]

Examples: phase flip Gate ((Subscript[σ, z]))

pf=(({{1, 0}, {0, -1}}));
pf\[ConjugateTranspose].pf//MatrixForm

GraphicsRow[{Show[bloch,blochvec@blochcoord[pf.gs],blochvec[blochcoord[gs],color->Red]],Show[bloch,blochvec@blochcoord[pf.es],blochvec[blochcoord[es],color->Red]]},ImageSize->600]

GraphicsRow[{Show[bloch,blochvec@blochcoord[pf.xps],blochvec[blochcoord[xps],color->Red]],
Show[bloch,blochvec@blochcoord[pf.xms],blochvec[blochcoord[xms],color->Red]]},ImageSize->600]

Projectors

Projection operators are simply formed by | ψ⟩⟨ψ|, i.e.

(* generate a projector matrix from a pure state *)
tps[x_]:=Outer[Times,Conjugate[x],x]
tps[x_,y_]:=Outer[Times,Conjugate[x],y]

tps[gs]

Properties :

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

Applications :

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.

tps[gs].gs
tps[gs].es
Conjugate[State[π/3,π/4]].tps[gs].State[π/3,π/4]

((Cos[π/3/2])^2)

Show[bloch,blochvec[{π/3,π/4}],Graphics3D[{Red,Thickness[0.015],Line[{{0,0,0},{0,0,Cos[π/3]}}],Green,Line[{{Sin[π/3]Cos[π/4],Sin[π/3]Sin[π/4],Cos[π/3]},{0,0,Cos[π/3]}}]}]]

Density Matrix

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.

ρarb=0.5tps[gs,gs]+0.3 tps[State[π/10,0],State[π/10,0]]+0.2 tps[State[π/10,π],State[π/10,π]]

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]

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

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]

Questions

Density Matrix for thermal state

A harmonic oscillator is desribed by the Hamiltonian H=(Subscript[ω, r])(SuperDagger[a]).a = (Subscript[ω, r]) n, 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.

Decomposition of a density matrix

An experimental setup produces  spin-1/2 particles polarized either in the state | ↑⟩  or in the state |↓⟩ with equal probability. A different experimental setup produces also spin-1/2 particles but this time polarized either in the state |+⟩  =(| ↑⟩  +| ↑⟩)/(Sqrt[2]) or in the state |-⟩  = (| ↑⟩ - | ↑⟩)/(Sqrt[2]). Can you find a measurement which can distinguish between these two setups?