Unprotect[scalar, hermitian];
Needs["Utilities`CleanSlate`"]
CleanSlate[]
Remove["Global`*"]
Protect[scalar, hermitian]

Last modified on August 10, 2015

Section 1: Vectors in Mathematica

Mathematica represents vectors as an ordered list of numbers (real or complex) called the coordinates or the components of the vector.  In quantum mechanics we use bra and ket notation.  I do not use this notation here.

Cell 1. Mathematica represents vectors as lists

Below I use Table to generate a 5-dimensional vector with real components. The function Random[Real,{-10,10}] generates random real numbers between -10 and 10:

I generate below several vectors which I use in the examples that follow

Two vectors having numerical components:

SeedRandom[4]; (* this function ensures that I get the same random numbers every time I generate the vectors*)
vector1=Table[RandomReal[{-10,10}],{i,1,5}]
vector2=Table[RandomReal[{-10,10}],{i,1,5}]

Use of SeedRandom[4] ensures that I get the same vectors whenever I run the program (and so do you)

I also generate two five-dimensional vectors whose components are symbols

ve1=Table[a[i],{i,1,5}]
ve2=Table[b[i],{i,1,5}]

Mathematicians represent sometimes vectors as a column. Mathematica can do that also by using the function MatrixForm.

www= ve3//MatrixForm

The function MatrixForm is used for display only.  www define above is not longer a vector, as you can see below:

www//FullForm

The components of the vectors in quantum mechanics are often complex numbers. Mathematica has no trouble generating them. Generate several vectors having complex components.

SeedRandom[3];
vecC1=Table[RandomComplex[{-10-I 10,10+I 10}],{i,1,4}]
vecC2=Table[RandomComplex[{-5-I 3,4+I 2}],{i,1,4}]

These lists of complex numbers can also be displayed as a column by using the function MatrixForm

vecC1//MatrixForm

Cell 2. Operations with vectors: addition of two vectros and the multiplication of a vector with a number

Addition

Display the vectros ve1 and ve2 and then add them. You can use the operator + for addition

ve1
ve2
ve1+ve2

You can use the same syntax for vectors that have numerical components

vector1
vector2
vector1+vector2

Multiplication of a vector with a number (real or complex)

The syntax is  number*vector of number vector

Mathematica multiplies each component with the number:

3 ve1

3*vecC1

Cell 3:  Scalar product. If the components of the vector are real numbers use the function Dot.  If they are complex numbers write your own function to perform the scalar product

ve1
ve2
Dot[ve1,ve2]

There is a short-hand notation for the function Dot:
	 v.w is the same as Dot[v, w]

ve1.ve2
Dot[ve1,ve2]

If you perform this operation with vectors having numerical components you get a number

vector1 .vector2

The definition of the scalar product for vectors whose components are complex numbers is given below:
	if 
		e1={a,b,c,..}						(1)
		e2= {A, B, C,..}					(2)
and any of a, b, c, A, B, C are complex numbers then the scalar product is defined by

	⟨e1|e2⟩ ≡ a* A + b* B + c* C + ...				(3)
	
here a* is the complex conjugate of a (complex conjugation means replacing i =(Sqrt[-1]) with -i ).
The function Dot does not give the correct scalar product if the component of the vector are complex numbers. 
We  have to create our own function which I call scalar.

To calculate the scalar product for vectors with complex components use:

Unprotect[scalar];
Clear[scalar]
scalar[x_?VectorQ,y_?VectorQ]:= 
Conjugate[x].y
Protect[scalar];

Test that the function scalar works correctly

Clear[a,A,b,B]
scalar[{a+ I A,c},{d,b+I B}]

Note that Mathematica does not know whether A or B are real or complex numbers. This is why the function Conjugate cannot be evaluated.

Calculate the scalar product for the vectors vecC1 and vecC2 whose components are complex numbers

vecC1
vecC2
scalar[vecC1,vecC2]

Note that Dot gives an incorrect result:

Dot[vecC1,vecC2]

The evaluation below is also incorrect:

Dot[{a+ I A,5},{2,b+I B}]

?Dot

Cell 4: The length of a vector a={(Subscript[a, 1]), (Subscript[a, 2]), ...,(Subscript[a, N])} is defined by

	(Sqrt[(∑_(i=1)^N)(Subscript[a, i]^( * ))(Subscript[a, i])])≡ (Sqrt[⟨a|a⟩])≡ ||a||

This is also called the norm of the vector a.

Calculate the norm of the vector ve1

ve1
Sqrt[scalar[ve1,ve1]]

The function scalar did not evaluate complex conjugation because the Mathematica does not know whether the symbols a[i] are real or complex numbers. Here is the norm of two vectors with numberical components

vector1
vector2
Sqrt[scalar[vector1,vector1]]

Mathematica provides the function Norm which works correctly with vectors having complex coordinates (components).

ve1
Norm[ve1]

Norm assumes that the symbolic components are complex variables.  I do the a calculation with Norm and the same calculation using our function scalar to check that they give the same result when the vectors have complex components:

{Norm[vecC1],Sqrt[scalar[vecC1,vecC1]]//Chop}

Cell 5: To normalize a vector you divide it by its length.

If the components of the vector are real numbers you can use the function Dot to calculate the scalar product (in the length of the vector).  A vector v is normalized by using

	 v/(Sqrt[v.v])≡v/||v||
Since (Sqrt[v.v]) is the length for the vector the normalized vector is a vector of unit length parallel to the original vector v.

Normalize the vector vector1 which has numerical components. Because the components of vector1 are real numbers I use the function Dot (the shorthand notation) to calculate the scalar product

vector1
normalizedVector1=vector1/Sqrt[vector1.vector1]

Check that the normalized vector has length 1:

Norm[normalizedVector1]

To normalize a vector whose coordinates are complex  use Norm

vCC=vecC1/Norm[vecC1]

Check that the vector is normalized.

scalar[vCC,vCC]//Chop

Norm[vCC]

Or use my the function scalar which I defined earlier. Below temp is the normalized vector vecC1 and the second line tests that the result is normalized

Cell 6: Basis set.

You have learned in your mathematics class that a vector can always be written as a linear combination of basis set vectors. The most common basis set consists of vectors of unit length that are orthogonal to each other.  By definition, two vectors are orthogonal if their scalar product is zero.  

When I write a vector like
	ve1 = {a[1], a[2], a[3], a[4], a[5]}
what is the basis set?  In other words if I want to write this as 

	ve1= a[1]  (Subscript[(OverVector[e]), 1])+a[2]  (Subscript[(OverVector[e]), 2])+a[3]  (Subscript[(OverVector[e]), 3])+a[4]  (Subscript[(OverVector[e]), 4])+a[5]  (Subscript[(OverVector[e]), 5]),
	
where  (Subscript[(OverVector[e]), 1]), (Subscript[(OverVector[e]), 2]), (Subscript[(OverVector[e]), 3]),  (Subscript[(OverVector[e]), 4]) and (Subscript[(OverVector[e]), 5]) are orthonormal vectors, what do I use for   (Subscript[(OverVector[e]), 1]), (Subscript[(OverVector[e]), 2]), (Subscript[(OverVector[e]), 3]),  (Subscript[(OverVector[e]), 4]) and (Subscript[(OverVector[e]), 5]) ?  I will show that the basis set consists of the vectors
	(Subscript[(OverVector[e]), 1])={1, 0, 0, 0, 0}
	(Subscript[(OverVector[e]), 2])={0, 1, 0, 0, 0}
	.
	.
	(Subscript[(OverVector[e]), 5])={0, 0, 0, 0, 1}

I can generate all these vectors in one line by using the function DiagonalMatrix

Unprotect[e];
Clear[e]
{e[1], e[2], e[3], e[4], e[5]}=DiagonalMatrix[{1,1,1,1,1}] ;
Protect[e];

Check that the assignment works as expected. The second vector in the basis set is

e[2]

I can now evaluate the expression

	a[1]  (Subscript[(OverVector[e]), 1])+a[2]  (Subscript[(OverVector[e]), 2])+a[3]  (Subscript[(OverVector[e]), 3])+a[4]  (Subscript[(OverVector[e]), 4])+a[5]  (Subscript[(OverVector[e]), 5])

For example calculate a[1] e[1]

a[1] e[1]

Similarly

a[2] e[2]

Now let us show that the sum

	(∑_(i=1)^5)a[i] (Subscript[(OverVector[e]), i]) 
is equal to the vector {a[1], a[2],...,a[5]}

Plus@@Table[a[i]e[i],{i,1,5}]

The syntax used above is very common and useful.  Here is what Table does

temp=Table[a[i]e[i],{i,1,5}]

Plus@@ applied to this expression adds the elements

Plus@@temp

This means that the expression 

	(∑_(i=1)^5)a[i] (Subscript[(OverVector[e]), i]) 
is the same as the vector ve1

ve1

Section 2: Matrices  in Mathematica: representation and simple operations

Cell 1. Mathematica represents a matrix as lists of lists. Here are some examples:

m=.
Clear[m1]
m1=Table[m[i][j],{i,1,5},{j,1,5}]

In mathematics we use the notation (Subscript[m, ij]) for the elements of a matrix. The symbols m[i][j] used above are an imitation of that. The order of the elements is important. If two matrices have the same elements but in a different order they are different matrices. The matrix {{2,1},{3,4}} is not the same as the matrix {{1,2},{3,4}}

Mathematicians display a small matrix as a table. You can do that in Mathematica by using the function MatrixForm:

m1//MatrixForm

You can check whether m1 is a matrix as follows:

MatrixQ[m1]

Cell 2.  Addition of two matrices:

If (Subscript[a, ij]) is the ij-th element of matrix A and (Subscript[b, ij]) is the ij-th element of matrix B than (Subscript[a, ij])+(Subscript[b, ij]) is the ij-th element of the matrix 
	C ≡ A +B. 
In other words
	(Subscript[C, ij]) = (Subscript[A, ij]) + (Subscript[B, ij])

Create two 3×3 matrices m1 and m2 to work with later:

m1=Table[(Subscript[a, i,j]), {i,1,3},{j,1,3}]

m1//MatrixForm

m2=Table[(Subscript[b, i,j]), {i,1,3},{j,1,3}]

m2//MatrixForm

To add these matrices we use the operator + (see below)

m1+m2//MatrixForm

If you feel like it you can use

Plus[m1,m2]//MatrixForm

Cell 3:  Multiplication with a scalar: 
	means multiply every element of the matrix with the scalar.

Below I multiply the matrix m1 with 3

m1//MatrixForm
3 m1//MatrixForm

Cell 4: The product of two matrices

The definition of 
	
	C= A B 
	
	is that the ij-th element of C is 

	(Subscript[C, i,j])=(∑_(n=1)^N)(Subscript[A, i,n]) (Subscript[B, n,j])

The function Dot performs the multiplication of two matrices. The shorthand is A.B

Multiply the matrices m1 and m2:

Dot[m1,m2]

m1.m2

In what follows we only work with square matrices (N rows and N columns). You can only multiply square matrices having the same dimension. See what happens if you try to multiply a 3×3 matrix with a 4×4 matrix:

A=Table[(Subscript[a, i,j]),{i,1,3},{j,1,3}];
B=Table[(Subscript[b, i,j]),{i,1,4},{j,1,4}];
A.B

There are matrices whose product is commutative but in general A.B ≠B.A. If it happens that 

	A.B =B.A 

we say that the matrices A and B commute.

Generate two 5×5 matrices and show that they are not commutative

SeedRandom[7];
A=Table[RandomReal[{-4,4}],{i,1,5},{j,1,5}];
B=Table[RandomReal[{-4,4}],{i,1,5},{j,1,5}];
A.B==B.A

The form X==Y is short-hand for Equal[X, Y] which returns True if X is equal to Y.

HoldForm[a==b]//FullForm

Cell 5: The product of a matrix with a vector or the action of a matrix on a vector or the transformation of a vector by a matrix (these phrases mean the same thing).

If A is a matrix having the components (Subscript[a, i,j]) and v is a vector having the components (Subscript[v, i]) then 

	w ≡ Av   

is the vector   with the components
	
	(Subscript[w, j]) = (∑_(i=1)^N)(Subscript[a, j,i])(Subscript[v, i])			(1)
	
and 
	y≡vA   
	
is the vector w having the components

	(Subscript[w, i]) = (∑_(i=1)^N)(Subscript[v, j])(Subscript[a, j,i])

In Mathematica these operations are performed by the function Dot. Let's see how this works. Generate first the vector vVector having the components (Subscript[v, i])

Clear[vVector,mMatrix,v,m, vVectorC]

vVector=Table[(Subscript[v, i]), {i,1,3}]

Furthermore, define a vector vVectorC whose components are complex numbers

SeedRandom[3]
vVectorC=Table[RandomComplex[{-3-I 3,3+ 3I}],{i,1,3}]

Define a 3×3 matrix mMatrix

mMatrix=Table[(Subscript[m, i,j]), {i,1,3},{j,1,3}]

Now act with the matrix mMatrix on the vector vVector and display the result by using MatrixForm

Dot[mMatrix,vVector]//MatrixForm

As you can see this syntax give the correct result

Test that Dot works correctly by using the definition  

	(Subscript[w, j]) = (∑_(i=1)^N)(Subscript[a, j,i])(Subscript[v, i])
I calculate below the first element of the resulting vector to compare to the result obtained above.

Sum[mMatrix[[1,i]]*vVectorC[[i]],{i,1,3}]

As a short hand for a matrix acting on a vector you can use:

mMatrix.vVector//MatrixForm

This operation is also defined when a vector acts on a matrix from the left: if 

	w ≡ v A (or v.A) 
then

(Subscript[w, i ])= (∑_(i=1)^N)(Subscript[v, j])(Subscript[a, j,i])

This operation is performed by the function Dot. Now the vector is the first argument and the matrix is the second one:

Dot[vVector,mMatrix]

Or in short form

vVector.mMatrix

One can now define the operation 

	w. A .v  
where w and v are vectors and A is a matrix, but you must be careful when the components of w are complex numbers. If you want to calculate 

	⟨w|A v⟩	

then you should use

	scalar[w, A.v]
 
 where the function scalar is the one I defined earlier in this notebook.   If you use Dot or its short hand you get the wrong result (if the components of the vector w are complex numbers):

Let's see how this work by defining two vectors aVector and xVector and the matrix mMatrix and then calculating 

	aVector.mMatrix.xVector

Clear[ aVector,xVector,mMatrix]
aVector=Table[(Subscript[a, i]),{i,1,3}]
xVector=Table[(Subscript[x, i]),{i,1,3}]
mMatrix=Table[(Subscript[m, i,j]),{i,1,3},{j,1,3}]

aVector.mMatrix.xVector

Note that  is A is a matrix and x is a vector than  A.x is a vector.  If y is a vector then y.A.x is the scalar product of vector A.x with the vector y and therefore y.A.x is a number (i.e. it is not a vector or a matrix)

If y is a vector whose components are complex numbers (or complex expressions (i.e. meaning they contain (Sqrt[-1])) do not use Dot: use the function scalar defined earlier:

scalar[aVector, mMatrix.xVector]

In quantum mechanics a.M.x this is denoted ⟨a|M|x⟩. Therefore we can identify a bra ⟨a| with acting with the vector a on M.x from the left .
	⟨a |x⟩ is equivalent to scalar[a,x]
	⟨a | M x⟩ is equivalent to scalar[a,M.x]

Section 3: Various kinds of matrices

Unitary and Hermitian matrices are routinely used in Quantum Mechanics. In this section I defined these matrices and write programs that will generate numerical examples of various types of matrices.  You will learn how to generate a unit matrix, the transpose of a matrix, the adjoint of a matrix and the inverse of a matrix. You will also learn how to generate numerical examples of Hermitian and unitary matrices.

Cell 1 : The unit matrix

The unit matris I is a matrix which causes no change when acting on a vector or on other matrix. The unit matrix is for matrix algebra what the number 1 is for multiplication arithmetics.

Mathematica has the function IdentityMatrix[n] that generates the n×n unit matrix:

IdentityMatrix[4]//MatrixForm

If we act with this matrix on a vector v the result is v
	I.v = v
Check that this is true:

IdentityMatrix[3].{a1, a2,a3}

Multiplying a matrix A with the identity matrix does not cause any change
	I.A=A

matrix=Table[(Subscript[a, i,j]),{i,1,3},{j,1,3}]
IdentityMatrix[3].matrix

The unit matrix commutes with any matrix.

Cell 2:  The transpose of a matrix

If the elements of a matrix M are (Subscript[m, i,j ]), then the transpose (M^T) of M has the elements 

	(Subscript[((M^T)), i,j]) = (Subscript[m, j,i])

In Mathematica the function Transpose[m] gives the transpose of the matrix m:

Clear[mMatrix]
mMatrix=Table[(Subscript[m, i,j]),{i,1,3},{j,1,3}]
{mMatrix//MatrixForm,Transpose[mMatrix]//MatrixForm}

Cell 3:  The adjoint of a matrix.

If the matrix M has the elements (Subscript[m, i,j])  then its adjoint (denoted (SuperDagger[M])) has the elements
	
	(Subscript[((SuperDagger[M]) ), i,j]) = (Subscript[m, j,i]^( * ))
The symbol * indicates complex conjugation (replace  i=(Sqrt[-1]) with -i).  To obtain the adjoint of M you transpose M and then you take the complex conjugate of all the elements

To experiment, generate a matrix having complex elements:

SeedRandom[4];
matrixM=Table[RandomComplex[{-3-4 I,3+4 I}],{i,1,3},{j,1,3}];
matrixM;
matrixM//MatrixForm

I define the function adjoint which given a matrix  M as its argument it returns the adjoint of M.

Unprotect[adjoint];
Clear[adjointM];
adjoint[M_?MatrixQ] :=Module[{adjointM,n},
n=Length[M];
adjointM=Table[  0,{i,1,n},{j,1,n}]; 
Table[adjointM[[i,j]]=Conjugate[ M[[j,i]]] ,{i,1,n},{j,1,n}];
adjointM];
Protect[adjoint];

Look at a simple example to see if this works properly

Clear[a,b,c,d,ma]
ma={{a,b},{c,d}};
ma//MatrixForm
  adjoint[ma]//MatrixForm

One more check that the function works properly

adjoint[matrixM]//MatrixForm

matrixM//MatrixForm

Cell 4:  The inverse of a matrix

The inverse of a matrix A is denotede (A^(-1)) and it is defined by the equation

	(A^(-1))A =(AA^(-1)) =I
where I is the unit matrix. Not all matrices have an inverse. A necessary and sufficient condition that  a matrix has an inverse is that its determinant differs from zero. The inverse is often used to solve  equations of the form 
	A x = y
where A is a matrix x is a vector with unknown components and y is a vector with known components. 
The solution is 
	x =(A^(-1))y

Mathematica calculates the invers of a matrix m by using the function Inverse[m] and the determinant by using the function Det[m]

Generate a matrix with complex elements:

m=.
SeedRandom[5]
m=Table[RandomComplex[{-3-4 I,3+4 I}],{i,1,3},{j,1,3}];
matrixM

Calculate its determinant

Det[m]

This is not zero therefore the matrix has an inverse:

mInv=Inverse[m];
mInv//MatrixForm

Test whether the inverse satisfies the definition 
	mInv.m = Identity

m.mInv//Chop//MatrixForm

Also

mInv.m//Chop//MatrixForm

Let us use matrix inversion to solve the system of equations

	m.x =y

Generate the components of the vector y.

SeedRandom[4];
y=Table[RandomComplex[{-2-4 I,2+4 I}],{i,1,3}];
y

Find x that solves equation   m.x = y

solution=Inverse[m].y

Use the function Solve to test whether you get the same result

Generate the unknown components of the vector x

x={x1, x2, x3}

The system of equations is:

eqs=Table[(m.x)[[i]]==y[[i]],{i,1,3}]

I can solve these equation by using solve.

Solve[eqs,{x1,x2,x3}]

Compare to the solution provided by using the inverse matrix

solution

The two procedures give the same result.

Cell 5:   Hermitian matrices

Matrices that are equal to their adjoint are called Hermitian matrices. Their elements satisfy

	(Subscript[M, i, j]) = (Subscript[M, j, i]^( * ))

These matrices are essential in quantum mechanics since all observables are represented by Hermitian matrices.  We are going to experiment with Hermitian matrices and therefore I will create the function hermitian which produces a hermitian matrix of desired dimension with real or complex elements.

Unprotect[hermitian];
Clear[hermitian];
hermitian::usage="The function hermitian[dimension][seed, c, {x,y}] generates a dimension×dimension
hermitian matrix. If c=Real the matrix has real elements. If c=Complex the matrix has complex elements. 
The arguments {x,y} give the limits for the Random function which generate the elements. If c=Real then x and y are
real numbers; if c=Complex then x and y are complex numbers. The argument seed sets the seed in SeedRandom. If you use the same seed in two different calculations you get the same matrix. To get 
a new matrix use a new seed.";
hermitian[dimension_Integer][seed_Integer,c_,{a_,b_}]:=Module[{m},
 SeedRandom[seed];
m=Table[0,{i,1,dimension},{j,1,dimension}];
Table[m[[i,j]]=Random[c,{a,b}],{i,1,dimension},{j,i+1,dimension}];
Table[m[[j,i]]=Conjugate[m[[i,j]]] ,{i,1,dimension},{j,i+1,dimension}];
Table[m[[i,i]]:=Random[Real,{Re[a],Re[b]}],{i,1,dimension}] ;
m
];
Protect[hermitian];

The function creates first a matrix m whose elements are all equal to zero. Then it fills the upper triangle in the matrix with random numbers (real or complex depending of the value of the argument c). Then it defines  the lower triangle by using the rule (Subscript[M, i, j]) = (Subscript[M, j, i]^( * )).  Then it fills the diagonal with real numbers (a Hermitian matrix must have real numbers on the diagonal).

Let's generate a 4×4 Hermitian matrix whose elements are real numbers between -5 and 5. I took the seed for the random number generator equal to 7.

MR=hermitian[4][7,Real,{-5,5}] ;
MR//MatrixForm

Create now a complex Hermitian matrix

MC=hermitian[4][7,Complex,{-5-I 5,5+ 3 I}];
MC//MatrixForm

Cell 6:   Unitary matrices

There are several equivalent definitions of an unitary matrix U. One of them is 

	 (SuperDagger[U])U = U(SuperDagger[U])=I
	
where (SuperDagger[U]) is the Hermitian conjugate of U and I is the unit matrix.

A simple way of generating an unitary matrix is to use a theorem that says that the eigenvectors of a hermitian matrix form an unitary matrix.  Write a function unitary that returns an unitary matrix of desired dimension with real or complex elements

Unprotect[unitary];
Clear[unitary];
unitary::usage="The function unitary[dimension][seed, a, {x,y}] generates a dimension×dimension
unitary matrix. If c==Real the matrix has real elements. If c==Complex the matrix has complex elements. 
The arguments {x,y} give the limits for the Random function which generate the elements. If c==Real then x and y are
real numbers; if c=Complex then x and y are complex numbers. seed sets the seed in SeedRandom. To get 
a new matrix use a new seed. The rows of the unitary matrix are the eigenvectors of a hermitian matrix";
unitary[dimension_Integer][seed_Integer,c_,{a_,b_}]:=Module[{m},
 SeedRandom[seed];
m=Table[0,{i,1,dimension},{j,1,dimension}];
Table[m[[i,j]]=Random[c,{a,b}],{i,1,dimension},{j,i+1,dimension}];
Table[m[[j,i]]=Conjugate[m[[i,j]]] ,{i,1,dimension},{j,i+1,dimension}];
Table[m[[i,i]]:=Random[Real,{Re[a],Re[b]}],{i,1,dimension}] ;
Eigenvectors[m]
];
Protect[unitary];

The line before the last create a Hermitian matrix m.  The last line finds the eigenvectors of m. The output is a matrix in which each row is an eigenvector of the Hermitian matrix m.  There is a theorem that shows that such a matrix is unitary.

Generate an unitary matrix whose elements are real numbers

uR=unitary[4][2, Real,{-4,6}];
uR//MatrixForm

Verify some of the known properties of  unitary matrices

1. The absolute value of the determinant of an unitary matrix is equal to 1

2. The elements (Subscript[m, i,j]) of an unitary matrix satisfy the equations
		 (∑_(m=1)^N)(Subscript[m, m,i]^( * ))(Subscript[m, m,j])=(Subscript[δ, ij])					(1)

	and

		(∑_(i=1)^N)(Subscript[m, im])(Subscript[m, jm]^( * ))=(Subscript[δ, ij])					(2)
These are equivalent to  (SuperDagger[U])U = U(SuperDagger[U])=I

3. (SuperDagger[U]) =(U^(-1))

4. The inverse (U^(-1)) of an unitary matrix U is also an unitary matrix

1.  Generate a unitary matrix with complex elements to test that the absolute value of its determinant is equal to 1. :

uC=unitary[3][7,Complex,{-6-I 4,6+ I 4}];
uC//MatrixForm

The absolute value of the deteminant of an unitary matrix  must be equal to 1. Check this:

Det[uC]//Abs

The determinant itself is a complex number (whose absolute value is equal to 1)

Det[uC]

2.  An unitary matrix M satisfy the equations  (SuperDagger[MM])= I and  (SuperDagger[M])M=I 
which are equivalent to the following conditions for the matrix elements  (Subscript[m, i,j])  

(∑_(m=1)^N)(Subscript[m, m,i]^( * ))(Subscript[m, m,j])=(Subscript[δ, ij])					(1)

and

(∑_(i=1)^N)(Subscript[m, im])(Subscript[m, jm]^( * ))=(Subscript[δ, ij])					(2)

Let's test them for the unitary matrices that we have created

Test first 
	(SuperDagger[(OverHat[M])]).(OverHat[M])  = (OverHat[I])
Use the function adjoint defined earlier

adjoint[uC].uC//Chop//MatrixForm

Check that this is equivalent to Eq. 1

g[i_Integer,j_Integer][u_?MatrixQ]:=Sum[Conjugate[u[[m,j]]]*u[[m,i]],{m,1,Length[u]}]

g[1,1][uC]//Chop

g[1,2][uC]//Chop

g[1,3][uR]//Chop

Check that if U is unitary (U^(-1))is also unitary

Inverse[uR]//Det//Abs

Check that the matrix we generate satisfies the theorem that the inverse of an unitary matrix is equal to its adjoint. Check this:

Inverse[uC]//Chop//MatrixForm
adjoint[uC]//Chop//MatrixForm

Section 4:  The eigenvalue problem

The eigenvalue problem for a matrix (or an operator) requires us to find the vector x and the scalar λ that satisfy the equation

	M.x =λ x

where M is a NXN matrix and x is an N-dimensional vector. 
I will only deal here with unitary and hermitian matrices, since these are the most common in quantum mechanics.

Mathematica provides several functions for solving the eigenvalue problem:
	Eigenvalues[m]  gives the eigenvalues of the matrix m
	Eigenvectors[m] gives the eigenvectors of the matrix m
	Eigensystem[m] gives the eigenvalues of m and the corresponding eigenvectors. 
	CharacteristicPolynomial[m, x] gives the characteristic polinomial of the matrix m, having the variable x.
	
These functions have various options that I will explain below.

Cell 1:  Create a hermitian matrix M, find its eigenvalues and test some of their properties

Quantum mechanics deals with Hermitian and Unitary matrices.  I study here the eigenvalues of these two kind of matrices.  I copy here the function hermitian from Cell 4 of Section 3, which creates a Hermitian matrix.

Unprotect[hermitian];
Clear[hermitian];
hermitian::usage="The function hermitian[dimension][seed, c, {x,y}] generates a dimension×dimension
hermitian matrix. If c=Real the matrix has real elements. If c=Complex the matrix has complex elements. 
The arguments {x,y} give the limits for the Random function which generate the elements. If c=Real then x and y are
real numbers; if c=Complex then x and y are complex numbers. The argument seed sets the seed in SeedRandom. To get 
a new matrix use a new seed.";
hermitian[dimension_Integer][seed_Integer,c_,{a_,b_}]:=Module[{m},
 SeedRandom[seed];
m=Table[0,{i,1,dimension},{j,1,dimension}];
Table[m[[i,j]]=Random[c,{a,b}],{i,1,dimension},{j,i+1,dimension}];
Table[m[[j,i]]=Conjugate[m[[i,j]]] ,{i,1,dimension},{j,i+1,dimension}];
Table[m[[i,i]]:=Random[Real,{Re[a],Re[b]}],{i,1,dimension}] ;
m
];
Protect[hermitian];

Exercise:  the matrix   A+(SuperDagger[A]) is hermitian. Prove that this is so and then use this formula to write a program which creates a Hermitian matrix.

Create a 4 × 4 hermitian matrix m with complex elements to experiment with:

SeedRandom[3];
M=hermitian[4][3,Complex,{-5-5 I,5+5 I}]

I use SeedRandom to ensure that I get the same random numbers when I execute this evaluation at some other time.   I wanted to make sure that if you try to repeat my calculation you get the same matrix as I did.

Cell 1:  Calculate the eigenvalues of the hermitian matrix M

The eigenvalues of a matrix M are calculated by the function Eigenvalue[M]

eigVal=Eigenvalues[M]//Chop

In the output of this function the eigenvalues are sorted so that the one with the largest absolute value is given first, the one with the second largest absolute value is second etc.  This makes absolutely no sense to me.  In quantum mechanics we are often interested in the smallest eigenvalue.  You can extract the smallest eigenvalue from the result given above by using the function Min

Min[eigVal]

You can sort the eigenvalues so that the smallest is the first, the second smallest is the second, etc.

Sort[eigVal]

The eigenvalues of a hermitian matrix must be real numbers and our example conforms to this theorem. We also expect that the trace of the matrix equals the sum of the eigenvalues (this is true for any matrix).  Check this

The trace of a matrix is (by definition) the sum of the diagonal elements. For our matrix

trace1= Sum[M[[i,i]],{i,1,4}]

Mathematica provides the function Tr[matrix] which gives the trace of a matrix

Tr[M]

There is a theorem that says that the trace is equal to the sum of the eigenvalues. Calculate that sum

trace2=Plus@@Eigenvalues[M]//Chop

The calculation agrees with the theorem

A theorem specifies that the product of the eigenvalues of a matrix is equal to the determinant of the matrix (this is true for all matrices not just the Hermitian ones). Check this property for M by using the function Det[A] to calculate the determinant of A

Det[M]//Chop

Calculate the product of the eigenvalues and compare its value to that of the determinant

Times@@Eigenvalues[M]//Chop

The calculation conforms with the theorem

If you need to do further work with the eigenvalues you can assign them to symbolic variables. 
In the example below I want to assign the eigenvalues to the variables eva[1], eva[2],  etc.  I will also sort the eigenvalues so that eva[1] is the smallest, eva[2] is the second smallest etc

{eva[1], eva[2], eva[3], eva[4]}=Eigenvalues[M]//Chop//Sort

You can check that the assignment is correct

eva[1]

Or you can make the following assignment:

eigVal=Eigenvalues[M]//Chop//Sort

You can retrieve any one of these eigenvalues by using the following syntax:

eigVal[[2]]

You can also turn these assignments into replacement rules if you wish

ruleEigenValues=MapThread[#1->#2&,{{a[1],a[2],a[3],a[4]},Eigenvalues[M]//Chop//Sort}]

The syntax of the above statement may be obscure to you. Use Mathematica help (or my lecture notes on Mathematica) to learn how MapThread works. Also you need to know how pure functions work. I give you a little help below:

Here is what MapThread does:

res=MapThread[f,{{a,b,c},{A,B,C}}]

Now replace f with the nameless function (#1 -> #2) &

res/.f->(#1->#2&)

You can find the eigenvalues of a matrix by finding the roots of the characteristic polynomial of the matrix. 
 The characteristic polynomial of the matrix M is defined to be the determinant of M-z I where I is the unit matrix.
 The function CharacteristicPolynomial[M,z] gives the characteristic polynomial of the matrix M and the variable in the polynomial is z.

cp=CharacteristicPolynomial[M,z]//Chop

You can verify that the characteristic polynomial given by this function is the same as the one  given by the definition

	determinant[m-z I]

where I is the unit matrix and z is the unknown variable in the polynomial. z can take real or complex numerical values (it is not a matrix or a vector)

Clear[z]
acp=Det[M-z IdentityMatrix[4]]//Chop

Check that the two expressions are the same

cp==acp

Solve the equation 
	characteristic polynomial=0
to test that the solutions are the eigenvalues:

Solve[cp==0,z]//Flatten//Sort

The eigenvalues calculated earlier are:

Eigenvalues[M]//Chop//Sort

The eigenvalues are indeed the roots of the characteristic polynomial.

Cell 2:  If you want only the eigenvectors of the matrix M use   Eigenvectors[M] ( I do not recommend this).

We define below the hermitian matrices matrix MC and MR. MC has complex elements and MR has real ones. I copy the function hermitian, defined earlier, which creates a Hermitian matrix.

Unprotect[hermitian];
Clear[hermitian];
hermitian::usage="The function hermitian[dimension][seed, c, {x,y}] generates a dimension×dimension
hermitian matrix. If c=Real the matrix has real elements. If c=Complex the matrix has complex elements. 
The arguments {x,y} give the limits for the Random function which generate the elements. If c=Real then x and y are
real numbers; if c=Complex then x and y are complex numbers. The argument seed sets the seed in SeedRandom. To get 
a new matrix use a new seed.";
hermitian[dimension_Integer][seed_Integer,c_,{a_,b_}]:=Module[{m},
 SeedRandom[seed];
m=Table[0,{i,1,dimension},{j,1,dimension}];
Table[m[[i,j]]=Random[c,{a,b}],{i,1,dimension},{j,i+1,dimension}];
Table[m[[j,i]]=Conjugate[m[[i,j]]] ,{i,1,dimension},{j,i+1,dimension}];
Table[m[[i,i]]:=Random[Real,{Re[a],Re[b]}],{i,1,dimension}] ;
m
];
Protect[hermitian];

Create a 4 × 4 hermitian matrix MC with complex elements to experiment with:

SeedRandom[3];
MC=hermitian[4][3,Complex,{-5-5 I,5+5 I}]

eigenVec=Eigenvectors[MC]//Chop

The output is a list of four lists, each having four elements.  The first list, namely {-0.313641 + 0.284894 i, -0.215179 + 0.298219 i, -0.0307692 + 0.423997 i, 0.710288}  is an eigenvector.  There are four eigeinvectors because the matrix is 4×4.    I don't like this function because there is  the possibility of getting confused regarding which eigenvector corresponds to which eigenvalue.

You can get the first eigenvector only by using

Eigenvectors[M,1]

Cell 3: If you want the eigenvalues and the eigenvectors of the matrix m use
		Eigensystem[m]

Here is a good way to assign the results

{eVa,eVe}=Eigensystem[MC]//Chop

The first list in the output, namely {9.80944, 8.83082, -7.46447, -0.628336}, is the list of the eigenvalues.  Note that Mathematica orders them according to the absolute values: the eigenvalue that has the largest absoulute value (namely 12.37) is the first, etc. The four lists that follow are the corresponding eigenvectors. Because it is tedious and prone to error to try to keep track of which list is what,  I instructed the computer to store the eigenvalues in the variable eVa and the eigenvectors in eVe.
This makes it easy to use the output. For example the first eigenvalue is

eVa[[1]]

The second one is

eVa[[2]]

The first eigenvector (corresponding to the first eigenvalue) is

eVe[[1]]

I will test that the results satisfy the eigenvalue equation
	MC. eigenvector(i) = eigenvalue(i). eigenvector(i)   
Here i labels the four eigenvectors and the four eigenvalues

(MC.eVe[[1]]//Chop)==(eVa[[1]]*eVe[[1]]//Chop)

I need the Chop function to eliminate roundoff errors that will trick the computer into thinking that the two numerical results are different. See what happens if I don’t chop

MC.eVe[[1]] == eVa[[1]]*eVe[[1]]

The equation m v = e v is satisfied for the first eigenvector and eigenvalue. You can check that it is satisfied for all other eigenvectors and eigenvalues.

Cell 4: The eigenvectors of a hermitian matrix that belong to different eigenvalues are orthogonal to each other and are normalized.

The scalar product of two N-dimensional vectors is defined as:

	⟨v|u⟩ = (∑_(i=1)^N)(Subscript[v, i]^( * ))(Subscript[u, i])

To evaluate it I created earlier the function scalar[v,u], which I copy below

Unprotect[scalar];
Clear[scalar]
scalar[x_?VectorQ,y_?VectorQ]:= 
Conjugate[x].y
Protect[scalar];

Denote the eigenvector of the matrix MC by v[1], v[2], v[3], v[4]. This assignment is made in the next line of code

{eigenvalues, {v[1],v[2],v[3],v[4]}}=Eigensystem[MC]//Chop

Check that the assignment worked

v[2]

This is indeed the second eigenvector.  Now calculate ⟨v[i]|v[j]⟩. We find that if i ≠j the scalar product is zero (meaning that v[i] and v[j] are orthogonal) and if i=j the scalar product is equal to 1 (meaning that each eigenvector is normalized). In previous version of Mathematica the eigenvectors were not normalized, so always check normalization.

Table[scalar[v[i],v[j]],{i,1,4},{j,1,4}]//Chop//MatrixForm

Cell 5: The eigenvectors of a hermitian matrix form an unitary matrix X

The eigenvalue problem for an N-dimensional matrix M is 

	M.x = λ x		(1)
There are m vectors that satisfy this equation. I denote the α-th eigenvectors by 

	 x(α) ≡{(Subscript[x, 1])(α), (Subscript[x, 2])(α), ..., (Subscript[x, N])(α) }			(2)
	
(Subscript[x, j])(α) is the j-th component of α-th eigenvector. A theorem says that the matrix X  

X =   (({{({{(Subscript[x, 1])(1), (Subscript[x, 2])(1)}, {(Subscript[x, 1])(2), (Subscript[x, 2])(2)}}), ({{(Subscript[x, 3])(1), (Subscript[x, 4])(1)}, {(Subscript[x, 3])(2), (Subscript[x, 4])(2)}})}, {({{(Subscript[x, 1])(3), (Subscript[x, 2])(3)}, {(Subscript[x, 1])(4), (Subscript[x, 2])(4)}}), ({{(Subscript[x, 3])(3), (Subscript[x, 4])(3)}, {(Subscript[x, 3])(4), (Subscript[x, 4])(4)}})}}))				(3)

is unitary. Note that each row of  matrix X is a different eigenvector of the Hermitian matrix M.

SeedRandom[3];
M=hermitian[4][3,Complex,{-5-5 I,5+5 I}]

Create X which is a matrix in which each row is the eigenvector of the Hermitian matrix M.
The matrix X is unitary (so the theory says). We can check this. Create X:

X=Eigenvectors[M]//Chop

One of the definitions of an unitary matrix is that the absolute value of any of its eigenvalues is equal to 1 . Let's check if this is true for the matrix X defined above

Eigenvalues[X]//Abs

Also the absolute value of the determinant of X should be equal to 1.

Det[X]//Abs

The expected relationship checks out.

Cell 6:  If M is Hermitian there is an unitary matrix S such that 

	(S^(-1))MS
	
is a diagonal matrix, whose diagonal elements are the eigenvalues of M.  We say that the matrix S diagonalizes matrix M. The matrix (S^(-1))MS has the same eigenvalues as M but its eigenvectors are different from those of M. Since only the eigenvalues are measured experimentally, the matrix M is physically equivalent to S.  Since the eigenvectors are not measurable quantities it does not matter that those of M are different from those of  (S^(-1))MS. However, sometimes the eigenvectors of (S^(-1))MS may be easier to interpret physically than those of M.  One example is the (sp^3) hybridization in which the hybrid orbitals are obtained from the s and p orbitals of carbon by an unitary transformation.

There is a theorem that states that the N×N matrix S is constructed by using the eigenvectors of M as columns in S (first eigenvector of M is the first column in S). The easiest way to obtain S in Mathematica is by  using the syntax Transpose[Eigenvectors[M]].  This works because Eigenvectors[M] is a matrix whose rows are the eigenvectors of M; Transpose[Eigenvectors[M]] takes these rows and constructs a new matrix in which the former rows are now columns.  
	The transformation that takes matrix M into the matrix (A^(-1))M A, where A is a matrix that has an inverse, is called a similarity transformation. Some sources use this term even when A is not unitary others confine the term to unitary matrices A .

Let us test the statements made above about the equivalence of M and  (S^(-1))M S  
Create a Hermitian matrix MC

Unprotect[hermitian];
Clear[hermitian];
hermitian::usage="The function hermitian[dimension][seed, c, {x,y}] generates a dimension×dimension
hermitian matrix. If c=Real the matrix has real elements. If c=Complex the matrix has complex elements. 
The arguments {x,y} give the limits for the Random function which generate the elements. If c=Real then x and y are
real numbers; if c=Complex then x and y are complex numbers. The argument seed sets the seed in SeedRandom. If you use the same seed in two different calculations you get the same matrix. To get 
a new matrix use a new seed.";
hermitian[dimension_Integer][seed_Integer,c_,{a_,b_}]:=Module[{m},
 SeedRandom[seed];
m=Table[0,{i,1,dimension},{j,1,dimension}];
Table[m[[i,j]]=Random[c,{a,b}],{i,1,dimension},{j,i+1,dimension}];
Table[m[[j,i]]=Conjugate[m[[i,j]]] ,{i,1,dimension},{j,i+1,dimension}];
Table[m[[i,i]]:=Random[Real,{Re[a],Re[b]}],{i,1,dimension}] ;
m
];
Protect[hermitian];

SeedRandom[3];
MC=hermitian[4][3,Complex,{-5-5 I,5+5 I}]

Create S so that its columns are the eigenvectors of MC

S=Transpose[Eigenvectors[MC]]

Test that (S^(-1)) MC S is a diagonal matrix:

Inverse[S].MC.S//Chop//MatrixForm

Form a diagonal matrix with the eigenvalues of MC on the diagonal, to compare with the above.

Eigenvalues[MC]//Chop//DiagonalMatrix//MatrixForm

The matrix S transforms the basis set vectors into the eigenvectors of MC

Orthonormal basis sets are essential for all computations in Quantum Mechanics. Let the orthonormal basis set of an N-dimensional ket-space consist of the kets |(Subscript[e, 1])⟩ , |(Subscript[e, 2])⟩, ... |(Subscript[e, N])⟩.   This means that we can write any ket as

	|ψ⟩ = (Subscript[a, 1])|(Subscript[e, 1])⟩ +.....(Subscript[a, N])|(Subscript[e, N])⟩			(1)
	
where  (Subscript[a, n]), n=1, 2, ... N are complex numbers. 
	In all calculations in Quantum Mechanics the basis set is known. This means that we can evaluate all operations in which the kets |(Subscript[e, n])⟩ are involved.  The equation 1 puts the ket |ψ⟩ in correspondence with the N-dimensional vector 

		{(Subscript[a, 1]), (Subscript[a, 2]),....,(Subscript[a, N])}			(2)
If we know the numbers (Subscript[a, n]), n=1, 2, ...N we know  |ψ⟩. There is a one-to-one correspondence between any ket and a specific N-dimensional vector (such as the one in Eq. 2) 
	In this section I will:
1. Derive the vectors that correspond to the kets forming the basis set |(Subscript[e, i])⟩
2. Show that the matrix S that diagonalizes a Hermitian matrix MC gives the eigenvectors of MC when it acts on the basis set.  In other words S|(Subscript[e, i])⟩ is an eigenvector of MC.

What N-dimensional vector represents the basis set ket |(Subscript[e, i])⟩?  The general formula for finding the component (Subscript[a, i]) of the vector corresponding to a ket |ψ⟩ is

	(Subscript[a, i]) = ⟨ (Subscript[e, i])|ψ ⟩					(3)
	
	Using this definition I can calculate what is the vector representing the ket  |(Subscript[e, 1])⟩.  The j-th component of this vector is (because the basis set kets are orthonormal)

	(Subscript[a, j]) = ⟨ (Subscript[e, j])|(Subscript[e, 1])⟩ = (Subscript[δ, j,1])

where (Subscript[δ, j,1]) is the Kronecker delta.  The vector representing |(Subscript[e, 1])⟩ is  therefore 

	e(1) ={1, 0, ...,0}
	
Similarly, the vector representing |(Subscript[e, 2])⟩  is 

	e(2)= {0, 1, 0, ...,0}
	
The vectors e(1), e(2), ... e(N) represent the basis set kets |(Subscript[e, 1])⟩ , ... |(Subscript[e, n])⟩ in the N-dimensional vector space.

In what follows we will do some calculations using these basis set vectors. First I need to define them and I use the function IdentityMatrix. I will use a 4-dimensional vector space.

Unprotect[e];
Clear[e]
{e[1],e[2],e[3],e[4]}=IdentityMatrix[4] 
Protect[e];

Show that  a[1] e[1] + a[2] (Subscript[e, 2])+...+a[N] e[N]  is the vector having the components a[1], ...,a[N].

Sum[a[i]*e[i],{i,1,4}]

Now let us calculate how S acts on the vectors representing the basis set . Is S e[1] and eigenvector of MC?  To see if this is true I calculate S e[1] and subtract from it the first eigenvector of MC. If they are the same I should get a vector with all components equal to zero.

( S.e[1] - (Eigenvectors[MC,1]//Flatten ))//Chop

You can do a similar calculation to see that S e[i] gives the i-th eigenvector of MC.

You can also get the basis set vectors by applying (S^(-1)) to the eigenvectors of M:

Inverse[S].(Eigenvectors[MC,1]//Flatten)//Chop

Inverse[S].(Eigenvectors[MC,{2}]//Flatten)//Chop

If you write out formally what these operations do, you will find that this is not miraculous at all. 
Summary:  the matrix S  constructed  by placing the eigenvectors of MC as columns in S has interesting properties:
	S is unitary
	(S^(-1)).M.S is a diagonal matrix with the eigenvalues of M on diagonal
	S.{1,0,0,...} gives the first eigenvector of M
	S.{0,1,0...} gives the second eigenvector of M etc. 
	Conversely:  if v(i) is an eigenvector of MC then (S^(-1))v(i) = e(i)  is a basis set vector.

Cell 7: The eigenstates and the eigenvectors of an unitary matrix

The eigenvalues (Subscript[λ, n]) of an unitary matrix satisfy
	|(Subscript[λ, n])| = 1.
and the eigenvectors calculated with Eigensystem are orthonormal

Create a 4×4 unitary matrix by using the function unitary

Unprotect[unitary];
Clear[unitary];
unitary::usage="The function unitary[dimension][seed, a, {x,y}] generates a dimension×dimension
unitary matrix. If c==Real the matrix has real elements. If c==Complex the matrix has complex elements. 
The arguments {x,y} give the limits for the Random function which generate the elements. If c==Real then x and y are
real numbers; if c=Complex then x and y are complex numbers. seed sets the seed in SeedRandom. To get 
a new matrix use a new seed. The rows of the unitary matrix are the eigenvectors of a hermitian matrix";
unitary[dimension_Integer][seed_Integer,c_,{a_,b_}]:=Module[{m},
 SeedRandom[seed];
m=Table[0,{i,1,dimension},{j,1,dimension}];
Table[m[[i,j]]=Random[c,{a,b}],{i,1,dimension},{j,i+1,dimension}];
Table[m[[j,i]]=Conjugate[m[[i,j]]] ,{i,1,dimension},{j,i+1,dimension}];
Table[m[[i,i]]:=Random[Real,{Re[a],Re[b]}],{i,1,dimension}] ;
Eigenvectors[m]
];
Protect[unitary];

Find the eigenvalue (store them in eval) and eigenvectors (store them in evec) of this matrix by using the function Eigensystem.

Clear[eval,evec]
{eval,evec}=Eigensystem[unitary[4][5,Complex,{-5-I 5,5+I 5}]]

Test whether the absolute values of the eigenvalues are equal to 1.

Table[Abs[eval[[i]]],{i,1,4}]

Test whether the eigenvectors are orthonormal

Table[Conjugate[evec[[i]]].evec[[j]],{i,1,4},{j,1,4}]//Chop//MatrixForm

Cell 8: If A is an arbitrary matrix that has an inverse then the matrix

	S ≡ (A^(-1))M A
	
has the same eigenvalues as M. Moreover if 

	S.x(i) =λ(i) x(i)
then  
	M .A .x(i) = (Subscript[λ, i]) A .x(i).  
	
If x(i) and λ(i) are eigenvectors and eigenvalues of S 
	then
A.x(i) is an eigenvector of M and λ(i) is an eigenvalue of M.  This is very trivial to prove and we verify it only to give you some practice in matrix manipulations in Mathematica

Create an arbitrary 4×4 matrix A with complex elements:

SeedRandom[3]
A=RandomComplex[{-3-3I,3+3 I},{4,4}]

Create an arbitrary matrix M (note that A does not have to be unitary and M does not have to be Hermitian)

SeedRandom[5]
M=RandomComplex[{-3-3I,3+3 I},{4,4}]

Calculate S =(A^(-1))M A

S=Inverse[A].M.A

Calculate the eigenvalues of S

{λ[1],λ[2],λ[3],λ[4]}=Eigenvalues[S]//Chop

Calculate the eigenvalues of M:

Eigenvalues[M]//Chop

The two sets of eigenvalues are obviously the same

Now calculate:
1. The eigenvectors of S and call them x[1], x[2], x[3] and x[6]

Clear[x]
{x[1],x[2],x[3],x[4]}=Eigenvectors[S]

2. Verify that 
	M A x[i] = (Subscript[λ, i]) A x[i]

M.A.x[1]-λ[1] A.x[1]//Chop

You can repeat this for the other eigenvectors and eigenvalues

Section 5:  The bra and ket representation of operators

A Hermitian or an unitary operator (OverHat[A]) in Quantum Mechanics is routinely represented by 

	(OverHat[A]) = (∑_(i=1)^∞)|(Subscript[a, i])⟩ (Subscript[a, i]),⟨(Subscript[a, i])|					(1)
Here (Subscript[a, i]) are the eigenvalues of (OverHat[A]), and |(Subscript[a, i])⟩ are the eigenkets, which means that these quantities satisfy the equation  

	(OverHat[A]) |(Subscript[a, i])⟩ = (Subscript[a, i]) |(Subscript[a, i])⟩					(2)
	
The quantity ⟨(Subscript[a, i])| is called the bra corresponding to the ket |(Subscript[a, i])⟩.  An arbitrary bra ⟨y| becomes meaningful only when it acts on a ket |x⟩ to give ⟨y|x⟩, which is a complex number that we can calculate. 

In this section we verify whether Eq. 1 does represent correctly an operator. We will use a 4-dimensional space for testing this. In this space the operator is represented by a 4×4 matrix. 

 Eq. 1, which contains a bra, is meaningful only when it acts on some ket |x⟩. Therefore we should verify whether 

	(OverHat[A])|x⟩ = (∑_(i=1)^N)|(Subscript[a, i])⟩ (Subscript[a, i]), ⟨(Subscript[a, i])|x⟩					(3)
	
In our 4-dimensional representation A is a 4×4 matrix, |(Subscript[a, i])⟩ is its i-th eigenvector and (Subscript[a, i]) is its i-th eigenvalue. This equation should be true for any n-dimensional vector and any n×n Hermitian matrix. 

Below I will:
1. Create an arbitrary 4-dimensional vector x
2. Create an arbitrary 4×4 Hermitian matrix A
3. Calculate the eigenvalues and the eigenvectors of A
4. Test whether Eq. 3 is correct.

Create the vector x

x=Table[Random[Complex,{-5.-5.I, 5.+ 5.I}],{i,1,4}]

Create a Hermitian matrix A. Use the function hermitian which I defined earlier and reproduce here

Unprotect[hermitian];
Clear[hermitian];
hermitian::usage="The function hermitian[dimension][seed, c, {x,y}] generates a dimension×dimension
hermitian matrix. If c=Real the matrix has real elements. If c=Complex the matrix has complex elements. 
The arguments {x,y} give the limits for the Random function which generate the elements. If c=Real then x and y are
real numbers; if c=Complex then x and y are complex numbers. The argument seed sets the seed in SeedRandom. If you use the same seed in two different calculations you get the same matrix. To get 
a new matrix use a new seed.";
hermitian[dimension_Integer][seed_Integer,c_,{a_,b_}]:=Module[{m},
 SeedRandom[seed];
m=Table[0,{i,1,dimension},{j,1,dimension}];
Table[m[[i,j]]=Random[c,{a,b}],{i,1,dimension},{j,i+1,dimension}];
Table[m[[j,i]]=Conjugate[m[[i,j]]] ,{i,1,dimension},{j,i+1,dimension}];
Table[m[[i,i]]:=Random[Real,{Re[a],Re[b]}],{i,1,dimension}] ;
m
];
Protect[hermitian];

Create the matrix A

A=hermitian[4][7,Complex,{-5-5 I,5+ 5 I}]

Find the eigenvalues and the eigenvectors of A

{eigenva,eigenve}=Eigensystem[A]//Chop

Test that the assignments work

eigenva[[2]]

eigenve[[1]]

I have to calculate scalar products so I need the function scalar which I defined earlier. I redefine it here

Unprotect[scalar];
Clear[scalar]
scalar[x_?VectorQ,y_?VectorQ]:= 
Conjugate[x].y
Protect[scalar];

Test that the eigenvectors of A are orthonormal (they must always be orthogonal but some computer programs do not normalize them).

Table[scalar[eigenve[[i]],eigenve[[j]]],{i,1,4},{j,1,4}]//Chop

They are orthonormal

Calculate A.x

A.x

Compare A.x with (∑_(i=1)^4)|(Subscript[a, i])⟩ (Subscript[a, i]), ⟨(Subscript[a, i])|x⟩ which is calculated below

Sum[eigenve[[i]]*eigenva[[i]]*scalar[eigenve[[i]],x],{i,1,4}]

Obviously the two expressions are equal, which is what we set out to prove.

Section 6:  Functions of matrices

We study some functions, provided by Mathematica, whose argument is a matrix and which return a matrix

Cell 1:  MatrixPower[m,n] where m is a matrix and n is an integer (negative or positive)

MatrixPower[m,n]  gives the n(^th) power of the matrix m.  This is very simple to use

m={{2., 3.4}, {2.1, 5.}}

MatrixPower[m, 2]

Test that the result is correct by multiplying the matrix with itself

m.m

MatrixPower[m,-n] calculates ((((m^(-1))))^n)

MatrixPower[m,-2]

Test

MatrixPower[Inverse[m],2]

The function is capable, in some cases, to calculate symbolic expressions of  (matrix^n)   (assumes that n is integer). 
Calculate the n-th power of the Pauli matrix (Subscript[σ, x])

σx={{0,1},{1,0}};

 MatrixPower[σx,n]

{f[n]/.n->10, MatrixPower[σx,10]}

Here is  an interesting application from Wolfram.

Define the matrix rotInf that rotates a vector by the “infinitesimal” angle θ/n (you’ll take n to infinity at the end of calculation.

Apply the “infinitesimal” rotation n times and call the result nInfRot

Take the limit n to infinity.  The net result is a continuous rotation with the angle θ

rotInf={{1,θ/n},{-θ/n,1}}

nInfRot=MatrixPower[rotInf,n]

nInfRot//MatrixForm

Limit[nInfRot[[1,1]],n->Infinity]//FullSimplify

Limit[nInfRot,n->Infinity]//FullSimplify

This is a rotation by the angle θ around the OZ axis.

Cell 2: MatrixExp[m] and MatrixExp[m,v] give  (e^m) and (e^m). v, respectively, where m is a matrix and v is a vector

MatrixExp[{{1.,2.3},{5.2,8.4}}]

MatrixExp[{{1.,2.3},{5.2,8.4}}].{1.2,3.1}

You can construct any function of a matrix by finding the eigenvalue and the eigenvectors and then using exp[A]= (∑_(i=1)^n)|(Subscript[a, n])⟩ exp[(Subscript[a, n])]⟨(Subscript[a, n])| where |(Subscript[a, n])⟩ are the eigenvectors of A and (Subscript[a, n ])the eigenvalues. Obviously it is handy to have a function like MatrixExp to avoid all these manipulations.

Let’s do some calculations with Pauli matrices

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

It is claimed that 
	exp[ i π (Subscript[σ, x])/4] (Subscript[σ, z])exp[-i π Subscript[σ, x]/4] =(Subscript[σ, y])

ma= FullSimplify/@ (MatrixExp[ I Pi σx/4].σz.MatrixExp[-I Pi σx/4])

The  equation is correct.  It has also been claimed that 
	(Cos[Pi/2] + I (Subscript[σ, x]) Sin[Pi/2] ).(Subscript[σ, z]) = σy
Let’s verify this

(Cos[Pi/2] + I σx Sin[Pi/2] ).σz==σy

The function MatrixExp is very useful for solving equations of the  form

	((dv(t))/dt)= A.v(t)
where v is a n-dimensional vector and A is a time independent, n×n matrix. The formal solution is 
	v(t) = exp[A t] v[0]
Let’s see if this is true. 

Wolfram gives an example

({{x'[t] == z[t]}, {y'[t] == y[t] - z[t]}, {z'[t] == -x[t]}})

The matrix of this system is

a = (({{0, 0, 1}, {0, 1, -1}, {-1, 0, 0}}));

Test that the matrix is correct (it gives the right hand side of the system of differential equations)

Clear[x1,y1,z1]
a.{x1,y1,z1}

It does.  The propagator is

MatrixExp[a t]

Act with this on a particular initial condition
	v[0]={1, 0, 0}

v[0]={1, 0, 0};
MatrixExp[a t ,v[0]]

Use DSolve to test that these solutions are correct (or course they are)

DSolve[{x'[t] == z[t],y'[t] == y[t] - z[t],z'[t] == -x[t],x[0]==1,y[0]==0,z[0]==0}, {x[t],y[t],z[t]},t]

One can find the solution when the matrix is time dependent by applying successively the propagator for short time.  This can also be used in Feit and Fleck method.

MatrixExp[-m] is the inverse of the matrix MatrixExp[m]. Verify this

m={{1.2,8.3},{6.7,9.5}}

Inverse[MatrixExp[m]]//MatrixForm
MatrixExp[-m]//MatrixForm

Cell 3: MatrixLog calculates the logarithm

MatrixLog[m]

The usual relation between the logarithm and the exponential is valid.

m
MatrixLog[MatrixExp[m]]

Cell 4:  MatrixFunction[f,m] uses the power series of the function f to calculate f[m] (where m is a matrix)

For example here is the cosine of the matrix m defined earlier

m

MatrixFunction[Cos,m]

We can play with this.  For example, is  ((cos[m])^2) +((sin[m])^2)= unit matrix ?

MatrixPower[MatrixFunction[Cos,m],2]+MatrixPower[MatrixFunction[Sin,m],2]//Chop

Yes it is! Isn’t this miraculous? No, it isn’t.  This functions are defined by their power series and satisfy any relationship that the power series satisfy,
Let’s try 
	Exp[ I m]= Cos[m] + I Sin[m]
Is this correct for a matrix m?

MatrixExp[ I m]==MatrixFunction[Cos,m] + I MatrixFunction[Sin,m]

MatrixFunction is very versatile: you can use a pure function as an argument

MatrixFunction[ 2. # + 4 #^2&,m]

Verify

2 m + 4 MatrixPower[m,2]

You can use Sin, Cos, Log etc with this syntax

Test whether the relationship
	((d exp[a t])/dt)= a exp[a t]
holds if a is a matrix (it does)

m={{2., 3.4}, {2.1, 5.}}

D[MatrixExp[ m t],t]

m.MatrixExp[m t] //Chop//FullSimplify