Quantum Systems cont'd

The problem of calculating the energy levels of bound states in quantum mechanics, and not by refining the Hi-Lo trial-and-error approach, is far from trivial. This is exemplified in the next-to-simplest potential, the symetric potential well, treated in section 5.3 (p. 595–601). Using Mathematica to do the en route calculations (derivatives, evaluations, determinants) is invaluable.

The difficulty of the problem is appreciated from the fact that the equation that determines the energy levels is a determinant of a system of a matrix such as in Eq. (5.3.19)

detEq4E  = 
      Det[{
{Exp[-κ a], -Cos[k a], 
            Sin[k a], 
            0}, 
{κ Exp[-κ a], -k Sin[k a], -k Cos[k a], 0}, 
{0, -Cos[k a], 
            Sin[k a], Exp[-κ a]}, 
{0, 
            k Sin[k a], -k Cos[
                k a], -κ Exp[-κ a]}}]

FullSimplify[detEq4E]

the joint effect of which is also obtainable by post-processing:

detEq4E  = 
      Det[{
{Exp[-κ a], -Cos[k a], 
              Sin[k a], 
              0}, 
{κ Exp[-κ a], -k Sin[k a], -k Cos[k a], 0}, 
{0, -Cos[k a], 
              Sin[k a], Exp[-κ a]}, 
{0, 
              k Sin[k a], -k Cos[
                  k a], -κ Exp[-κ a]}}] // 
        FullSimplify

For every “interface”, where the potential jumps discontinuously, you have two matching conditions, \(TraditionalForm\`ψ(\(x\_i\^-\)) = ψ(\(x\_i\^+\))\) and \(TraditionalForm\`ψ' \((\(x\_i\^-\))\) = ψ' \((\(x\_i\^+\))\)\) — or, if the discontinuity in the potential is a δ-function, λ δ(x - (Subscript[x, 0])), the latter condition becomes \(TraditionalForm\`ψ' \((\(x\_0\^-\))\) - ψ' \((\(x\_i\^+\))\) ∝ \ λ\ \(ψ(x\_0)\)\). (Having enforced \(TraditionalForm\`ψ(\(x\_0\^-\)) = ψ(\(x\_0\^+\))\), this central value is well-defined.)

Since the solution in every region (between two adjacent “interfaces”) has two integration constants and there are two more (true) boundary conditions for bound states: ψ(x -> ±∞) = 
        0, for n interfaces, you have 2n+2 integration constants and the same number of boundary/matching conditions. Whence a determinant of an (2  n + 2) × (2  n + 2) matrix, and machine computation is definitely easy.

Now, this approach applies whenever the potential function crosses the energy level discontinuously.

That was not the case in the models we studied last time.

For potentials where the potential function crosses the energy level continuously—and no analytic (symbolic) solution can be obtained—you use numerical solutions (as we did last time), or WKB. (More about that later)

The section 5.5 in the text discusses AN anharmonic oscillator—the Pöschl-Teller potential,

-(Subscript[V, 0]) (sech^2)(x)

that, in fact and rather miraculously, admits an exact, analytic (symbolic) solution. (See text.)

3D

Motion in a central field is very well studied, since (1) many potentials in Nature do exhibit spherical symmetry, (2) at least in the first approximation, so that corrections can be dealt with subsequently.

Mathematica recognizes SphericalHarmonicY[], LegendreP[], Laguerre[], and more. The textbook also provides a package, “CentralField`” (p. 698–701) that defines more direct functions such as Radial[], Angle[], AnglePlot[], Orbital[].

Of course, even without this convenience package, you can plot things “directly”, using only the readily installed packages:

?SphericalPlot3D

is not a built-in function, but will be defined upon reading in the package:

<< Graphics`ParametricPlot3D`

?SphericalPlot3D

Soo...

Let's define H-atom wave-functions.

psi[n_, l_, m_, r_, θ_, ϕ_, Z_, a0_] := 
      Sqrt[((2 Z/a0)^3) (n - l - 
                    1)!/(2 n (n + 1)!)] 
        Exp[-(2 Z/(a0 n)) r/2] ((2 Z/(a0 n)) r)^
          l LaguerreL[n - l - 1, 2 l + 1, (2 Z/(a0 n)) r] 
        SphericalHarmonicY[l, m, θ, ϕ]

psiStar[n_, l_, m_, r_, θ_, ϕ_, Z_, a0_] := 
      Sqrt[((2 Z/a0)^3) (n - l - 
                    1)!/(2 n (n + 1)!)] 
        Exp[-(2 Z/(a0 n)) r/2] ((2 Z/(a0 n)) r)^
          l LaguerreL[n - l - 1, 2 l + 1, (2 Z/(a0 n)) r] 
        SphericalHarmonicY[l, -m, θ, ϕ]

{psi[3, 1, 1, r, Pi/3, Pi/5, 1, 1], 
      psiStar[3, 1, 1, r, Pi/3, Pi/5, 1, 1]}

SphericalPlot3D[
      Re[psi[5, 4, 2, 1, θ, ϕ, 1, 1]] , {θ, Pi/2, 
        Pi}, {ϕ, 0, Pi}]

Note: this takes considerable time (which was the point of the demonstration), so I cut down the tanges of the angles to do only a quarter of the job.