A3 Quantum Mechanics and Mathematical Methods (updated 2023)

Definitions for plots

These style directives are used later on in order to produce nice plots

ts={FontFamily->"Times",FontSize->12,FontWeight->"Bold"};
ts2={FontFamily->"Times",FontSize->14,FontWeight->"Bold"};
tsB={FontFamily->"Times",FontSize->18,FontWeight->"Bold"};
tsBB={FontFamily->"Times",FontSize->24,FontWeight->"Bold"};

Operator representation

The following sections illustrate the representation of operators by matrices. Since we work with orthonormal bases, the matrix representation of an operator can be found via the scalar product.

Scalar product

Define the scalar product between two vectors.
Note that vector/matrix multiplication in Mathematica is denoted by a “.”

SP[vec1_,vec2_]:=(Conjugate[vec1]).(vec2)

Mathematica is capable of symbolic and numerical calculations. Let’s define three vectors in a four-dimensional space:

vec1=Table[u[k],{k,1,4}];
vec2=Table[v[k],{k,1,4}];
vec3=Table[w[k],{k,1,4}];

Inspect the components of vec1:

MatrixForm[vec1]

Calculate the scalar product of vec1 and vec2:

SP[vec1,vec2]

Calculate the scalar product of vec1 and vec1:

SP[vec1,vec1]
%//Simplify
%//FullSimplify

Check the axiom that SP is linear in the second argument

SP[vec1,a*vec2 + b*vec3] - (a*SP[vec1,vec2]+b*SP[vec1,vec3])
%//Simplify

Verify the axiom regarding the interchange of the two arguments in SP:

SP[vec1, vec2] - Conjugate[ SP[vec2,vec1]]

Verify SP[vec,vec] > 0 if vec !=  0

SP[vec1,vec1]
%//Simplify
%//FullSimplify

Numerical example:

vec1 = {I,0,1};
vec2 = {2,1,1};

SP[vec1,vec2]
SP[vec1,vec1]
SP[vec2,vec2]

Pauli matrices

Define the Pauli matrices:

sx

sx={{0,1},{1,0}};
sx//MatrixForm

sy = {{0,-I},{I,0}};
sy//MatrixForm

sz = {{1,0},{0,-1}};
sz//MatrixForm

dim is the dimension of our state space, and bz[k] is the canonical (orthonormal) basis:

dim = 2;
bz[1]={1,0};
bz[2]={0,1};

The basis vectors bz[1] and bz[2] are the eigenvectors of sz with eigenvalues 1 and -1, respectively:

sz.bz[1]-bz[1]

sz.bz[2]-(-bz[2])

For any orthonormal basis b[k], the matrix representation of a linear operator A with respect to that basis can be found via the formula

A[i,j] = SP[  b[i], A.b[j]   ] 

The Pauli matrices as defined above are representations of the corresponding operators in the basis bz[k] of eigenvectors of sz. Let’s verify this explicitly:

Table[SP[bz[i],sx.bz[j]],{i,1,dim},{j,1,dim}];
%-sx//MatrixForm

Table[SP[bz[i],sy.bz[j]],{i,1,dim},{j,1,dim}];
%-sy//MatrixForm

Table[SP[bz[i],sz.bz[j]],{i,1,dim},{j,1,dim}];
%-sz//MatrixForm

Next let’s choose a different basis, the basis of eigenstates corresponding to sx:

evx=Eigenvectors[sx];
bx[1]=Normalize[evx[[2]]]
bx[2]=Normalize[evx[[1]]]

Check that these are the correct eigenstates with eigenvalues 1 and -1:

sx.bx[1]-bx[1]

sx.bx[2]-(-bx[2])

Let’s work out the matrix representation of the Pauli matrices in the basis bx[k]:

sxBX=Table[SP[bx[i],sx.bx[j]],{i,1,dim},{j,1,dim}];
%//MatrixForm

syBX=Table[SP[bx[i],sy.bx[j]],{i,1,dim},{j,1,dim}];
%//MatrixForm

szBX=Table[SP[bx[i],sz.bx[j]],{i,1,dim},{j,1,dim}];
%//MatrixForm

All we have done is a change of basis between orthonormal bases. The basis change from bx[k] to bz[k] is done via the unitary matrix U:

U = Transpose[{bx[1],bx[2]}];
Udag = Transpose[Conjugate[U]];

U.Udag//MatrixForm

Check that the Pauli matrices in the basis bx[k] can be obtained from those in the basis bz[k] via the basis transformation U:

sxBX - Udag.sx.U
syBX - Udag.sy.U
szBX - Udag.sz.U

Random matrix

Build a random matrix of dimension 10x10

dimR=10;
Rmat = RandomReal[{-5,5},{dimR,dimR}] + I RandomReal[{-5,5},{dimR,dimR}] ;

Make Rmat hermitian:

Rmat = Rmat + Transpose[Conjugate[Rmat]];

Find the eigenvalues and eigenstates of this matrix:

ESrmat=Eigensystem[Rmat];

Eigenvalues:

EvalRmat=ESrmat[[1]]

Eigenstates:

EvecRmat=ESrmat[[2]];
EvecRmat=Normalize[#]&/@ EvecRmat;

Find the representation of Rmat in the basis of eigenstates. This is just a diagonal matrix with the eigenvalues on its diagonal:

RmatDiag = Table[SP[EvecRmat[[i]],Rmat.EvecRmat[[j]]],{i,1,dimR},{j,1,dimR}]//Chop;
RmatDiag//MatrixForm

Diagonal[RmatDiag]-EvalRmat//Chop

We can also find the diagonal form of Rmat using the unitary transform UR:

UR = Transpose[EvecRmat];
URdag = Transpose[Conjugate[UR]];

UR.URdag//Chop

URdag.Rmat.UR-RmatDiag//Chop//MatrixForm

Functions of operators

Sine of a 2x2 matrix

In Quantum Mechanics we often use the exponential of operators in different contexts. Let’s calculate the Sine of the following matrix:

mat={{2,-I},{I,2}};

Find the eigenvalues and eigenvectors of this matrix:

ESmat=Eigensystem[mat];

Evalmat=ESmat[[1]]

Evecmat=ESmat[[2]];
Evecmat=Normalize[#]&/@ Evecmat;

The unitary U allows us to write mat in the form mat = U.DiagonalMatrix(Evalmat).Udag:

U = Transpose[Evecmat];
Udag = Transpose[Conjugate[U]];

U.DiagonalMatrix[Evalmat].Udag-mat//Simplify

Functions of operators are defined via the power series of the function. By taking the unitaries U and Udag outside the sum, we can just apply the function to the eigenvalues:
(Note that Sin threads over lists, i.e., Sin[Evalmat] = {Sin[Evalmat[[1]]],Sin[Evalmat[[2]]]})

SQmat=U.DiagonalMatrix[Sin[Evalmat]].Udag//Simplify

In Mathematica, you can calculate functions of matrices with “MatrixFunction”. Compare this with our result:

SQmat-MatrixFunction[Sin,mat]//Simplify

Time evolution

Build a random matrix of dimension 10x10

dimR=10;
Rmat = RandomReal[{-5,5},{dimR,dimR}] + I RandomReal[{-5,5},{dimR,dimR}] ;

Make Rmat hermitian:

Rmat = Rmat + Transpose[Conjugate[Rmat]];

Build an initial state:

Rvec = Table[0,{k,1,dimR}];
Rvec[[1]]=1;
Rvec//MatrixForm

Calculate time evolution via MatrixExp (time evolution operator):

Calculate the probability to find the system in the initial state at as a function of time:

TimeEvolution[t_,init_]:=MatrixExp[-I Rmat*t,init]

The variable TimeVal contains a list with time steps from 0 to 1:

TimeVal=Table[k,{k,0,1,0.001}];

RvecTime contains the evolved state at these time steps:

RvecTime=Table[TimeEvolution[TimeVal[[k]],Rvec],{k,1,Length[TimeVal]}];

Finally, we calculate the probability to find the system in the initial state at these time steps:

ProbTime=Table[{TimeVal[[k]],((Abs[SP[Rvec,RvecTime[[k]]]])^2)},{k,1,Length[TimeVal]}];

Plot the resulting function:

MexpPlot=ListPlot[ProbTime,BaseStyle->ts,Axes->False,Frame->True,PlotRange->All,Joined->True,FrameLabel->{"Time","Probability"},PlotStyle->{Red,Thick}]

Calculate time evolution using stationary states:

Find the eigenvalues and eigenstates of this matrix:

ESrmat=Eigensystem[Rmat];

Eigenvalues:

EvalRmat=ESrmat[[1]]

Eigenstates:

EvecRmat=ESrmat[[2]];
EvecRmat=Normalize[#]&/@ EvecRmat;

Expand the initial state Rvec (see above) in terms of eigenstates:

ExpansionCoeff= Table[SP[EvecRmat[[k]],Rvec],{k,1,dimR}];
SP[ExpansionCoeff,ExpansionCoeff]

Sum[EvecRmat[[k]]*ExpansionCoeff[[k]],{k,1,dimR}]-Rvec//Chop

Build the time-evolved state:

TimeEvolution2[t_]:=Sum[EvecRmat[[k]]*ExpansionCoeff[[k]]*Exp[-I EvalRmat[[k]]*t],{k,1,dimR}]

RvecTime2=Table[TimeEvolution2[TimeVal[[k]]],{k,1,Length[TimeVal]}];

ProbTime2=Table[{TimeVal[[k]],((Abs[SP[Rvec,RvecTime2[[k]]]])^2)},{k,1,Length[TimeVal]}];

StatPlot=ListPlot[ProbTime2,BaseStyle->ts,Axes->False,Frame->True,PlotRange->All,Joined->True,FrameLabel->{"Time","Probability"},PlotStyle->{Black,Dashed}]

Compare our result to the one obtained with the Matrix Exponential:

Show[MexpPlot,StatPlot]