Section 1:  The eigenfunctions of commuting matrices.

I work with two dimensional matrices because it saves me typing. There is nothing specific to two-dimensional matrices in what follow.  

I define two matrices A and B, that have the following properties:
	A and B commute
	A has degenerate eigenvalues. 
I picked A to be the unit matrix I because it commutes with every matrix and its eigenvalues are 1 and 1 and are therefore degenerate.  Food for thought: how can you create two matrices that have these properties without using the unit matrix as one of them?

Define matrices A and B.

Remove["Global`*"]

A={{1,0},{0,1}}

B={{1,I},{-I,1}}

{A//MatrixForm, B//MatrixForm}

Show that A and B commute

A.B===B.A

Find the eigenvalues and eigenvectors of A

{{a[1],a[2]},{av[1],av[2]}}=Eigensystem[A]

The eigenvectors of A are orthonormal

Table[av[i].av[j],{i,1,2},{j,1,2}]//MatrixForm

Find the eigenvalues and the eigenvectors of B

{{b[1],b[2]},{bv[1],bv[2]}}=Eigensystem[B]

The eigenvalues of B are not degenerate. The theorem given in the text of the lecture says that  the eigenvectors of B must also be eigenvectors of A, corresponding to the eigenvalues a[1] and a[2]. Test that this is true.

Table[A.bv[i]==a[i] bv[i],{i,1,2}]

It is an accident that the first eigenvector bv[1] corresponds to the eigenvalue a[1]. It could have happened that A.bv[1] ==a[2] bv[1]

The eigenvalues of A are degenerate. The theorem says that the  eigenvectors of A do not have to be eigenvectors of B. Test whether they are eigenvectors of B or not:

{B.av[1], b[1] av[1], b[2] av[1]}

{B.av[2], b[1] av[2], b[2] av[2]}

We see that B.av[1] is not equal to either b[1] av[1] or b[2] av[2].  Therefore av[1] is not an eigenvector of B.  The theorem does not say that it is impossible for a degenerate eigenvector of A to be an eigenvector of B. It only says that unless you are very lucky, the degenerate eigenvectors of A are not eigenvectors of B.

The theorem also tells us that we can construct eigenvectors of B that are also degenerate eigenvectors of A. We can perform this construction because we have some flexibility. Any linear combination of the form

	q =  q1 av[1] + q2 av[2]           (q1 and q2 are arbitrary numbers)
	
is an eigenstate of A corresponding to the degenerate eigenvalue a[1] = a[2]=1. The question is: can we pick q1 and q2 so that q is also an eigenvalue of B? The answer is yes. In the text of the lectures I showed that to do this we need to construct the matrix Q whose elements are:

(Subscript[Q, ij]) = ⟨(Subscript[av, i])|(OverHat[B])|(Subscript[av, j])⟩ = av[i].B.av[j]
	
This is the matrix of the operator B in the basis set {av[1], av[2]}

Construct Q

Q=Table[av[i].B.av[j],{i,1,2},{j,1,2}]

Find the eigenvalues and the eigenvectors of Q.

{{q[1],q[2]},{qv[1],qv[2]}}=Eigensystem[Q]

We can now go back in the original space. The eigenvectors of B in the basis av[i] (eigenvectors of A) are given by

Table[vecB[j]=Sum[qv[j][[i]]*av[i],{i,1,2}],{j,1,2}]

We can test whether vecB[1] and vecB[2] are eigenvectors of B:

Table[B.vecB[j]==b[j]*vecB[j],{j,1,2}]

Check that the vectors  vecB[i] are also eigenvectors of A corresponding to the degenerate eigenvalue (they are by construction)

Table[A.vecB[i]==a[i] vecB[i],{i,1,2}]

Test whether these new eigenvectors are orthonormal. To perform this test I need to define the scalar product for a space that includes vectors whose components are complex numbers.

Clear[scalar]
scalar[a_?VectorQ,b_?VectorQ]:=(a/. Complex[x_,y_]:> Complex[x,-y]).b

Table[scalar[vecB[j],vecB[i]],{i,1,2},{j,1,2}]//MatrixForm

They are orthogonal but they need to be normalized:

Table[newVecB[i]= vecB[i]/Sqrt[scalar[vecB[i],vecB[i] ]],{i,1,2}]

These calculations verified the theorem. They show that because A and B commute their eigenfunctions have the following properties:
	1. The eigenfunctions of B are not degenerate and therefore they are eigenfunctions of A.
	2. The eigenfunctions of A are degenerate. Therefore they do not have to be eigenfunctions of B. However, we can use linear combinations of the degenerate eigenvectors of A to construct two vectors that are eigenvectors of both A and B. 3

Section 2:  Standard deviation in general

Cell 1: Die throw

Consider a die throw.
The probability that the number i will appear is

	P(i)= (1/6)

Write a function that gives this probability

p[i_Integer]:=1/6

The a average value of the results obtained for a large numbers of die throws is 

	⟨ i ⟩ = (∑_(i=1)^6)i P(i)

Calculate the average

av=Sum[i p[i],{i,1,6}]

The standard deviation is given by:

	σ = (Sqrt[⟨ (((i-⟨ i ⟩))^2)⟩])=  (Sqrt[(∑_(i=1)^6)P[i] (((i -⟨i⟩))^2)])
	
Calculate it:

σ1=Sqrt[Sum[p[i]( i-av)^2,{i,1,6}]]//N

Look at the standard deviation for a rigged die.

If we have a rigged die so that 1 and 6 don't ever show up then the probability is

p(1)=0
p(6)=0
p(2)=p(3)=p(4)=p(5)=1/4

Define the function p[i] that gives this probability distribution

Table[p[i]=0,{i,1,6}];
Table[p[i]=1/4,{i,2,5}];
Table[p[i],{i,1,6} (*display the p[i] values*)

Calculate the average:

av=Sum[i p[i],{i,1,6}]

The average for the rigged die is the same as that for the honest one only because of the way we chose the values that have zero probability. The standard deviation is

σ2=Sqrt[Sum[(i-av)^2 p[i],{i,1,6}]]//N

Compare the two standard deviations:

{σ1,σ2}//N

Obviously the standard deviation of the rigged die is smaller and as is our uncertainty regarding the outcome of a throw.

Let us look now at a die that is so rigged that it will always give 3.

Define the probability distribution:

Table[p[i]=0,{i,1,6}];
 p[3]=1 ;
Table[p[i],{i,1,6}]

Calculate the average:

av=Sum[i p[i],{i,1,6}]

The standard deviation is:

σ3=Sqrt[Sum[(i-av)^2 p[i],{i,1,6}]]

Compare the three results:

{σ1,σ2,σ3}//N

Section 3: The standard deviation of coordinate and momentum for a Gaussian wave packet

The system is in the state 
	|η⟩ =  C(∫_(-∞)^∞)exp[-((k^2)/(a^2))] |k⟩ dk  				(1)
The probability that the particle has a wave vector between k and k + dk is 

	(Subscript[P, η])(k) dk =  |⟨k|η⟩((|)^2) dk				(2)
	
Determine C by normalizing the probability. 

You know that when the state is |k⟩ the momentum is ℏk.  Therefore the average momentum is

	(Subscript[⟨p⟩, η]) = ℏ (∫_(-∞)^∞)k (Subscript[P, η])(k)dk				(3)

and 

	(Subscript[⟨(p^2)⟩, η]) = (ℏ^2) (∫_(-∞)^∞)(k^2) (Subscript[P, η])(k)dk				(4)
	
Therefore the standard deviation is:

	(Subscript[σ, η])(p) = (Sqrt[(Subscript[⟨(p^2)⟩, η])-((Subscript[⟨p⟩, η])^2)])				(5)
	
The probability that the coordinate takes values between x and x + dx is 

	(Subscript[P, η])(x)dx  =  |⟨x|η⟩((|)^2) dx				(6)

and 

	⟨x|η⟩ =   (∫_(-∞)^∞) ⟨x|k⟩ ⟨k|η⟩ dk  = C(∫_(-∞)^∞)exp[-((k^2)/(a^2))] ⟨x|k⟩ dk  			(7)
	
with 

	⟨x|k⟩ =(1/(Sqrt[2 π])) exp[ i k x]				(8)
	
We can combine (6) - (8) and calculate 

	(Subscript[⟨x⟩, η]) = (∫_(-∞)^∞) x (Subscript[P, η])(x) dx  					(9)
	
	(Subscript[⟨(x^2)⟩, η]) = (∫_(-∞)^∞) (x^2) (Subscript[P, η])(x) dx  				(10)
	
	(Subscript[σ, η])(x) = (Sqrt[(Subscript[⟨(x^2)⟩, η])-((Subscript[⟨x⟩, η])^2)]) 				(11)
Using these results we show that 

	(Subscript[σ, η])(x) (Subscript[σ, η])(p)  = (ℏ/2)

Cell 1: The Gaussian wave packet is defined by (in momentum representation)

	|η⟩ ≡ C(∫_(-∞)^∞)exp[-((k^2)/(a^2))] |k⟩ dk 				(1)
	
where C and a are constants and a has the same units as the wave vector (1/length). The probability distribution that a particle in this state has the momentum ℏk is

	(Subscript[P, η])(k) =|⟨k|η⟩((|)^2) = B exp[ - ((2(k^2))/(a^2))]			(2)
	
Determine B from:

	(∫_(-∞)^∞)(Subscript[P, η])(k)dk = 1 					(3)
	
After that calculate ⟨p⟩, ⟨(p^2)⟩ and (Subscript[σ, η])(p)

Define the function pη[k] which gives (Subscript[P, η])(k)

pη[k_,a_]:=B Exp[- 2 k^2/a^2]

Here B ≡ (C^2) .  Since C is unknown I might as well denote (C^2) as B and calculate it from the normalization condition 
The equation 

(∫_(-∞)^∞)(Subscript[P, η])(k)dk = 1 	

is:

eq=(Integrate[pη[k,a],{k,-Infinity,Infinity},Assumptions->Re[a^2]>0]//PowerExpand)==1

Solve for B

sol=Solve[eq,B] //Flatten

The probability distribution for k is therefore

Pη[k_,a_]:=pη[k,a]/.sol

Pη[k,a]

The mean momentum is 

	 (Subscript[⟨p⟩, η]) = ℏ (∫_(-∞)^∞)k (Subscript[P, η])(k)dk

pav= ℏ Integrate[k Pη[k,a],{k,-Infinity,Infinity},Assumptions->Re[a^2]>0]

The mean value of momentum squared is 

	(Subscript[⟨(p^2)⟩, η]) = (ℏ^2) (∫_(-∞)^∞)(k^2) (Subscript[P, η])(k)dk

psav= ℏ^2 Integrate[k^2 Pη[k,a],{k,-Infinity,Infinity},Assumptions->Re[a^2]>0] //PowerExpand

The square root of (Subscript[⟨(p^2)⟩, η]) is

σpη=Sqrt[psav]//PowerExpand

Cell 2: Evaluate 
	
	(Subscript[P, η])(x) = |⟨x|η⟩((|)^2)

where

⟨x|η⟩ =C(∫_(-∞)^∞)exp[-((k^2)/(a^2))] ⟨x|k⟩ dk=
F (∫_(-∞)^∞)dk  exp[-((k^2)/(a^2))] exp[i k x] 

I determine the constant factor in (Subscript[P, η])(x) or ⟨x|η⟩ by using

	(∫_(-∞)^∞)(Subscript[P, η])(x)dx = (∫_(-∞)^∞)|⟨x|η⟩((|)^2)dx = 1

Create the function ψ[x] that evaluates

 ψ[x] = (∫_(-∞)^∞)dk  exp[-((k^2)/(a^2))] exp[i k x]

Clear[x,A,k,a,ψ]
ψ[x_]= Integrate[Exp[-k^2/a^2 ]Exp[I k x],{k,-Infinity,Infinity},Assumptions->Re[a^2]>0]//PowerExpand

The probability is 

	(Subscript[P, η])(x)  = C | ψ[x]((|)^2)

and the constant C is determine from the normalization 

 (∫_(-∞)^∞)(Subscript[P, η])(x)dx  = 1

so=Solve[C Integrate[ψ[x]^2,{x,-Infinity,Infinity},Assumptions->Re[a^2]>0]==1,C]//Flatten//PowerExpand

Define the function pη[ x ] which gives the probability distribution   (Subscript[P, η])(x)

pη=.
pη[x_]:= C ψ[x]^2/.so

pη[y]

Check that it is normalized:

Integrate[pη[x] , {x,-Infinity,Infinity},Assumptions->Re[a^2]>0]//PowerExpand

Plot it:

Plot[ pη[x]^2 /.a->1 ,{x,-4,4}]

Cell 3: Calculate the mean value of x and the standard deviation of x in the state η[x]

The average value of x:

(Subscript[⟨x⟩, η]) = (∫_(-∞)^∞)x (Subscript[P, η])  (x)dx

avx=Integrate[ x pη[x],{x,-Infinity,Infinity},Assumptions->Re[a^2]>0]

The average value of (x^2)

(Subscript[⟨(x^2)⟩, η]) = (∫_(-∞)^∞)((x )^2)(Subscript[P, η])  (x)dx

avx2 =Integrate[pη[x] x^2 ,{x,-Infinity,Infinity},Assumptions->Re[a^2]>0]//PowerExpand

Calculate (Subscript[σ, ϕ])(x) which is the standard deviation of position in the state ϕ[x]

σxη=Sqrt[ avx2-avx^2]//PowerExpand

Summarize the results:

(Sqrt[(Subscript[⟨(p^2)⟩, η])])=((a ℏ)/2)
and 

(Sqrt[(Subscript[⟨(x^2)⟩, η])]) =(1/a)

Therefore 
(Sqrt[(Subscript[⟨(p^2)⟩, η])])(Sqrt[(Subscript[⟨(x^2)⟩, η])]) = (ℏ/2)
This is Heisenberg uncertainty principle (remember that the average values of k and x are zero.

Section 4:  Use the general theory to determine the uncertainty of position and momentum.

I want to see whether the general equation for the uncertainty in the measurement of two observables A and B whose operators do not commute gives the same result as Heisenberg uncertainty result if A is position, B is momentum and the state |ψ⟩ is the Gaussian wave packet 

	ψ(x) ≡ ⟨x|ψ⟩ = (1/(((2π))^(1/4)))(Sqrt[a ]) Exp[ -(((a^2)(x^2))/4)]
 For this calculation I need to evaluate the average values of the commutator

	⟨ψ| [(OverHat[p]), (OverHat[x])] |ψ⟩ = (∫_(-∞)^∞)dx (SuperStar[ψ(x)]) [ (ℏ/i) (∂/(∂x)), x ]  ψ(x) 
	and the anticommutator
	⟨ψ| {(OverHat[p]), (OverHat[x])}|ψ⟩ = (∫_(-∞)^∞)dx (SuperStar[ψ(x)]) { (ℏ/i) (∂/(∂x)), x }  ψ(x) 
 where

[A,B]=AB-BA
and 
{A,B}=AB+BA

Create a function that generates the Gaussian wave packet in coordinate space (this is ⟨x|ψ⟩ from lecture notes)

ψ=.
ψ[x_]:= (1/(2 Pi)^(1/4))Sqrt[a] Exp[- a^2 x^2/4]

Test it:

ψ[y]

Calculate the commutator ⟨ψ|[(OverHat[p]),(OverHat[x])]|ψ⟩ where (OverHat[p]) and (OverHat[x]) are the momentum and the position operators
⟨ψ | [(OverHat[p]), (OverHat[x])] |ψ⟩ = (∫_(-∞)^∞)ψ(x) [(ℏ/i) (∂/(∂x)), x ]ψ(x)dx =(∫_(-∞)^∞)ψ(x) [(ℏ/i) (∂/(∂x)) x ψ(x)-x(ℏ/i) (∂/(∂x))ψ(x)]dx

compxψ=Integrate[ ψ[x](( ℏ/I) D[x ψ[x],x]-(ℏ/I)x D[ψ[x],x]),{x,-Infinity,Infinity},Assumptions->Re[a^2]>0]//PowerExpand

Calculate the average value of the anticommutator

⟨ψ | {(OverHat[p]), (OverHat[x])} |ψ⟩ = (∫_(-∞)^∞)ψ(x) { (ℏ/i) (∂/(∂x)), x }ψ(x)dx
where
	
	 { (ℏ/i) (∂/(∂x)), x }ψ(x) = (ℏ/i) (∂/(∂x))(x ψ(x) + x(ℏ/i) (∂/(∂x)) ψ(x)

anticompxψ=Integrate[ ψ[x](( ℏ/I) D[x ψ[x],x]+(ℏ/I)x D[ ψ[x],x]),{x,-Infinity,Infinity},Assumptions->Re[a^2]>0]

Check each term separately to see why they cancel

Integrate[ ψ[x]( ℏ/I) D[x ψ[x],x],{x,-Infinity,Infinity},Assumptions->Re[a^2]>0]//PowerExpand

Integrate[ ψ[x] (ℏ/I)x D[ ψ[x],x] ,{x,-Infinity,Infinity},Assumptions->Re[a^2]>0]//PowerExpand

Heisenberg says that (((Subscript[σ, ψ])(x))^2) (((Subscript[σ, ψ])(p))^2) > ((ℏ^2)/4). Calculate the right hand side of the equation below

	(((Subscript[σ, ψ])(x))^2) (((Subscript[ϕ, ψ])(p))^2)≥ (1/4)| ⟨ψ|[p,x]|ψ⟩((|)^2) + (1/4)⟨ψ|{p,x}|ψ⟩((|)^2)

Abs[compxψ^2/4]+ anticompxη^2/4

We get the same result as Heisenberg.