WorkBookQM. 18  The hydrogen atom

This is the corrected version of the WorkBookQM.18 given on the CD with Physical Chemistry Quantum Mechanics

WorkBookQM.18.1.   The Coulomb interaction energy

Define a function that gives Coulomb interaction energy between the proton and the electron in eV, when r is given in Å
Use
	V(r) = ((e^2)/(4 π (Subscript[ϵ, 0]) (r^2))) 
	
Use SI units to get V in joule. The conversion from meter to Å is

	1 m = (10^(-10))Å
and 
	1 joule = 6.242*(10^18)eV

Off[Remove::rmnsm]
Remove["Global`*"]

V[r_ (*in A*) ] :=  
      With[{e = 
            1.602177*10^-19  (*Coulomb*) , 
ϵ0 = 
            8.854187*10^-12  (*C^2 /(N m^2)*) , 
convertJtoeV = 6.242*10^18, 
AtoM = 
            10^-10  (*converts A to meter*) }, 
-convertJtoeV*
          e^2/(4 Pi ϵ0 r*AtoM )]

Check that the function works:

V[1]

Plot[V[r], {r, 0, 30}, PlotRange -> {-10, 0}, 
      AxesLabel -> {"r Å", "V(r) eV"}]

Note two things: the interaction energy goes to -∞ when r →0 and the interaction is very long ranged. It is not zero even when r = 30 Å

WorkBookQM.18. 2.  The energy of the hydrogen atom

The energy is given by :

\(TraditionalForm\`\(\(	\)\(E\_n\)\)\)= \(TraditionalForm\`μe\^4\/\(2\ \(\((4  πϵ\_0)\)\^2\) ℏ\^2\)\)(Z/(n^2)) 	n=1, 2, 3, ......		SI units

with the reduced mass 

	μ =(((Subscript[m, e]) (Subscript[m, N]))/((Subscript[m, e]) + (Subscript[m, N])))
	
The values of the constants are  
	me = 9.109389*10^-31  kg
mp = 1.672623*10^-27 kg
e = 
          1.602177*10^-19  Coulomb, 
ϵ0 = 
          8.854187*10^-12  Coulomb (Newton^(-1)) 
            (meter^(-2)), 
ℏ = 
          1.054572*10^-34  J sec, 
1 joule = 
          10^7*6.2420*10^11  eV

Define a function that gives the energy of the hydrogen atom in eV

Needs["Utilities`CleanSlate`"]
CleanSlate[]
Remove["Global`*"]

Cell 1: Calculate the energy and make an energy level diagram for H

\(\(\(\(Clear[energyH];\)

    energyH[
        n_\ ] := 
Module[{me = 
            9.109389*10^\(-31\)\  (*kg*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ mp = 
            1.672623*10^\(-27\) (*kg*) , 
	\ \ \ \ \ e = 
            1.602177*10^\(-19\) (*Coulomb*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ ϵ0 = 
            8.854187*10^\(-12\)\  (*C\ \(Newton\^\(-1\)\) 
              meter\^\(-2\)*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ℏ = 
            1.054572*10^\(-34\)\  (*J\ sec*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ JToeV = 10^7*6.2420*10^11, 
          R\ , μ}, 
μ = 
          me*mp/\((me + mp)\); 
R = 
          JToeV*μ*
            e^4/\((2\ \((4\ Pi\ ϵ0\ ℏ)\)^2)\); 
\(-R\)/n^2]\)\(	\ \ \ \ \ \ \)\)\)

Test the function:

energyH[n]

Make a table with the energies for various values of n:

\(Table[{n, energyH[n]}, {n, 1, 30}] // 
      TableForm[#, 
          TableHeadings -> {None, {"n", \*""(Subscript[E, n]) eV""}}] &\)

Note that a n increases the energies get closer and closer together and closer and closer to zero.  The states having high values of n are called Rydberg states and. All atoms and all molecules have Rydberg states.

Make an energy level diagram. Create horizontal lines at coordinates equal to the energy of the atom.

energies = 
      Table[ln[n] = Line[{{0, energyH[n]}, {1, energyH[n]}}], {n, 1, 
          10}]

Show[Graphics[#] & /@ energies , 
        Axes -> True];

Cell 2:  The effect of mass

Calculate the energy in the case when the nuclear mass is infinity (in this case the reduced mass equals the electron mass)

Clear[energyHInfNucl];
\(\(\(energyHInfNucl[
        n_\ ]\)\(:=\)\(
\)\(Module[{me = 
            9.109389*10^\(-31\)\  (*kg*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ mp = 
            1.672623*10^\(-27\) (*kg*) , 
	\ \ \ \ \ e = 
            1.602177*10^\(-19\) (*Coulomb*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ ϵ0 = 
            8.854187*10^\(-12\)\  (*C\ \(Newton\^\(-1\)\) 
              meter\^\(-2\)*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ℏ = 
            1.054572*10^\(-34\)\  (*J\ sec*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ JToeV = 10^7*6.2420*10^11, R}, 
μ = 
          me\ ; 
R = 
          JToeV*μ*
            e^4/\((2\ \((4\ Pi\ ϵ0\ ℏ)\)^2)\); 
\(-R\)/
          n^2]\)\(\ \)\( (*gives\ the\ energy\ in\ eV*) \)\(	\ \ \ \ \ \ \)\)\)

energyHInfNucl[n]

energyH[n]

Calculate the energy of the photon causing the transition
	{n=1}→{n=2}
for infinite nuclear mass and for the real mass. For the real mass:

eH = energyH[2] - energyH[1]

For infinite nuclear mass:

eInf = energyHInfNucl[2] - energyHInfNucl[1]

The difference is:

eInf - eH

This energy difference is detectable in the laboratory.

WorkBookQM.18. 3. The radial wave functions R\(TraditionalForm\`\_\(n, ℓ\)\)(r)  and (Subscript[Φ, n,  ℓ])(r): the definitions and expressions for various values of n and ℓ

The definition of the radial wave function is:


(Subscript[R, n, ℓ])(r)= (((((2 Z)/(n a))))^(3/2))\(\(Φ\_\(n, ℓ\)\)(ρ)\) 							(1)

(Subscript[Φ, n, ℓ])(r) = (Sqrt[(((n - ℓ - 1)!)/((n + ℓ)! 
                2 n))])\(\((\(2\ Z\ ρ\)\/n)\)\^\(\(\ \)\(ℓ\)\)\)  Exp[-(( Z ρ)/n) ]\(\(L\_\(n - ℓ - 1\)\%\(2  ℓ + 1\)\)(\(2\ Z\ ρ\)\/n)\)			(2)

with 
	ρ=(( r)/( a))									(3)
and 
   	a=\(TraditionalForm\`\(\(ℏ\^2\)(4\ π\ ϵ\_0)\)\/\(μ\ e\^2\)\) 								(4)
In many books the general equations (1)-(3) are erroneous: they do not generate the examples given in the same books (i.e. the expression for (Subscript[R, 2, 1])(r) obtained from their general formula do not agree with the function (Subscript[R, 2, 1])(r) given in the same book).

Define a function that generates (Subscript[Φ, n, ℓ])(ρ)

Unprotect[Φ];
Clear[Φ];  Unprotect[a, Z]; Clear[a, Z];  
    Protect[a, Z];
\(\(Φ[n_Integer, ℓ_Integer]\)[ρ_] := 
      Module[{const}, 

        const = Sqrt[\(\(\ \)\(\((n - ℓ - 1)\)!\)\)\/\(\(\((n + ℓ)\)!\)\ 2\ n\)]; 
const\ \((2\ Z\ ρ/
                  n)\)^ℓ\ \ LaguerreL[n - ℓ - 1, 
            2\ ℓ + 
              1, \(2\ Z\ ρ\)\/n]\ Exp[\(-\(\(\(\ \)\(Z\ ρ\)\)\/n\)\)]]\)
Protect[Φ];

Check that Φ works properly:

Φ[3, 2][ρ]

Define a function that generates (Subscript[R, n, ℓ])(r)

Unprotect[R];
Clear[R];
R[n_, ℓ_][
        r_] := (2  Z/ (n a))^(3/2) Φ[
            n, ℓ][r/a]
Protect[R];

See if (Subscript[R, n, ℓ]) (r)works

R[1, 0][r]

I made below a table with the radial function (Subscript[R, n, ℓ])(r)  and I compared my results with those given by W. H. Flygare, Molecular Structure and Dynamics, Prentice Hall, Inc., Englewood Cliffs, New Jersey, 1978, p.144. My functions agree with Flygare's. They also agree with the functions given by Pauling and Wilson for specific values of n and ℓ. However their general formula  Eq. 20-11 (in the Dover edition of Pauling and Wilson) does not agree with the examples given by them.

\(\(Table[{n, ℓ, \((\(\(R[n, ℓ]\)[r] // Simplify\) // 
                Cancel\ )\)
}, {n, 1, 4}, {ℓ, 0, 
            n - 1}] // Flatten[#, 1] &\) // 
      TableForm[#, 
          TableHeadings -> {None, {"n", "ℓ", \*""(Subscript[R, n, ℓ])(r)""}}] &\)

To compare to Pauling and Wilson I use their variable
	ρ ≡((2 Z r)/(n a))
Make a table of (Subscript[R, n])ℓ(ρ) with this variable

\(\(Table[{n, ℓ, \((\(R[n, ℓ]\)[r] /. 
                  r -> ρ\ n\ a/\((2\ Z)\))\) // 
              Cancel\ 
}, {n, 1, 4}, {ℓ, 0, 
            n - 1}] // Flatten[#, 1] &\) // 
      TableForm[#, 
          TableHeadings -> {None, {"n", "ℓ", \*""(Subscript[R, n, ℓ])(r)""}}] &\)

The functions in the table are the same as those of Pauling and Wilson, p. 135 (Dover edition)

In the textbook I prefer to use the dimensionless function (Subscript[Φ, n, ℓ])(ρ) which I display below (for Z=1):

\(\(\(\(radf = \ \ Table[{n, ℓ, \ \((\ \(Φ[
                          n, ℓ]\)[ρ] /. Z -> 1\ \  // 
                    Simplify)\)\  // Cancel\ }, {n, 1, 4}, {\ ℓ, 0, 
              n - 1}] // Flatten[#, 1] &;\)

    TableForm[radf, TableSpacing -> {5, 3}, 
      TableHeadings -> {None, {"n", "ℓ", "Φ[n,ℓ][ρ]"}}, 
TableSpacing -> {3, 5}]\)\(\ \)\)\)

WorkBookQM.18. 4. Test orthonormality 

	\(TraditionalForm\`\(\(∫\_0\%∞\( 
                r\^2\) \(\(\(R\_\(n, ℓ\)\)(
                  r)\)\^*\) \(R\_\(n', ℓ'\)\) dr\)\(\ \ \)\(=\)\)\) (Subscript[δ, n, n'])

Copy the function R[n,ℓ][r], which defines the radial wave function, from WorkBookQM.18.3 :

Off[
      Remove::rmptc]   
Off[
      Remove::relex]   
Needs["Utilities`CleanSlate`"]
CleanSlate[]
Remove["Global`*"]
Off[General::spell1]
Off[General::spell]

Redefine the functions Φ and R:

Unprotect[Φ];
Clear[Φ];  Unprotect[a, Z]; Clear[a, Z];  
    Protect[a, Z];
\(\(Φ[n_Integer, ℓ_Integer]\)[ρ_] := 
      Module[{const}, 

        const = Sqrt[\(\(\ \)\(\((n - ℓ - 1)\)!\)\)\/\(\(\((n + ℓ)\)!\)\ 2\ n\)]; 
const\ \((2\ Z\ ρ/
                  n)\)^ℓ\ \ LaguerreL[n - ℓ - 1, 
            2\ ℓ + 
              1, \(2\ Z\ ρ\)\/n]\ Exp[\(-\(\(\(\ \)\(Z\ ρ\)\)\/n\)\)]]\)
Protect[Φ];

Define a function that generates (Subscript[R, n, ℓ])(r)

Unprotect[R];
Clear[R];
R[n_, ℓ_][
        r_] := (2  Z/ (n a))^(3/2) Φ[
            n, ℓ][r/a]
Protect[R];

Test that  

	(∫_0^∞) 
                (r^(2 )) (Subscript[R, n,  ℓ])(
              r) dr  = 1

meaning that the radial functions are normalized

Integrate[Expand[r^2 R[2, 1][r]^2] // Evaluate, {r, 0, Infinity}, 
      GenerateConditions -> False]

Table[
      Integrate[Expand[r^2 R[n, ℓ][r]^2] , {r, 0, Infinity}, 
        GenerateConditions -> False] , {n, 1, 5}, {ℓ, 0, 
        n - 1}]

Now let me check the orthogonality:

	\(TraditionalForm\`\(\(∫\_0\%∞\( 
                r\^2\) \(\(R\_\(n, \ ℓ\)\)(
              r)\)\ \(R\_\(m, \ ℓ\)\) dr\)\(\ \)\(=\)\(\ \)\(δ\_\(n, m\)\)\(\ \)\)\)
	
 The integral must be zero if n differs from m and the two functions have the same value for ℓ.

Below I examine the case when n differs from m and ℓ = 0 in both functions.

totable = 
      Table[{n, m, 
            Integrate[
              Expand[r^2 R[n, 0][r]*R[m, 0][r]] , {r, 0, Infinity}, 
              GenerateConditions -> False]} , 
{n, 1, 
            3}, {m, n + 1, 4}] // Flatten[#, 1] &

TableForm[totable, 
      TableHeadings -> {None, {"n", "m", "integral"}}]

Below I calculate the integrals for different  principal quantum numbers in the two functions, but  with ℓ = 1 in both functions:

totable = 
      Table[{n, m, 
            Integrate[
              Expand[r^2 R[n, 1][r]*R[m, 1][r]], {r, 0, Infinity}, 
              GenerateConditions -> False]} , 
{n, 2, 
            4}, {m, n + 1, 4}] // Flatten[#, 1] &

TableForm[totable, 
      TableHeadings -> {None, {"n", "m", "integral"}}]

Integrals in which n≠n' and ℓ≠ℓ' do not have to be orthogonal; the orthogonality of the total wave function is ensured by the integral over the spherical harmonics. Below the principal quantum numbers n and m are different from each other and the ℓ  values are different (in one function ℓ = 1 and in the other ℓ = 0)

totable = 
      Table[{n,  m, 
            Integrate[
              Expand[r^2 R[n, 0][r]*R[m, 1][r]], {r, 0, Infinity}, 
              GenerateConditions -> False]} , 
{n, 2, 
            4}, {m, n + 1, 4}]  // Flatten[#, 1] &

TableForm[totable, 
      TableHeadings -> {None, {"n", "m", "integral"}}]

WorkBookQM.18.5.  A radial function for numerical calculations for hydrogen.

Define a function RHNum[n,ℓ][x] which gives the radial wave function for the H atom in a form appropriate for numerical calculations

Clear[RHNum]
RHNum[
          n_Integer /; 
            n > 0, ℓ_Integer /; ℓ ≥ 0][
        r_  ]:=Module[\({const, ρ, \ 
            me = 9.109389*10^\(-31\)\  (*kg*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ mp = 
              1.672623*10^\(-27\) (*kg*) , 
	\ \ \ \ \ e = 
              1.602177*10^\(-19\) (*Coulomb*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ ϵ0 = 
              8.854187*10^\(-12\)\  (*C\ \(Newton\^\(-1\)\) 
                meter\^\(-2\)*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ℏ = 
              1.054572*10^\(-34\)\  (*J\ sec*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ JToeV = 10^7*6.2420*10^11, Z = 1\ , μ, a}\),
 μ = me*mp/(me + mp);
a =  
              10^10*ℏ^2*(4 Pi ϵ0)/(μ  e^2);
If[ℓ < n, Null, 
              Print[" ℓ must be less than n"]]; 
const = 
              Sqrt[((( 2 Z/(n a))^3 (n - ℓ - 1)!)/((n + ℓ)! 2 n))];
ρ = ((2 Z r)/(n a));
const ρ^ℓ LaguerreL[
                n - ℓ - 1, 
                2 ℓ + 1, ρ] Exp[-(ρ/2)]]

Check that it works

RHNum[2, 1][r]

WorkBookQM.18.6.  Plots of (Subscript[R, n, ℓ])(r) and of the radial probability (r^2)(((Subscript[R, n, ℓ])(r))^2) as a function of r, for various values of n and ℓ.

Make plots of (Subscript[R, n, ℓ])(r)and of (r^2)(Subscript[R, n, ℓ])(r) as a function of r, for n = 1, ℓ = 0 and n = 2, ℓ = 0, 1

I plot the two functions for various values on n and ℓ. However, I use dynamic labelling (the values of n and ℓ are automatically updated by the computer in the labels of the axes) and I also changed the fonts because otherwise the labels would be hard to read.

Off[
      Remove::rmptc]   
Off[
      Remove::relex]   
Needs["Utilities`CleanSlate`"]
CleanSlate[]
Remove["Global`*"]
Off[General::spell1]
Off[General::spell]

Copy from WorkBookQM.18.5 the function RHNum[n,ℓ][r] which gives the radial wave function for the hydrogen atom

Unprotect[RHNum];
Clear[
        RHNum];
RHNum[
            n_Integer /; 
              n > 0, ℓ_Integer /; ℓ ≥ 0][
          r_  ]:=Module[\({const, ρ, \ 
              me = 9.109389*10^\(-31\)\  (*kg*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ mp = 
                1.672623*10^\(-27\) (*kg*) , 
	\ \ \ \ \ e = 1.602177*10^\(-19\) (*Coulomb*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ ϵ0 = 
                8.854187*10^\(-12\)\  (*C\ \(Newton\^\(-1\)\) 
                  meter\^\(-2\)*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ℏ = 
                1.054572*10^\(-34\)\  (*J\ sec*) , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ JToeV = 10^7*6.2420*10^11, Z = 1\ , μ, a}\),
 μ = me*mp/(me + mp);
a =  
                10^10*ℏ^2*(4 Pi ϵ0)/(μ  e^2);
If[ℓ < n, Null, 
                Print[" ℓ must be less than n"]]; 
const = 
                Sqrt[((( 2 Z/(n a))^3 (n - ℓ - 1)!)/((n + ℓ)! 2 n))];
ρ = ((2 Z r)/(n a));
const ρ^ℓ LaguerreL[
                  n - ℓ - 1, 
                  2 ℓ + 1, ρ] Exp[-(ρ/2)]];
Protect[RHNum];

The function plR[n, ℓ, rmax] plots (Subscript[R, n, ℓ])(r) versus r.   rmax determines the largest value of r in the plot. The function is calculated with RHNum which was defined in WorkBookQM.18.5, and gives the wave function for hydrogen.

Clear[plR]
plR[n_, ℓ_, rMax_]:= Plot[RHNum[n, ℓ][r] ,{r, 0, rMax},
PlotRange -> All,AxesLabel -> {StyleForm["r", 
                        FontSize -> 14, FontWeight -> "Bold", 
                        FontFamily -> "Geneva" ], 
                      StyleForm[
                        ToString[
                          StringForm[" R[``,``][r]", n, ℓ]], 
                        FontSize -> 14, FontWeight -> "Bold", 
                        FontFamily -> "Geneva" ]},DisplayFunction -> Identity ];

Let us look at the graph for n=2 and ℓ=1:

Show[plR[2, 1, 8] , 
      DisplayFunction -> $DisplayFunction]

The function plP[n, ℓ, rmax] plots (r^2) (((Subscript[R, n, ℓ])(r))^2) versus r.   rmax determines the largest value of r in the plot. The function is calculated with the help of RHNum, which was defined in WorkBookQM.18.5, and gives the probability function for hydrogen.

Clear[plP]
plP[n_, ℓ_, rMax_] :=  
      Plot[r^2 *RHNum[n, ℓ][r]^2 , {r, 0, 
          rMax}, 
PlotRange -> All, 
        AxesLabel -> {StyleForm["r", FontSize -> 14, 
              FontWeight -> "Bold", 
              FontFamily -> "Geneva" ], 
            StyleForm[
              ToString[StringForm["r^2  R[``,``][r]^2", n, ℓ]], 
              FontSize -> 14, FontWeight -> "Bold", 
              FontFamily -> "Geneva" ]}, 
        DisplayFunction -> Identity ]

Display the probability for n=2 and ℓ =1:

Show[plP[2, 1, 8], DisplayFunction -> $DisplayFunction]

Make a table with the plots of the two functions for n=1 and n=2 and the corresponding values of ℓ

tb12 = 
        Table[{plR[n, ℓ, 8], plP[n, ℓ, 8]}, {n, 1, 
              2}, {ℓ, 0, n - 1}] // Flatten[#, 1] &;

Show the plots. Use GraphicsArray to put all the graphs in one plot.

Show[
      GraphicsArray[tb12, 
        DisplayFunction -> $DisplayFunction]]

Make the same plots for n=3 and ℓ=0, 1, 2

tb3 = 
        Table[{plR[3, ℓ, 14], 
            plP[3, ℓ, 14]},  {ℓ, 0, 2}];

Show[
      GraphicsArray[tb3, 
        DisplayFunction -> $DisplayFunction]]

Make the same plots for n=4 and ℓ=0, 1, 2, 3

tb4 = 
        Table[{plR[4, ℓ, 22], 
            plP[4, ℓ, 22]},  {ℓ, 0, 3}];

Show[
      GraphicsArray[tb4, 
        DisplayFunction -> $DisplayFunction]]

Note various peaks in the probability that indicate that the electron seems to reside in  "shells".

WorkBookQM.18. 7. Mean position, mean position squared, mean Coulomb interaction energy, mean rotational energy and the mean kinetic energy.

Off[
      Remove::rmptc]   
Off[
      Remove::relex]   
Needs["Utilities`CleanSlate`"]
CleanSlate[]
Remove["Global`*"]
Off[General::spell1]
Off[General::spell]

Cell 0: The radial wave function

Unprotect[Φ];
Clear[Φ];  Unprotect[a, Z]; Clear[a, Z];  
    Protect[a, Z];
\(\(Φ[n_Integer, ℓ_Integer]\)[ρ_] := 
      Module[{const}, 

        const = Sqrt[\(\(\ \)\(\((n - ℓ - 1)\)!\)\)\/\(\(\((n + ℓ)\)!\)\ 2\ n\)]; 
const\ \((2\ Z\ ρ/
                  n)\)^ℓ\ \ LaguerreL[n - ℓ - 1, 
            2\ ℓ + 
              1, \(2\ Z\ ρ\)\/n]\ Exp[\(-\(\(\(\ \)\(Z\ ρ\)\)\/n\)\)]]\)
Protect[Φ];

Define a function that generates (Subscript[R, n, ℓ])(r)

Unprotect[R];
Clear[R];
R[n_, ℓ_][
        r_] := (2  Z/ (n a))^(3/2) Φ[
            n, ℓ][r/a]
Protect[R];

Check that R works

R[3, 2][r]

Φ does not work if n and ℓ are symbols

Φ[n, ℓ][r/a]

Φ works if n and ℓ have numerical values.

Φ[3, 2][(r/a)]

Cell 1: 	Mean position, mean position squared and the dispersion σ of the position 


	⟨ n, ℓ | r | n, ℓ ⟩ ≡ \(TraditionalForm\`∫\_0\%∞ 
          r\ \(\(\(R\_\(n, ℓ\)\)(
                r)\)\^*\) \(\(R\_\(n, ℓ\)\)(
            r)\) \(r\^2\) \(\(dr\)\(\ \ \ \ \ \ \ \ \ \ \)\)\)     		(1)
	

\(TraditionalForm\`\(\(	\)\(⟨r\^2⟩\  ≡ \ ⟨n, ℓ | r\^2 | 
              n, ℓ⟩\  = \ ∫\_0\%∞\( r\^2\) {r\^2\ \(\(\(R\_\(n, ℓ\)\)(
                      r)\)\^*\) \(\(R\_\(n, ℓ\)\)(
                  r)\)} dr\)\)\)	     	(2)


	σ = \(TraditionalForm\`\@\(\(⟨r\^2⟩ - ⟨r⟩\^2\)\/⟨r\^2⟩\)\)                                             			(3)

Calculate the mean position ⟨n, ℓ | r | n, ℓ⟩  (see Eq. 1) and ⟨(r^2)⟩ (see Eq. 2) for various values of n and ℓ; make a table with the results:

tbMeanr = 
      Table[{n, ℓ, rmean[n][ℓ] = 
              Integrate[r^3 R[n, ℓ][r]^2, {r, 0, Infinity}, 
                GenerateConditions -> False ]}, {n, 1, 4}, {ℓ, 
            0, n - 1}] // Flatten[#, 1] &

TableForm[tbMeanr, 
      TableHeadings -> {None, {"n", "ℓ", "\
⟨r⟩"}}]

Clear[rMean2]
tb = 
        Table[{n, ℓ,   rMean2[n][ℓ] = 
                Integrate[r^4 R[n, ℓ][r]^2, {r, 0, Infinity}, 
                  GenerateConditions -> 
                    False ] }, 
{n, 1, 4}, {ℓ, 
              0, n - 1}] // Flatten[#, 1] &;
\(TableForm[tb, 
      TableHeadings -> {None, {"n", "ℓ", \ \*""⟨(r^2)⟩""\ }}]\)

We see that, as expected, ⟨r⟩ anc \(TraditionalForm\`\(\(⟨r\^2⟩\)\(\ \)\)\)are proportional to the length scale a and (a^2)respectively and inverse proportional to Z or (Z^2)(respectively). Higher the nuclear charge closer is the electron to the nucleus.  Since the probability function is rather complicated the magnitude of σ is not meaningful.

Calculate numerical values of the same quantitites, for hydrogen. For this I need the wave function of the hydrogen atom, which I copy (see below) from WorkBookQM.18.5

tb = 
      Table[{n, ℓ, rmean[n][ℓ], rMean2[
                  n][ℓ], 
              Sqrt[rMean2[
                        n][ℓ] - rmean[
                          n][ℓ]^2]/rmean[
                    n][ℓ]} /. {Z -> 1,  
              a -> 0.529   }, {n, 1, 4}, {ℓ, 0, 
            n - 1}] // Flatten[#, 1] &

\(TableForm[tb, \ 
      TableHeadings -> {None, {"n", "ℓ", "⟨r⟩", \*""⟨(r^2)⟩"", \*""(Sqrt[((⟨(r^2)⟩ - ⟨r⟩^2)/⟨)r⟩^2])""}}]\)

If  n increases, at fixed ℓ, the mean radius increases;  if ℓ increases, at fixed n, the mean radius decreases.  That seems strange: the mean distance between the electron and the nucleus decreases as  the rotational energy is increased. This happens because of the manner in which the Coulomb and the centrifugal energies combine to give the total energy.

Cell 2:  The average 
	radial kinetic energy ⟨K⟩, 
	rotational energy \(TraditionalForm\`\(\(⟨E\_r⟩\)\(\ \)\)\)and 
	Coulomb energy ⟨C⟩ 
	

\(TraditionalForm\`\(\(⟨K⟩\)\(\ \)\(=\)\(\ \)\)\)⟨n,ℓ|-((ℏ^2)/(2 μ (r^2)))\(TraditionalForm\`\(∂\/∂r\) \(r\^2\) ∂\/∂r\)|n,ℓ⟩ = (∫_0^∞)(r^2)(Subscript[R, n, ℓ])(r) {-((ℏ^2)/(2 μ (r^2)))\(\(∂\/∂r\) \(r\^2\) ∂\/∂r\)  } (Subscript[R, n, ℓ])(r)dr	(1)


\(TraditionalForm\`\(\(⟨E\_r⟩\)\(\ \)\(=\)\(\ \)\)\)⟨n,ℓ| \(TraditionalForm\`\(\(ℓ(ℓ + 1)\) ℏ\^2\)\/\(2\ μ\ r\^2\)\)|n,ℓ⟩ = (∫_0^∞)(r^2)(Subscript[R, n, ℓ])(r) {-(((ℏ^2) ℓ(ℓ + 
                              1))/(2 μ (r^2)))  } (Subscript[R, n, ℓ])(r)dr			(2)


\(TraditionalForm\`\(\(⟨C⟩\)\(\ \)\(=\)\(\ \)\)\)⟨n,ℓ|- ((e^2)/(4 π (Subscript[ϵ, 0]) r))|n,ℓ⟩ = (∫_0^∞)(r^2)(Subscript[R, n, ℓ])(r) {-((e^2)/(4 π (Subscript[ϵ, 0]) r))  } (Subscript[R, n, ℓ])(r)dr			(3)

Calculate the mean kinetic energy: define the function KinR[n, ℓ][r] that gives 

	-((ℏ^2)/(2 μ (r^2)))\(\(∂\/∂r\) \(r\^2\) ∂\/∂r\)  (Subscript[R, n, ℓ])(r)

Clear[KinR];
KinR[n_, ℓ_][r_]:=(-ℏ^2/(2 μ))*(1/x^2)*D[x^2 D[R[n, ℓ][x], x],x]/.x -> r;

Calculate ⟨K⟩ from Eq. 1:

K =. 

    K[n_, ℓ_] := 
      Integrate[
        r^2 R[n, ℓ][r]*KinR[n, ℓ][r], {r, 0, 
          Infinity}, GenerateConditions -> False]

Check that it works

K[3, 1]

Calculate the mean "rotational energy":

Er =. 

    Er[n_, ℓ_] := 
      Integrate[
        r^2 R[n, ℓ][
              r]^2*ℏ^2 ℓ (ℓ + 
                1)/(2 μ r^2), {r, 0, Infinity}, 
        GenerateConditions -> False]

Check that it works

Er[3, 2]

Compare the numerical results to the equation

	\(TraditionalForm\`⟨E\_r⟩\  = ℏ^2\/\(2\ μ\)\) \(TraditionalForm\`\(ℓ\ \((ℓ + 1)\) Z\^2\)\/\(\(a\^2\) \(\(n\^3\)(ℓ\  + \ 1\/2)\)\)\)

ErFormula[
        n_, ℓ_] :=  ℏ^2 Z^2 ℓ (ℓ + 
              1)/(2 μ a^2 n^3 (ℓ + 
                  1/2))

ErFormula[3, 2]

The two calculations give the same result.

Calculate the mean Coulomb energy ⟨C⟩

c =. 

    c[n_, ℓ_] := 
      Integrate[
        r^2 R[n, ℓ][r]^2*(-1)*
          Z  e^2/(4 π (Subscript[ϵ, 0]) r), {r, 0, Infinity}, 
        GenerateConditions -> False]

Calculate the mean Coulomb  energy for n=3 and ℓ =2, to see if it works

c[3, 2]

Compare to the formula:

⟨C⟩ = - ((Ze^2)/(4 π (Subscript[ϵ, 0]))) (Z/(an^2))

cFormula =. 
cFormula[
        n_] := -Z^2 e^2/(4 π (Subscript[ϵ, 0]) a n^2)

cFormula[3]

Since all other equations are given in terms of a and μ and  do not contain (Subscript[ϵ, 0]) I eliminate (Subscript[ϵ, 0]) from this equation by using 
	(Subscript[ϵ, 0])= ((μ (e^2) a)/(4 π (ℏ^2)))
which is obtained from the definition of a in §5.

rule(Subscript[ϵ, 0]) =  (Subscript[ϵ, 0]) -> μ e^2 a/(4 π ℏ^2)

cFormula[n] /. rule(Subscript[ϵ, 0])

Make a function giving the total energy

energy =. 
energy[
        n_ ] := 
Module[{abc}, 
-μ*Z^2* e^4/(2  n^2 (4 Pi (Subscript[ϵ, 0]) ℏ)^2)]

energy[n]

Replace (Subscript[ϵ, 0]) with its expression:

energy[n] /. rule(Subscript[ϵ, 0])

Calculate the mean kinetic energy ⟨K⟩, the mean rotational energy ⟨(Subscript[E, r])⟩, the mean Coulomb energy ⟨C⟩ and the total energy E, for various values on n and ℓ and make a table with the results:

tb = 
        Table[{n,  ℓ,  K[n, ℓ], ErFormula[n, ℓ], 
              cFormula[n] /. rule(Subscript[ϵ, 0]), 
              energy[n] /. rule(Subscript[ϵ, 0])}, {n, 1, 4}, {ℓ, 0, 
              n - 1}] // Flatten[#, 1] &;
\(TableForm[tb, 
      TableHeadings -> {None, {"n", "ℓ", "⟨K⟩", \*""⟨(Subscript[E, r])⟩"", "⟨C⟩", "E"}}]\)

We can replace a with its expression (ℏ^2) 4 π (Subscript[ϵ, 0])/(μ (e^2))to get results in terms of the fundamental constants.

Note an interesting  result: the mean Coulomb energy is equal to the total energy divided by two. Now we know that the total kinetic energy must be equal to the total energy minus the Coulomb energy and it is also equal to the radial kinetic energy plus the rotational energy.  You can check that these equality are true, which indicates that our calculations are correct.

As a test let us check that the three mean energies add up to the total energy:

ttb = 
        Table[{n, ℓ, 
                K[n, ℓ] + (ErFormula[n, ℓ] /. 
                      rule(Subscript[ϵ, 0])) + (cFormula[n] /. 
                      rule(Subscript[ϵ, 0])), energy[n]} /. 
              rule(Subscript[ϵ, 0]), {n, 1, 4}, {ℓ, 0, n - 1}] // 
          Flatten[#, 1] &;

\(ttb // 
      TableForm[#, 
          TableHeadings -> {None, {"n", "ℓ", "sum", "E"}}] &\)

WorkBookQM.18. 9.  The total wave function (radial and angular)

The wave function of the one-electron atom is 

	(Subscript[Ψ, n,  ℓ,  m])= (Subscript[R, n,  ℓ])(
          r) (Subscript[Y, ℓ]^m)(θ, ϕ)					(1)
	
The radial wave function is:

(Subscript[R, n, ℓ])(r)= (((((2 Z)/(n a))))^(3/2))\(\(Φ\_\(n, ℓ\)\)(ρ)\) 						(2)

with

(Subscript[Φ, n, ℓ])(r) = (Sqrt[(((n - ℓ - 1)!)/(2 n (n + ℓ)!))])\(\((\(2\ Z\ ρ\)\/n)\)\^\(\(\ \)\(ℓ\)\)\)  Exp[-((Z ρ)/n)](Subscript[L, n - ℓ - 
                  1]^(2  ℓ + 1))(((2 Z ρ)/n))		(3)

You may find that Eq. 2  differs from what some books give. Many books give an erroneous formula for this quantity.
We have
	ρ=(( r)/( a))								(4)
and 
   	a=\(TraditionalForm\`\(\(ℏ\^2\)(4\ π\ ϵ\_0)\)\/\(μ\ e\^2\)\) 							(5)

Plan: 
1. Define the radial wave function (Subscript[R, n,  ℓ])
2. Define the total wave function (Subscript[Ψ, n,  ℓ,  m]) and display the formulae for various values of n, ℓ and m
3. Define the wave functions 1s, 2s, 3s, 2  (Subscript[p, x]), 2  (Subscript[p, y]), 2 (Subscript[p, z]),  etc. and display the formulae and plot them

Needs["Utilities`CleanSlate`"]
CleanSlate[]
Unprotect[a,  R,  Z, Φ];
Remove["Global`*"]
Off[General::spell1]
Off[General::spell]
Protect[a,  R,  Z, Φ];

Cell 1:  The radial wave function.  This is copied from WorkBookQM.18.3

The function Φ:

Unprotect[Φ];
Clear[Φ];  Unprotect[a, Z]; Clear[a, Z];  
    Protect[a, Z];
\(\(Φ[n_Integer, ℓ_Integer]\)[ρ_] := 
      Module[{const}, 

        const = Sqrt[\(\(\ \)\(\((n - ℓ - 1)\)!\)\)\/\(\(\((n + ℓ)\)!\)\ 2\ n\)]; 
const\ \((2\ Z\ ρ/
                  n)\)^ℓ\ \ LaguerreL[n - ℓ - 1, 
            2\ ℓ + 
              1, \(2\ Z\ ρ\)\/n]\ Exp[\(-\(\(\(\ \)\(Z\ ρ\)\)\/n\)\)]]\)
Protect[Φ];

Define a function that generates (Subscript[R, n, ℓ])(r)

Unprotect[R];
Clear[R];
R[n_, ℓ_][
        r_] := (2  Z/ (n a))^(3/2) Φ[
            n, ℓ][r/a]
Protect[R];

Test that the function works

R[3, 2][r]

Cell 2: Define the total wave function (Subscript[Ψ, n,  ℓ,  m])(
        r, θ,  ϕ)

Use Eq. 1 to define ψ[n, ℓ, m][r, θ, ϕ]

ψ[n_, ℓ_, m_][
        r_, θ_, ϕ_] := R[n, ℓ][r]*
        SphericalHarmonicY[ℓ, m, θ, ϕ]

See if it works:

ψ[3, 2, -1][r, θ, ϕ]

Make a table with (Subscript[Ψ, n,  ℓ,  m])(
        r, θ,  ϕ)  for various values of n, ℓ and m

tb = 
        Table[{n, ℓ, 
              m, ψ[n, ℓ, m][r, θ, ϕ]}, {n, 1, 
              3}, {ℓ, 0, n - 1}, {m, -ℓ, ℓ}] // 
          Flatten[#, 2] &;
\(TableForm[tb, 
      TableHeadings -> {None, {"n", "ℓ", "m", \*""\!\(ψ\_\(\(n\)\(,\)\(ℓ\)\(,\)\(m\)\(\\ \)\)\)(r,θ,ϕ""}}]\)

Note the following:
	If m ≠ 0 the function is complex (this means that its values are complex numbers, not that the function is complicated) 
	When m ≠ 0 the dependence on the angle ϕ is given by exp[±i m ϕ]
	When m = 0 the wave function is independent of ϕ.

Cell 3: The wave functions 

	ns(r,θ,ϕ) =  (Subscript[R, n, 0])(
          r) (Subscript[Y, 0]^0)(θ, ϕ) 
	        
This generates the functions 1s, 2s, 3s, .....

The spherical harmonic (Subscript[Y, 0]^0)(θ, ϕ)  is:

SphericalHarmonicY[0, 0, θ, ϕ]

This is independent of angles θ and ϕ

Create the function s[r, θ,ϕ][n] which gives 
	ns[r, θ, ϕ]= (Subscript[R, n, 0])(
          r) (Subscript[Y, 0]^0)(θ, ϕ)

Unprotect[s];
s =. 
s[r_, θ_, ϕ_][n_] := 
      Module[{}, 
R[n, 0][r]* 
          SphericalHarmonicY[0, 
            0, θ, ϕ]]
Protect[s];

Make a table with the functions 1s, 2s, 3s, .....

tbs = 
        Table[{ToString[
              StringForm["``s(r, θ,ϕ)", n]],  s[
                r, θ, ϕ][n] }, {n, 1, 
            5}] ;
TableForm[tbs ]

These s-orbitals are the (Subscript[R, n, 0])(r) functions multiplied with (Subscript[Y, 0]^0)(θ,ϕ) =1/(2 (Sqrt[π])).

All the ns(r,θ,ϕ) functions are real (take real values) and are independent of ϕ or θ. This means that they have spherical symmetry.  The probability of finding the electron is the same on any point on a sphere whose center is on the nucleus. The electron has no angular preference.

For making contour plots I want to have the ns function in cartesian coordinates.  Use the replacement rule

convCarts = { r -> Sqrt[x^2 + y^2 + z^2]}

For the s functions the conversion is easy since they depend on r only.

The clearest way of representing these functions is to make contour plots which show the function in one of the planes of the coordinate system.
Below I plot 1s(x, y, z=0).

ContourPlot[
        Evaluate[(s[r, θ, ϕ][1] /. 
                  convCarts) /. {a -> 0.5292, Z -> 1} /. 
            z -> 0], {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
        PlotPoints -> 200,  Contours -> 15, PlotRange -> All, 
        FrameLabel -> {"x", "y"},  
        PlotLabel -> "1s(x,y,z=0)"];

The spherical symmetry is evident in the graph.  In the regions that are whiter than the background the function takes positive valuesThe 2s, 3s and 4s functions look the same except that they extend further away from the nucleus.

Cell 4: The wave functions

	(Subscript[np, x])(r,  θ, ϕ) = (1/(Sqrt[2])){(Subscript[Y, 1]^(-1))(θ, ϕ)-(Subscript[Y, 1]^1)(θ, ϕ)} (Subscript[R, n,  1])(r)

The angular part of all  (Subscript[np, x]) function is the same for all n, and it is given below

Clear[pxAng]
pxAng[θ, ϕ] := (SphericalHarmonicY[
                1, -1, θ, ϕ] - 
              SphericalHarmonicY[1, 1, θ, ϕ])/
          Sqrt[2]//ComplexExpand

pxAng[θ, ϕ]

To convert to cartesian coordinates I use:

	cos(θ)= (z/r)
	sin(θ) =((Sqrt[(x^2) + (y^2)])/r)
	sin(ϕ) =(y/(Sqrt[(x^2) + (y^2)]))
	cos(ϕ) =(x/(Sqrt[(x^2) + (y^2)]))
	
These are obtained by simple geometry (see any calculus book).  
For calculations I use also 

	r =(Sqrt[(x^2) + (y^2) + (z^2)])
	
but when I give the expression of the functions in cartesian coordinate I leave r in the formulae to make them shorter.

convCarts = { 
        r -> Sqrt[x^2 + y^2 + z^2]}
convCartp = {Cos[θ] -> z/r, 
        Sin[θ] -> Sqrt[x^2 + y^2]/r,  
        Cos[ϕ] -> x/Sqrt[x^2 + y^2],  
        Sin[ϕ] -> y/Sqrt[x^2 + y^2] }

pxAng[θ, ϕ] /. convCartp

Define the function

	px[r, θ, ϕ ][n] 
	
which gives (Subscript[np, x])(r,  θ,  ϕ) (these are 2  (Subscript[p, x])(r,  θ,  ϕ),  
      3  (Subscript[p, x])(
          r,  θ,  ϕ),   etc. )

Unprotect[px];
Clear[px]
px[r_, θ_, ϕ_][n_] := 
      Module[{ }, 
int = (SphericalHarmonicY[
                    1, -1, θ, ϕ] - 
                  SphericalHarmonicY[1, 1, θ, ϕ])/Sqrt[2] // 
            ComplexExpand  ; 
R[n, 1][r] *
          int ]
Protect[px];

Check whether the function works

px[r, θ, ϕ][3]

Convert the trigonometric functions to cartesian coordinates

px[r, θ, ϕ][3] /. convCartp

Make a table with the (Subscript[np, x]) orbitals:

\(\(tbp = 
        Table[{ToString[
              StringForm[\*""\!\(``p\_x\)(r, θ,ϕ)"", 
                n]], \(px[r, θ, ϕ]\)[
              n], \(px[r, θ, ϕ]\)[n] /. convCartp}, {n, 2, 
            5}];\)\)

\(TableForm[tbp, 
      TableHeadings -> {None, {\ "", "spherical", "cartesian"}}]\)

These functions  are all proportional to x and decay exponentially with r.  These two elements controls their shape.

Make a contour plot of 2  (Subscript[p, x])(x, y, z = 0)

\(pl2C = 
      ContourPlot[
        Evaluate[\(\(\(\(px[r, θ, ϕ]\)[2] /. convCartp\) /. 
                convCarts\) /. z -> 0\) /. {a -> 0.5292, 
              Z -> 1}], {x, \(-7\), 7}, {y, \(-5\), 
          5}, \ 
PlotPoints -> 200, \ 
        Contours -> 20, PlotRange -> All, 
        FrameLabel -> {"x", "z"}, \ 
        PlotLabel -> \*""2(Subscript[p, x])(x,y,z=0)""\ , 
ContourSmoothing -> True, 
        DefaultColor -> RGBColor[1, 0, 0]]\)

The function is elongated along the OX axis. The graph is whiter than the background where the function is positive and darker where the function is negative. The darker or the whiter a region, the higher the value of the function at those points.

Look at the function 2  (Subscript[p, x])(x,  y = 0, z = 0)

pl2x =  
      Plot[Evaluate[ px[r, θ, ϕ][2] /. convCartp /. 
                convCarts /. {y -> 0, z -> 0} /. {a -> 
                0.5292, Z -> 1}], {x, -12, 12}, 
        PlotPoints -> 100,  
PlotRange -> 
          All ]

You see that 2  (Subscript[p, x])(x,  y = 0,  z = 0) is negative when x ≤0 and positive otherwise

Look at the probability 
	(r^2)     2 (( (Subscript[p, x])(x,  y,  z))^2) for z=0

pl2p = 
      ContourPlot[
        Evaluate[ r^2 px[r, θ, ϕ][2]^2 /. 
                  convCartp /. convCarts /. 
              z -> 0 /. {a -> 0.5292, Z -> 1}], {x, -6, 
          6}, {y, -4, 4}, PlotPoints -> 100, 
        Contours -> 10, 
PlotRange -> All, 
        AxesLabel -> {"x", "z"}, ContourSmoothing -> True, 
        DefaultColor -> RGBColor[1, 0, 0]  ]

Plot 3  (Subscript[p, x])(x,  y,  z = 0)

\(pl3C = 
      ContourPlot[
        Evaluate[\(\(\(\(px[r, θ, ϕ]\)[3] /. convCartp\) /. 
                convCarts\) /. z -> 0\) /. {a -> 0.5292, 
              Z -> 1}], {x, \(-15\), 15}, {y, \(-10\), 10}, 
        PlotPoints -> 100, 
        Contours -> 10, 
PlotRange -> All, 
        FrameLabel -> {"x", "z"}\ , 
        PlotLabel -> \*""3\!\(p\_x\)(x, y,z=0)"", \ 
        DefaultColor -> RGBColor[1, 0, 0]]\)

Plot
	(r^2) 3 (( (Subscript[p, x])(x,  y, z))^2) as a function of x, for y=z=0

pl3x = 
      Plot[Evaluate[  px[r, θ, ϕ][3]  /. 
                  convCartp /. convCarts /. {z -> 0, 
                y -> 0} /. {a -> 0.5292, 
              Z -> 1}], {x, -15, 15}]

Plot the probability 
	
	(r^2) 3 (( (Subscript[p, x])(x,  y, z))^2) 
	
as a function of x, for  z=0

\(pl3p = 
      ContourPlot[
        Evaluate[\(\(\(r^2\ \(px[r, θ, ϕ]\)[3]^2 /. convCartp\) /. 
                convCarts\) /. z -> 0\) /. {a -> 0.5292, 
              Z -> 1}], {x, \(-15\), 15}, {y, \(-10\), 10}, 
        PlotPoints -> 100, 
        Contours -> 10, 
PlotRange -> All, 
        FrameLabel -> {"x", "z"}\ , 
        PlotLabel -> \*""3\!\(p\_x\)(x, y,z=0)"", \ 
        DefaultColor -> RGBColor[1, 0, 0]]\)

Plot
	(r^2) 3 (( (Subscript[p, x])(x,  y, z))^2) 
	
as a function of x, for y=z=0

Plot[
      Evaluate[r^2 px[r, θ, ϕ][3]^2 /. convCartp /. 
              convCarts /. {z -> 0, y -> 0} /. {a -> 
              0.5292, Z -> 1}], {x, -15, 15}]

Plot 
	4  (Subscript[p, x])(x,  y,  z = 0)

\(pl4C = 
      ContourPlot[
        Evaluate[\(\(\(\(px[r, θ, ϕ]\)[4] /. convCartp\) /. 
                convCarts\) /. z -> 0\) /. {a -> 0.5292, 
              Z -> 1}], {x, \(-20\), 20}, {y, \(-20\), 20}, 
        PlotPoints -> 100, 
        Contours -> 10, 
PlotRange -> All, 
        FrameLabel -> {"x", "z"}\ , 
        PlotLabel -> \*""4\!\(p\_x\)(x, y,z=0)"", \ 
        DefaultColor -> RGBColor[1, 0, 0]]\)

Plot
	3  (Subscript[p, x])(x,  y,  z) versus x, when y=z=0

pl4x = 
      Plot[Evaluate[  px[r, θ, ϕ][4]  /. 
                  convCartp /. convCarts /. {z -> 0, 
                y -> 0} /. {a -> 0.5292, 
              Z -> 1}], {x, -20, 20}]

Plot
	(r^2)  3 (( (Subscript[p, x])(x,  y,  z))^2) as a function of x and y, for z=0

\(pl4p = 
      ContourPlot[
        Evaluate[\(\(\(r^2\ \(px[r, θ, ϕ]\)[4]^2 /. convCartp\) /. 
                convCarts\) /. z -> 0\) /. {a -> 0.5292, 
              Z -> 1}], {x, \(-20\), 20}, {y, \(-20\), 20}, 
        PlotPoints -> 100, 
        Contours -> 10, 
PlotRange -> All, 
        FrameLabel -> {"x", "z"}\ , 
        PlotLabel -> \*""4\!\(p\_x\)(x, y,z=0)"", \ 
        DefaultColor -> RGBColor[1, 0, 0]]\)

Plot
	5  (Subscript[p, x])(x,  y,  z = 0)

\(pl5C = 
      ContourPlot[
        Evaluate[\(\(\(\(px[r, θ, ϕ]\)[5] /. convCartp\) /. 
                convCarts\) /. z -> 0\) /. {a -> 0.5292, 
              Z -> 1}], {x, \(-30\), 30}, {y, \(-20\), 20}, 
        PlotPoints -> 100, 
        Contours -> 10, 
PlotRange -> All, 
        FrameLabel -> {"x", "y"}\ , 
        PlotLabel -> \*""5\!\(p\_x\)(x, y,z=0)"", \ 
        DefaultColor -> RGBColor[1, 0, 0]]\)

Show 2  (Subscript[p, x]),  3  (Subscript[p, x]),  4 (Subscript[p, x]) and 5  (Subscript[p, x])on the same graph:

Show[GraphicsArray[{{pl2C, pl3C}, {pl4C, pl5C}}]]

Cell 5:  The functions 

	 (Subscript[np, y])(r,  θ, ϕ)=(i/(Sqrt[2])){(Subscript[Y, 1]^1)(θ, ϕ)+(Subscript[Y, 1]^(-1))(θ, ϕ)}(Subscript[R, n, 1])(r)

The angular part of (Subscript[np, y])is:

Clear[pyAng]
pyAng[θ_, ϕ_] := (I/
              Sqrt[2])*(SphericalHarmonicY[1, 1, θ, ϕ] + 
              SphericalHarmonicY[
                1, -1, θ, ϕ])//ComplexExpand

pyAng[θ, ϕ]

Unprotect[py];

    py =. 

    py[r_, θ_, ϕ_][n_] := 
      Module[{pyang}, 
pyang  = (I/
                    Sqrt[2])*(SphericalHarmonicY[1, 
                      1, θ, ϕ] + 
                    SphericalHarmonicY[1, -1, θ, ϕ]) // 
              ComplexExpand // 
            FullSimplify; 
pyang*R[n, 1][
            r]
]

py[r, θ, ϕ][n]

Convert to cartesian:

py[r, θ, ϕ][n] /. convCartp

Compare to (Subscript[np, x]):

px[r, θ, ϕ][n] /. convCartp

Note that you can obtain (Subscript[np, y]) by   interchanging x with y in (Subscript[np, x]). All plots on (Subscript[np, y])are the same as those of (Subscript[np, x])except that we exchange the labels x and y on the axes

Exercise:  Check the orthogonality of  2  (Subscript[p, x]) to 2  (Subscript[p, y])

Cell 6: The functions  

	 (Subscript[np, z])(θ, ϕ)= (Subscript[Y, 1]^0)(θ, ϕ) (Subscript[R, n, 1])(r)

Clear[pz]
pzAng[θ_, ϕ_] :=  
      SphericalHarmonicY[1, 0, θ, ϕ]

pzAng[θ, ϕ]

Unprotect[pz];
pz =. 
 pz[r_, θ_, ϕ_][n_] := 
      Module[{pz}, 
pz  =  
          SphericalHarmonicY[1, 
            0, θ, ϕ];  
pz*R[n, 1][
            r]]
Protect[pz];

Convert to cartesian coordinates:

pz[r,  θ,  ϕ][n] /. convCartp

Compare to (Subscript[np, x]):

px[r, θ, ϕ][n] /. convCartp

You can obtain (Subscript[np, z]) by interchanging x with z in (Subscript[np, x])

Cell 7: The function

	n(Subscript[d, xy]) (r,  θ, ϕ)= (Subscript[R, n, 2])(r)  \(TraditionalForm\`\(1\/\(i \@ 2\)\) {\(Y\_2\%2\)(θ, \ ϕ)\  - \(Y\_2\%\(-2\)\)(θ, \ ϕ)\ }\)

The angular part of (Subscript[nd, xy])

Clear[dxyAng]
dxyAng[θ_, ϕ_] := (1/( 
                I Sqrt[2])) *(SphericalHarmonicY[2, 
                2, θ, ϕ] -  
              SphericalHarmonicY[
                2, -2, θ, ϕ])//ComplexExpand

Test that it works:

dxyAng[θ, ϕ]

The conversion to cartesian coordinates is delicate, if I want to obtain the form usually given in textbooks, because of the presence of sin(2ϕ) for which we have not defined a conversion rule. However Mathematica can be used to write this expression in a convenient form. I'll do this  in several steps:

temp1 = (dxyAng[θ, ϕ] // TrigExpand) /. convCartp /. 
          convCarts // FullSimplify

I want to get rid of (x^2) + (y^2) + (z^2) and replace it with (r^2).  I use ReplacePart in combination with Position:

dxyCart=ReplacePart[temp1,r^2,Position[temp1, x^2 + y^2 + z^2] // Flatten]

Define the function dxy[r, θ, ϕ][n] that gives 

	(Subscript[nd, xy])(r,  θ, ϕ)

Clear[dxy]
dxy[r_, θ_, ϕ_][n_] := R[n, 2][r]*
        dxyAng[θ, ϕ]

Check that it works

dxy[r, θ, ϕ][3]

dxy[r, θ, ϕ][4]

Make a function that uses the coordinates x, y, z, and r (r is redundant, but I want to avoid using (Sqrt[(x^2) + (y^2) + (z^2)])(which is equal to r) so that the expressions are shorter.  When I use the function I convert r to its cartesian expression by using the rule convCarts).

dxyCartesian[n_] := dxyCart* R[n, 2][r]

dxyCartesian[3]

dxyCartesian[4]

The (Subscript[nd, xy]) functions are proportional to xy and this makes them large on the diagonals of the quadrants of the XOY plane.  They also decay exponentially with r. These two factors control the shape of the functions (Subscript[nd, xy]).

Plot 3  (Subscript[d, xy])  as a function of x and y for z=0.

\(pl3dxy = 
      ContourPlot[
        Evaluate[\ \(\(dxyCartesian[3]\  /. convCarts\) /. 
              z -> 0\) /. {a -> 0.5292, Z -> 1}], {x, \(-10\), 
          10}, {y, \(-10\), 10}, 
PlotPoints -> 50, \ 
        FrameLabel -> {"x", "y"}, 
        PlotLabel -> \*""3\!\(d\_xy\)""]\)

\(pl4dxy = 
      ContourPlot[
        Evaluate[\ \(\(dxyCartesian[4]\  /. convCarts\) /. 
              z -> 0\) /. {a -> 0.5292, Z -> 1}], {x, \(-20\), 
          20}, {y, \(-20\), 20}, 
PlotPoints -> 50, \ 
        FrameLabel -> {"x", "y"}, 
        PlotLabel -> \*""4\!\(d\_xy\)""]\)

\(pl5dxy = 
      ContourPlot[
        Evaluate[\ \(\(dxyCartesian[5]\  /. convCarts\) /. 
              z -> 0\) /. {a -> 0.5292, Z -> 1}], {x, \(-30\), 
          30}, {y, \(-30\), 30}, 
PlotPoints -> 50, \ 
        FrameLabel -> {"x", "y"}, 
        PlotLabel -> \*""5\!\(d\_xy\)""]\)

\(pl6dxy = 
      ContourPlot[
        Evaluate[\ \(\(dxyCartesian[6]\  /. convCarts\) /. 
              z -> 0\) /. {a -> 0.5292, Z -> 1}], {x, \(-40\), 
          40}, {y, \(-40\), 40}, 
PlotPoints -> 50, \ 
        FrameLabel -> {"x", "y"}, 
        PlotLabel -> \*""6\!\(d\_xy\)""]\)

Show[GraphicsArray[{{pl3dxy, pl4dxy}, {pl5dxy, pl6dxy }}]]

Cell 8: The functions 

	n(Subscript[d, yz])(r, θ, ϕ) =(Subscript[R, n,  ℓ])(r) \(TraditionalForm\`\(i\/\@2\) {\(Y\_2\%\(-1\)\)(θ, \ ϕ)\  + \ \(Y\_2\%1\)(θ, \ ϕ)}\)

Define the angular part of
	 n(Subscript[d, yz])(r, θ, ϕ)  
which is 
	\(TraditionalForm\`\(i\/\@2\) {\(Y\_2\%\(-1\)\)(θ, \ ϕ)\  + \ \(Y\_2\%1\)(θ, \ ϕ)}\)

Unprotect[dyzAng];
Clear[dyzAng]
dyzAng[θ_, ϕ_] := (I/ ( 
                  Sqrt[2]) )*(SphericalHarmonicY[
                  2, -1, θ, ϕ] +  
                SphericalHarmonicY[2, 1, θ, ϕ]) // 
          ComplexExpand // FullSimplify
Protect[dyzAng];

dyzAng[θ, ϕ]

Define the function dyz[r, θ, ϕ][n], which gives (Subscript[nd, yz])(r,  θ,  ϕ)

Unprotect[dyz];
Clear[dyz]
dyz[r_, θ_, ϕ_][n_] := 
      Module[{ }, 
dyzAng[θ, ϕ]*R[n, 2][
            r]]
Protect[dyz];

dyz[r, θ, ϕ][3]

Write the functions in cartesian coordinates since plotting is easier

dyz[r, θ, ϕ][n] /. convCartp

Compare this to (Subscript[d, xy])given in Cartesian coordinates.

dxyCartesian[n]

(Subscript[nd, yz]) can be obtained from (Subscript[nd, xy]) by exchanging x with z.  The plots of (Subscript[d, yz])in the YOZ plane will look like the plots of (Subscript[nd, xy]) in the XOY plane, with the OZ axis replacing the OX axis.

Cell 9: The functions 

	(Subscript[nd, xz]) (θ, ϕ)= (Subscript[R, n,  ℓ])(r) \(TraditionalForm\`\(1\/\@2\) {\(Y\_2\%\(-1\)\)(θ, \ ϕ)\  - \(Y\_2\%1\)(θ, \ ϕ)\ }\)

Unprotect[dxzAng];
Clear[dxzAng]
dxzAng[θ_, ϕ_] := (1/ 
                  Sqrt[2]) *(SphericalHarmonicY[
                    2, -1, θ, ϕ] -  
                  SphericalHarmonicY[2, 1, θ, ϕ]) // 
            ComplexExpand // FullSimplify;
Protect[dxzAng];

dxzAng[θ, ϕ]

In cartesian coordinates this is:

dxzAng[θ, ϕ] /. convCartp

Unprotect[dxz];
Clear[dxz];
dxz[r_, θ_, ϕ_][
        n_] := 
  dxzAng[θ,  ϕ]*R[n, 2][r] 
Protect[dxz];

dxz[r, θ, ϕ][n]

dxz[r, θ, ϕ][n] /. convCartp

This can be obtained from  (Subscript[d, xy]) by replacing y with z and keeping r unchanged. The plots of n(Subscript[d, xz])in the XOZ plane look like those of (Subscript[nd, xy])in the OXY plane; the OZ axis replaces the OY axis.

Cell 10:  The functions 

	n(Subscript[d, (x^2) - (y^2)]) (θ, ϕ)= (1/(Sqrt[2]))\(TraditionalForm\`{\(Y\_2\%2\)(θ, \ ϕ)\  + \ \(Y\_2\%\(-2\)\)(θ, \ ϕ)\ }\)

Define the function 
	dxSmySAng[θ, ϕ]
which gives the angular dependence of  
	 n (Subscript[d, (x^2) - (y^2)])

Unprotect[dxSmySAng];
Clear[dxSmySAng]
dxSmySAng[θ_, ϕ_] := (1/ 
              Sqrt[2]) *(SphericalHarmonicY[2, 
                2, θ, ϕ]
 + 
              SphericalHarmonicY[
                2, -2, θ, ϕ])//ComplexExpand 
Protect[dxSmySAng];

dxSmySAng[θ, ϕ]

Convert to cartesian coordinates:

temp1 = (dxSmySAng[θ, ϕ] /. convCartp // TrigExpand) /. 
        convCartp

temp2 = temp1 // FullSimplify

Give the result the name andx2my2

angdx2my2 = (x^2 - y^2) Sqrt[15/Pi]/(4 r^2)

Create a function that gives 
	(Subscript[nd, (x^2) - (y^2)])

Unprotect[dxSmyS];
Clear[dxSmyS]
dxSmyS[r_, θ_, ϕ_][n_] :=  
      dxSmySAng[θ,  ϕ]*R[n, 2][r]
Protect[dxSmyS];

dxSmyS[r, θ, ϕ][n]

dxSmyS[r, θ, ϕ][3]

dxSmyS[r, θ, ϕ][4]

Convert to cartesian coordinates x, y, z and r (as variables)

temp1 = ((dxSmyS[r, θ, ϕ][n] /. convCartp) // 
            TrigExpand) /. convCartp

temp2 = temp1 /. convCarts // FullSimplify

Replace (x^2) + (y^2) + (z^2)with (r^2):

temp3 = 
      ReplacePart[temp2, r^2, 
        Position[temp2,  x^2 + y^2 + z^2]]

Make a function out of this expression:

dxSmySCart[
        p_] :=  (1/r^2)*
          Sqrt[15/(2 Pi)] (x^2 - 
              y^2) (Z/(a n))^(3/2) Φ[n, 2][
            r/a ] /. n -> p

Look at two examples, for n=3 and n=4

dxSmySCart[3]

dxSmySCart[4]

Plot the function 3  (Subscript[d, (x^2) - (y^2)])in the OXY plane (z=0)

\(pl3dxSmyS = 
      ContourPlot[
        Evaluate[\ \(\(dxSmySCart[3] /. convCarts\) /. 
              z -> 0\) /. {a -> 0.5292, Z -> 1}], {x, \(-15\), 
          15}, {y, \(-15\), 15}, 
PlotPoints -> 50, \ 
        FrameLabel -> {"x", "y"}, 
        PlotLabel -> \*""3\!\(d\_\(x\^2 - y\^2\)\)""]\)

Plot the function 3  (Subscript[d, (x^2) - (y^2)])in the OXY plane (z=0)

\(pl4dxSmyS = 
      ContourPlot[
        Evaluate[\ \(\(dxSmySCart[4] /. convCarts\) /. 
              z -> 0\) /. {a -> 0.5292, Z -> 1}], {x, \(-25\), 
          25}, {y, \(-25\), 25}, 
PlotPoints -> 50, \ 
        FrameLabel -> {"x", "y"}, 
        PlotLabel -> \*""4\!\(d\_\(x\^2 - y\^2\)\)""]\)

Show[GraphicsArray[{pl3dxSmyS, pl4dxSmyS}]]

Cell 11:  The functions 

	n(Subscript[d, (z^2)]) (θ, ϕ) = (Subscript[R, n,  2])(
          r) (Subscript[Y, 2]^0)(θ,  ϕ)

Unprotect[dz2Ang];
Clear[dz2Ang]
dz2Ang[θ_, ϕ_] :=  
      SphericalHarmonicY[2, 0, θ, ϕ] 
Protect[dz2Ang];

dz2Ang[θ, ϕ]

dz2Ang[θ, ϕ] /. convCartp

Unprotect[dz2];
Clear[dz2 ]
 dz2 [r_, θ_, ϕ_][n_] :=  
 
      SphericalHarmonicY[2,  0, θ, ϕ] *R[n, 2][
          r]
Protect[dz2];

dz2 [r, θ, ϕ][3]

dz2[r, θ, ϕ][4]

Write the functions in cartesian coordinates since plotting is easy

dz2 [r, θ, ϕ][3] /. convCartp

dz2 [r, θ, ϕ][4] /. convCartp

\(pl3dz2 = 
      ContourPlot[
        Evaluate[\ \(\(\(\(dz2[r, θ, ϕ]\)[3] /. convCartp\) /. 
                convCarts\) /. z -> 0\) /. {a -> 0.5292, 
              Z -> 1}], {x, \(-25\), 25}, {y, \(-25\), 
          25}, 
PlotPoints -> 50, \ 
        FrameLabel -> {"x", "y"}, 
        PlotLabel -> \*""3\!\(d\_\(z\^2\)\)""]\)

\(pl3dz2 = 
      ContourPlot[
        Evaluate[\ \(\(\(\(dz2[r, θ, ϕ]\)[3] /. convCartp\) /. 
                convCarts\) /. y -> 0\) /. {a -> 0.5292, 
              Z -> 1}], {x, \(-25\), 25}, {z, \(-25\), 
          25}, 
PlotPoints -> 50, \ 
        FrameLabel -> {"x", "z"}, 
        PlotLabel -> \*""3\!\(d\_\(z\^2\)\)""]\)

\(pl4dz2 = 
      ContourPlot[
        Evaluate[\ \(\(\(\(dz2[r, θ, ϕ]\)[4] /. convCartp\) /. 
                convCarts\) /. y -> 0\) /. {a -> 0.5292, 
              Z -> 1}], {x, \(-25\), 25}, {z, \(-25\), 
          25}, 
PlotPoints -> 50, \ 
        FrameLabel -> {"x", "z"}, 
        PlotLabel -> \*""4\!\(d\_\(z\^2\)\)""]\)

Show[
      GraphicsArray[{{pl3dxSmyS, pl4dxSmyS}, {pl3dz2, pl4dz2}}]]

Cell 12: Collect the functional forms and make a table. This cell works only after you have run all other cells

\(\(tb = {{\*""(Subscript[np, x])"", 
            pxAng[θ, ϕ] R\_\(n, 1\)[r], 
            pxAng[θ, ϕ] R\_\(n, 1\)[r] /. 
              convCartp}, {\*""(Subscript[np, y])"", 
            pyAng[θ, ϕ] R\_\(n, 1\)[r], 
            pyAng[θ, ϕ] R\_\(n, 1\)[r] /. 
              convCartp}, {\*""(Subscript[np, z])"", 
            pzAng[θ, ϕ] R\_\(n, 1\)[r], 
            pzAng[θ, ϕ] R\_\(n, 1\)[r] /. 
              convCartp}, 
{\*""(Subscript[nd, xy])"", 
            dxyAng[θ, ϕ] R\_\(n, 2\)[r], 
            dxyCart*R\_\(n, 2\)[r]}, {\*""(Subscript[nd, yz])"", \ 
            dyzAng[θ, ϕ] 
              R\_\(n, 2\)[
                r], \((dyzAng[θ, ϕ] /. 
                  convCartp\ )\)\ R\_\(n, 2\)[
                r]}, \ {\*""(Subscript[nd, xz])"", 
            dxzAng[θ, ϕ] 
              R\_\(n, 2\)[
                r], \((dxzAng[θ, ϕ] /. 
                  convCartp\ )\)\ R\_\(n, 2\)[
                r]}, 
{\*""(Subscript[nd, (x^2) - (y^2)])"", dxSmySAng[θ, ϕ]\ \ R\_\(n, 2\)[r], 
            angdx2my2}, 
{\*""(Subscript[nd, (z^2)])"", 
            dz2Ang[θ, ϕ]\ R\_\(n, 2\)[
                r], \((dz2Ang[θ, ϕ] /. 
                  convCartp\ )\)\ R\_\(n, 2\)[r]}};\)\)

TableForm[tb]

WorkBookQM.18. 10. Spectroscopy: the transition moment

I have shown that if the light is linearly polarized, the only allowed transitions from the ground state are those in which

	{ℓ=0, m=0} →{ ℓ=1, m=0}
	
The transition dipole matrix element, which controls the intensity of these transitions is proportional to

	⟨n,  ℓ = 1 |  r  | 
              n = 1,  ℓ = 0⟩= (∫_0^∞) 
              (r^(2 )) (Subscript[R, n, 1])(
            r) r (Subscript[R, 1, 1])(r) dr		(1)
	
The analytical evaluation of these integrals is complicated and leads to expressions that are hard to evaluate. I prefer to use Mathematica to evaluate the individual integrals as needed. 
	The radial wave function is:

(Subscript[R, n, ℓ])(r)= (((((2 Z)/na)))^(3/2))\(\(Φ\_\(n, ℓ\)\)(ρ)\) 						(2)

with

(Subscript[Φ, n, ℓ])(r) = (Sqrt[(((n - ℓ - 1)!)/(2 n (n + ℓ)!))])\(\((\(2\ Z\ ρ\)\/n)\)\^\(\(\ \)\(ℓ\)\)\)  Exp[-((Z ρ)/n) ]\(\(L\_\(n - ℓ - 1\)\%\(2  ℓ + 1\)\)(\(2\ Z\ ρ\)\/n)\)		(3)

with 
	ρ=(( r)/( a))								(4)
and 
   	a=\(TraditionalForm\`\(\(ℏ\^2\)(4\ π\ ϵ\_0)\)\/\(μ\ e\^2\)\) 							(5)

Plan: 
1. Define the radial wave function. (Subscript[R, n,  ℓ])
2. Calculate the integral Eq. (1).

Needs["Utilities`CleanSlate`"]
CleanSlate[]
Unprotect[a,  R,  Z, Φ];
Remove["Global`*"]
Off[General::spell1]
Off[General::spell]
Protect[a,  R,  Z, Φ];

Cell 1:  The radial wave function.  This is copied from WorkBookQM.18.3

The function Φ:

Unprotect[Φ];
Clear[Φ];  Unprotect[a, Z]; Clear[a, Z];  
    Protect[a, Z];
\(\(Φ[n_Integer, ℓ_Integer]\)[ρ_] := 
      Module[{const}, 

        const = Sqrt[\(\(\ \)\(\((n - ℓ - 1)\)!\)\)\/\(\(\((n + ℓ)\)!\)\ 2\ n\)]; 
const\ \((2\ Z\ ρ/
                  n)\)^ℓ\ \ LaguerreL[n - ℓ - 1, 
            2\ ℓ + 
              1, \(2\ Z\ ρ\)\/n]\ Exp[\(-\(\(\(\ \)\(Z\ ρ\)\)\/n\)\)]]\)
Protect[Φ];

Define a function that generates (Subscript[R, n, ℓ])(r)

Unprotect[R];
Clear[R];
R[n_, ℓ_][
        r_] := (2  Z/ (n a))^(3/2) Φ[
            n, ℓ][r/a]
Protect[R];

Test that the function works

R[3, 2][r]

Cell 2. The transition dipole

To perform the integral defined by Eq. (1) I note that the integrand

(r^3) (Subscript[R, n, 1])(r) (Subscript[R, 1, 0])(
          r)

is a polynomial in r multiplied with an exponential that decays as r → ∞.   This means that I will not have trouble with the limit r →∞. Here is what I do to have control of the execution.  I denote by F(r) the definite integral 

	F(r) ≡∫(r^3) (Subscript[R, n, 1])(
            r) (Subscript[R, 1, 0])(r) dr
	
Then I recognize that

	 (∫_0^∞) 
              (r^3) (Subscript[R, n, 1])(r) (Subscript[R, 1, 0])(
            r) dr = F(∞)-F(0) = - F(0)
	 
In this way I do not give Mathematica a chance to worry about convergence and mess things up.  Here is a function transDip[n] that performs the integral

transDip[
        n_Integer] := -Integrate[r^3 R[n, 1] [r] R[1, 0][r], 
            r] /. r -> 0

transDip[2]

transDip[3]

Let us make a table with the results

toTable = Table[ transDip[n] // N, {n, 2, 10}]

Reduce the number of significant figures to four (this is too much but it costs us nothing)

toTb = NumberForm[#, 4] & /@ toTable

Make a table with the results:

TableForm[Transpose[{Table[n, {n, 2, 10}], toTb}], 
      TableHeadings -> {None, {"n", "⟨n,1|r|1,0\
⟩"}}]

We can calculate any transition dipole

	⟨(Subscript[n, f]), (Subscript[ℓ, f]) | r | 
          (Subscript[n, i,]) (Subscript[ℓ, i])⟩
	
but we have to make sure that  (Subscript[ℓ, f])≤(Subscript[n, f])-1 and (Subscript[ℓ, i])  ≤  (Subscript[n, i]) - 1 
The calculation is performed by the function shown below

transitionDipole[nf_Integer, ℓf_Integer][
        ni_Integer, ℓi_Integer] := -Integrate[
            r^3 R[nf, ℓf] [r] R[ni, ℓi][r], r] /. 
        r -> 0

For example
	⟨ 3, 2 | r |1, 0 ⟩ is

transitionDipole[3, 2][1, 0]