Calculation and conversion of  hydrogen atom wavefunctions and properties
Michael A. Morrison
Notebook prepared for Eric Abraham et al. 
May 31, 2000: Version 1.0

Introduction:
This notebook illustrates the use of some of the commands in the package HydrogenicAtom.m and associated Standard Mathematica packages which are useful in atomic physics. These packages will be used  for  test calculations in our study of perturbation theory energy shifts to hyperfine levels of Rb and Cs resulting from a detuned external monochromatic electromagnetic field. It also illustrates the simplest form of programming in Mathematica: constructing a Module

BookKeeping

ClearAll["Global`*"]

Load  Mathematica standard packages

Note: 
For illustrative purposes, I will load, then illustrate each package one at a time. Ordinarily, you should load all packages you will need later in the notebook in the BookKeeping section. If, as you work, you discover that you need another package, go ahead and load it where you are working, but when you revise the notebook move the Needs command to the BookKeeping section.

To calculate values of physical constants:

Needs["Miscellaneous`PhysicalConstants`"]

BohrRadius

ElectronMass

SpeedOfLight

To convert between different systems of units:

Needs["Miscellaneous`Units`"]

Convert[ElectronMass,  Gram]

Convert[BohrRadius, Fermi]

To obtain information about the chemical elements:

Needs["Miscellaneous`ChemicalElements`"]

AtomicNumber[Hydrogen]

Abbreviation[Helium]

AtomicWeight[Helium]

ElectronConfiguration[Rb]

ElectronConfigurationFormat[Rb]

StableIsotopes[Rb]

Load HydrogenicAtom.m, my package for generating hydrogen atom wave functions and properties

Note:
The package HydrogenicAtom.m must either be in the  directory that contains the notebook you're using or in a directory that is on your $Path. See Tip Sheets for ways to set up $Path most conveniently.

Needs["HydrogenicAtom`"]

?HydrogenicAtom

?HydrogenicAtom`*

?a0

?PhiFn

?Y

?Z

?uncrH

Generate radial functions, energies, and mean radii

Define a command that generates the reduced radial function for a given \(TraditionalForm\`\(\(\(|\)\(n\ ℓ\)\)\(⟩\)\)\) state in the position representation:

Clear[radialFunction]
radialFunction[n_, ℓ_] :=  uFnH[n, ℓ][
        r]

radialFunction[5, 4]

uFnH[5, 4][r]

(I checked these and many other cases against tables in the literature during package development.)

?atomicUnits

radialFunction[5, 4]  /.  atomicUnits

Calculate the energy of the (n  =  5) manifold:

?eH

Note that since the default value is what we want, we don't have to specify the second argument:

{eH[5], eH[5, 1],  eH[5, 2]}

To get a decimal number, we must feed the energy into the command N

eH[5] // N

Convert[(2 eH[5] // N) Rydberg, ElectronVolt]

Note:
To convert to eV we had to multiply the energy, which by default is in Hartree, by 2 to get it into Rydberg, since the Units` package doesn't recognize the unit Hartree.

n2EnergyeV  =  
      Convert[(2 eH[1] // N) Rydberg, ElectronVolt]

Generate one of the Hamiltonian eigenfunctions (Subscript[ψ, n])ℓ(r, θ, ϕ) for the (n, ℓ)  =  (5, 4) level:

psiH[5, 4, 3][r, θ, ϕ]

psiH[5, 4, 3][r, θ, ϕ]  /.  atomicUnits

Calculate the mean radius of the (n, ℓ)  =  (5, 4) level and convert the results to Angstroms, then to Meters:

meanrH[5, 4]  // N

Convert[(meanrH[5, 4]  // N) Bohr,  Angstrom]

Convert doesn't recognize the unit Bohr, so we must do the conversion to Angstrom "by hand"

BohrRadius

meanr54  =  (meanrH[5, 4] /. atomicUnits) * 
        BohrRadius

Convert[meanr54,  Angstrom]

Calculate the classical turning points for the 5g orbital

rctpH[5, 4]

There are two classical turning points, so the command rctpH returns a list of two values

rctpH[5, 4] // N

(rctpH[5, 4] /. atomicUnits // N)*BohrRadius

((rctpH[5, 4] /. atomicUnits // N)*
          BohrRadius)[[1]]

Convert only accepts one argument at a time, so we must extract each value from the returned list and feed it into Convert.

Convert[((rctpH[5, 4] /. atomicUnits // N)*
            BohrRadius)[[1]],  Angstrom]

Convert[((rctpH[5, 4] /. atomicUnits // N)*
            BohrRadius)[[2]],  Angstrom]

Showing off: Evaluate the mean radius of the 1100 p Rydberg state of atomic hydrogen:

meanrH[1100, 1]

meanrH[1100, 1]  /. a0 -> 1*BohrRadius  // 
      ScientificForm[#] &

Write a command to calculate the uncertainty in the radial position in a particular (n, ℓ) state

Description of the calculation:
The radial position uncertainty is calculated from the wave function \(TraditionalForm\`\(ψ\_\(n ℓ m\_ℓ\)\)(r, θ, ϕ)\) as
				\(TraditionalForm\`Δr\  = \ \@\(⟨r\^2⟩\  - \ ⟨r⟩\^2\)\)
where ⟨r⟩ is the mean radial position, which the command meanrH generates. I want to write this command so that it will work not only for hydrogenic systems but also for other analytic wave functions. See accopanying notes of 6/1/00 for details

?radialFunction

Note: I must use the atomicUnits replacement rules to get rid of the constants which HydrogenicAtom.m gives

Clear[uncertaintyIntegral]
\(uncertaintyIntegral[n_, ℓ_, p_]\  := 
	
      Integrate[
        radialFunction[n, ℓ]\^2*\ r\^p /. atomicUnits, {r, 
          0, ∞}]\)

Check the case where we know the right answer:

uncertaintyIntegral[1, 0, 1]

?meanrH

meanrH[1, 0]

Make sure that the p = 2 case executes:

uncertaintyIntegral[1, 0, 2]

Now define the command as described on the attached notes:

Clear[radialUncertainty]
radialUncertainty[
        n_, ℓ_]  :=  
   Module[{meanRadius, 
          meanRadiusSquared, uncertainty}, 
 
        meanRadius  =  
          uncertaintyIntegral[1, 0, 
            1]; 
meanRadiusSquared  =  
          uncertaintyIntegral[1, 0, 
            2]; 
uncertainty  = (Sqrt[meanRadiusSquared - (meanRadius^2) ])
]

Important Note:
Local variables defined in a Module are known by name to Mathematica but have no value associated with them.

?meanRadiusSquared

meanRadiusSquared

This contrasts with a global variable defined during a session, such as n2EnergyeV, which we defined above

?n2EnergyeV

Now we test our command on a case where we know or can find the answer: the (n, ℓ) = (1, 0) state of hydrogen:

radialUncertainty[1, 0]

Check my work using the command uncrH in HydrogenicAtom.m which evaluates (Subscript[((Δr)), n ℓ]) analytically for hydrogenic atoms:

uncrH[1, 0]

Aside: practical programming in Mathematica

I wrote the above command to illustrate the use of Module. It's actually more complicated than is necessary for this simple calculation. In practice, I would probably define the intermediate command uncertaintyIntegral as above, but I would just then define the function radialUncertainty directly, as:

Clear[radialUncertainty]
radialUncertainty[
        n_Integer?Positive, ℓ_Integer?
          NonNegative]  :=  
(Sqrt[uncertaintyIntegral[n, ℓ, 2]  -  uncertaintyIntegral[n, ℓ, 1]^2])

Note in the above that I have put in tests to "idiot-proof" the function commands:

{radialUncertainty[1, 0], radialUncertainty[1, -1],  
      radialUncertainty[0, 1]}