Mathematica program for Problem 12.12

  date : 28/11/2013
  author: Jorge C. Romão
  email : jorge.romao@tecnico.ulisboa.pt

1. Definitions

The Radial Functions R[n,L,r]

R = Function[{n,L,r},(2/(n a0))^(3/2) Sqrt[(n-L-1)!/(2 n)/((n +L)!)]
	     (2 r/(n a0))^L LaguerreL[n-L-1,2 L+1, 2 r/(n a0)] E^(-r/(n a0))];

The radial momentum

pr = Function[f, -I hbar ( D[f, r] + f/r)];

The operators p_r^2 R_{nl} and p_r^4 R_{nl}

p2RnL=Function[{n,L,r}, pr[pr[R[n,L,r]]] +1/r^2 hbar^2 L(L+1) R[n,L,r]];

p4RnL=Function[{n,L,r}, pr[pr[pr[pr[R[n,L,r]]]]] 
		 + 2/r^2 hbar^2 L(L+1) pr[pr[R[n,L,r]]]
	       +1/r^4 hbar^4 (L(L+1))^2 R[n,L,r]];

The Hamiltonian of the H Atom

H0=Function[{f,L},pr[pr[f]]/(2 m)+hbar^2/(2 m r^2) L(L+1) f-alf hbar c/r f];

Functions to Check Hermiticity of p_r, p_r^2, p_r^3 and p_r^4 (two definitions)

Checkpr=Function[{n,m,L},aux=r^2 (R[n, L, r] pr[R[m, L, r]] + 
	R[m, L, r] pr[R[n, L, r]]) // Simplify ; 
		 Integrate[aux,{r,0,Infinity},Assumptions->a0>0]]; 

Checkpr2=Function[{n,m,L},aux=r^2 (R[n, L, r] pr[pr[R[m, L, r]]] - 
	R[m, L, r] pr[pr[R[n, L, r]]]) // Simplify ; 
		 Integrate[aux,{r,0,Infinity},Assumptions->a0>0]];  

Checkpr3=Function[{n,m,L},aux=r^2 (R[n, L, r] pr[pr[pr[R[m, L, r]]]] + 
	R[m, L, r] pr[pr[pr[R[n, L, r]]]]) // Simplify ; 
		 Integrate[aux,{r,0,Infinity},Assumptions->a0>0]];  

Checkpr4A=Function[{n,m,L},aux=r^2 (R[n, L, r] pr[pr[pr[pr[R[m, L, r]]]]] -
	R[m, L, r] pr[pr[pr[pr[R[n, L, r]]]]]) // Simplify ; 
		 Integrate[aux,{r,0,Infinity},Assumptions->a0>0]];  

Checkpr4B=Function[{n,m,L},aux=r^2 (pr[pr[R[n, L, r]]] pr[pr[R[m, L, r]]] -
	R[m, L, r] pr[pr[pr[pr[R[n, L, r]]]]]) // Simplify ; 
		 Integrate[aux,{r,0,Infinity},Assumptions->a0>0]];

Question 1

aux=r pr[f[r]] -pr[ r f[r]] //Simplify ;
aux = ℏ aux/f[r]/hbar ;
Print["[r,p_r] = ",aux]

aux=pr[pr[f[r]]]/hbar^2 //Simplify;
aux=-ℏ^2 (-aux);
Print["p_r^2 f(r) = ",aux]

Questions 2 and 3

aux110=Checkpr[1,1,0];
aux210=Checkpr[2,1,0];
aux320=Checkpr[3,2,0];
aux321=Checkpr[3,2,1];

Print["<10| p_r 10> - <p_r 10|10> = ",aux110]
Print["<20| p_r 10> - <p_r 20|10> = ",aux210]
Print["<30| p_r 20> - <p_r 30|20> = ",aux320]
Print["<31| p_r 21> - <p_r 31|21> = ",aux321]

Question 4

aux110=Checkpr2[1,1,0];
aux210=Checkpr2[2,1,0];
aux320=Checkpr2[3,2,0];
aux321=Checkpr2[3,2,1];

Print["<10| p_r^2 10> - <p_r^2 10|10> = ",aux110]
Print["<20| p_r^2 10> - <p_r^2 20|10> = ",aux210]
Print["<30| p_r^2 20> - <p_r^2 30|20> = ",aux320]
Print["<31| p_r^2 21> - <p_r^2 31|21> = ",aux321]

Question 5

aux110=Checkpr3[1,1,0];
aux110=ℏ^3 aux110/hbar^3;
aux210=Checkpr3[2,1,0];
aux210=ℏ^3 aux210/hbar^3;
aux320=Checkpr3[3,2,0];
aux320=ℏ^3aux320/hbar^3;
aux321=Checkpr3[3,2,1];
aux432=Checkpr3[4,3,2];

Print["<10| p_r^3 10> - <p_r^3 10|10> = ",aux110]
Print["<20| p_r^3 10> - <p_r^3 20|10> = ",aux210]
Print["<30| p_r^3 20> - <p_r^3 30|20> = ",aux320]
Print["<31| p_r^3 21> - <p_r^3 31|21> = ",aux321]
Print["<42| p_r^3 32> - <p_r^3 42|32> = ",aux432]

Question 6

aux110=Checkpr4A[1,1,0];
aux210=Checkpr4A[2,1,0];
aux320=Checkpr4A[3,2,0];
aux321=Checkpr4A[3,2,1];
aux432=Checkpr4A[4,3,2];

Print["<10| p_r^4 10> - <p_r^4 10|10> = ",aux110]
Print["<20| p_r^4 10> - <p_r^4 20|10> = ",aux210]
Print["<30| p_r^4 20> - <p_r^4 30|20> = ",aux320]
Print["<31| p_r^4 21> - <p_r^4 31|21> = ",aux321]
Print["<42| p_r^4 32> - <p_r^4 42|32> = ",aux432]

Question 7

aux110=Checkpr4B[1,1,0];
aux110=ℏ^4 aux110/hbar^4;
aux210=Checkpr4B[2,1,0];
aux210=ℏ^4 aux210/hbar^4;
aux320=Checkpr4B[3,2,0];
aux320=ℏ^4 aux320/hbar^4;
aux321=Checkpr4B[3,2,1];
aux432=Checkpr4B[4,3,2];

Print["< p_r^2 10| p_r^2 10> - <10|p_r^4 10> = ",aux110]
Print["< p_r^2 20| p_r^2 10> - <20|p_r^4 10> = ",aux210]
Print[" p_r^2 <30| p_r^2 20> - <30|p_r^4 20> = ",aux320]
Print["< p_r^2 31| p_r^2 21> - <31|p_r^4 21> = ",aux321]
Print["< p_r^2 42| p_r^2 32> - <42|p_r^4 32> = ",aux432]

Question 8

aux100=r^2 pr[pr[R[1,0,r]]] pr[pr[R[1,0,r]]];
res100=-1/(8 m^3 c^2 )Integrate[aux100,{r,0,Infinity},Assumptions->a0>0];
res100=res100 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <p_r^2 100|p_r^2 100> = ",res100]

aux100=r^2 R[1,0,r] pr[pr[pr[pr[R[1,0,r]]]]];
res100=-1/(8 m^3 c^2 )Integrate[aux100,{r,0,Infinity},Assumptions->a0>0];
res100=res100 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <100|p_r^4 100> = ",res100]

aux300=r^2 pr[pr[R[3,0,r]]] pr[pr[R[3,0,r]]];
res300=-1/(8 m^3 c^2 )Integrate[aux300,{r,0,Infinity},Assumptions->a0>0];
res300=res300 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <p_r^2 300|p_r^2 300> = ",res300]

aux300=r^2 R[3,0,r] pr[pr[pr[pr[R[3,0,r]]]]];
res300=-1/(8 m^3 c^2 )Integrate[aux300,{r,0,Infinity},Assumptions->a0>0];
res300=res300 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <300|p_r^4 300> = ",res300]

Question 9

aux210=r^2 pr[pr[R[2,1,r]]] pr[pr[R[2,1,r]]];
res210=-1/(8 m^3 c^2 )Integrate[aux210,{r,0,Infinity},Assumptions->a0>0];
res210=res210 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <p_r^2 210|p_r^2 210> = ",res210]

aux210=r^2 R[2,1,r] pr[pr[pr[pr[R[2,1,r]]]]];
res210=-1/(8 m^3 c^2 )Integrate[aux210,{r,0,Infinity},Assumptions->a0>0];
res210=res210 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <210|p_r^4 210> = ",res210]

aux210=r^2 p2RnL[2,1,r] p2RnL[2,1,r];
res210=-1/(8 m^3 c^2 )Integrate[aux210,{r,0,Infinity},Assumptions->a0>0];
res210=res210 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <p^2 210|p^2 210> = ",res210]

aux210=r^2 R[2,1,r] p4RnL[2,1,r];
res210=-1/(8 m^3 c^2 )Integrate[aux210,{r,0,Infinity},Assumptions->a0>0];
res210=res210 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <210|p^4 210> = ",res210]

aux320=r^2 pr[pr[R[3,2,r]]] pr[pr[R[3,2,r]]];
res320=-1/(8 m^3 c^2 )Integrate[aux320,{r,0,Infinity},Assumptions->a0>0];
res320=res320 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <p_r^2 320|p_r^2 320> = ",res320]

aux320=r^2 R[3,2,r] pr[pr[pr[pr[R[3,2,r]]]]];
res320=-1/(8 m^3 c^2 )Integrate[aux320,{r,0,Infinity},Assumptions->a0>0];
res320=res320 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <320|p_r^4 320> = ",res320]

aux320=r^2 p2RnL[3,2,r] p2RnL[3,2,r];
res320=-1/(8 m^3 c^2 )Integrate[aux320,{r,0,Infinity},Assumptions->a0>0];
res320=res320 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <p^2 320|p^2 320> = ",res320]

aux320=r^2 R[3,2,r] p4RnL[3,2,r];
res320=-1/(8 m^3 c^2 )Integrate[aux320,{r,0,Infinity},Assumptions->a0>0];
res320=res320 /. a0->hbar/(m c α) // Simplify;
Print["-1/(8 m^3 c^2 ) <320|p^4 320> = ",res320]

Question 10

aux=H0[pr[pr[f[r]] + 1/r^2 L (L + 1)], L] - pr[pr[H0[f[r], L]]] - 
  1/r^2 L (L + 1) H0[f[r], L] /. L->0// Simplify ;
Print["[H_0,p_r^2] f(r) = ",aux]