HW #11

1. Variational Method

Plot of the potential

The potential in the problem is V = 50 ((((e^(-x)) - 1))^2).

V = 50 ((((E^(-x)) - 1))^2);
Plot[V, {x, -1, 1}];

It is basically quadratic around zero and is very steep.  One may guess that the harmonic oscillator is a pretty good approximation.

Initial guess

Series[V, {x, 0, 3}]

If we try to identify the first term with the harmonic oscillator potential (1/2) m ω (x^2), we find ω = 10 because m = 1. Then one may hope that it would give the ground-state energy (1/2) ℏ ω = 5.  However, this hope is not quite fulfilled.  Using the groundstate wave function of the harmonic oscillator,

ψ1raw = ((((( m ω)/(π ℏ))))^(1/4)) 
          (E^(-m ω (x^2)/(2 ℏ)));
ψ1 = ψ1raw /. {m -> 1, ω -> 
            10, ℏ -> 1}

we compute the total energy using the Hamiltonian in the position space:

K1 = ((-(ℏ^2))/(2  m)) 
          Integrate[ψ1 D[ψ1, {x, 
                  2}], {x, -∞, ∞}] /. {ℏ -> 1, m -> 1}

P1 = 
      Integrate[ψ(1^2) 
          V, {x, -∞, ∞}]

Ebar1 = K1 + P1

N[Ebar1]

The error is

((N[Ebar1] - 39/8)/(39/8))

and is bigger than 5%, but is pretty good, accurate within 7.22%.

Therefore, we are motivated to try a non-SHO Gaussian as a trial function.

Gaussian trial function

ψ2 = (1/((π^(1/4)) (Δ^(1/2)))) 
          (E^(-(x^2)/(2  (Δ^2))));

K2 = ((-(ℏ^2))/(2  m)) 
          Integrate[ψ2 D[ψ2, {x, 
                  2}], {x, -∞, ∞}, 
            Assumptions -> Δ > 0] /. {ℏ -> 1, 
          m -> 1}

P2 = 
      Integrate[ψ(2^2) V, {x, -∞, ∞}, 
        Assumptions -> Δ > 0]

Ebar2 = K2 + P2

Ebar2min = 
      FindMinimum[Ebar2, {Δ, (10^(-1/2))}]

((Ebar2min[[1]] - 39/8)/(39/8))

It has improved, but not yet within 5%.

Linear times Gaussian

One way to improve it further is the following.  Note that the potential is not parity symmetric, and hence the groundstate wavefunction is not expected to be an even function.  Because the potential is lower on the right, we expect the wave function is skewed towards the right.  Therefore, we can try

ψ3 = (1 + k x) 
          (E^(-(x^2)/(2  (Δ^2))));

It is not normalized at this point, and we calculate its norm

norm2 = 
      Integrate[ψ(3^2), {x, -∞, ∞}, 
        Assumptions -> Δ > 0]

K3 = ((-(ℏ^2))/(2  m)) 
          Integrate[ψ3 D[ψ3, {x, 
                  2}], {x, -∞, ∞}, 
            Assumptions -> Δ > 0] /. {m -> 
            1, ℏ -> 1}

P3 = 
      Integrate[ψ(3^2) V, {x, -∞, ∞}, 
        Assumptions -> Δ > 0]

Ebar3 = ((K3 + P3)/norm2)

Ebar3min = 
      FindMinimum[
        Ebar3, {Δ, (10^(-1/2))}, {k, 0}]

((Ebar3min[[1]] - 39/8)/(39/8))

This is correct at the 1% level!  The wave function is therefore

ψf = ψ(3/(Sqrt[norm2])) /. Ebar3min[[2]]

Plot[ψf, {x, -1, 1}];

As expected, it is more or less a Gaussian, but skewed to the right.

Obviously, there are many ways to improve Gaussian.  I hope you found one successfully.

The analytic solution

This problem is actually a special case of the Morse potential

Vm = D ((((E^(-a u)) - 1))^2);

This potential is a form of the inter-atomic potential in diatomic molecules proposed by P.M. Morse, Phys. Rev. 34, 57 (1929).  The variable u is the distance between two atoms minus its equilibrium distance u = r - (Subscript[r, 0]).  The Schrödinger equation can be solved analytically.  Expanding it around the minimum,

Series[Vm, {u, 0, 4}]

Harmonic oscillator approximation gives ω = (Sqrt[((2  (a^2) D)/m)]).  The correct energy eigenvalues are known to be (Subscript[E, n]) = ℏ ω(
              n + (1/2)) - (((ℏ^2) (a^2))/(2  
                  m)) (((n + (1/2)))^2), where the second term is called the anharmonicanharmonic correction.  For large n, the second term will dominate and the energy appears to become negative.  Clearly, the bound state spectrum does not go forever, and ends at a certain value of n, quite different from the harmonic oscillator.  But this is expcted because the potential energy asymptotes to D for u -> +∞ and hence states for E > D must be unbound and have a continuous spectrum.

The ground-state wave function is known to have the form

ψm = (E^(-d (E^(-a u)))) 
          (E^(-b a u/2));

where b = 2  d - 1, d = ((Sqrt[2  D m])/(a ℏ)).  Let us see that it satisfies the Schrödinger equation.

Simplify[(((-((ℏ^2)/(2  m)) D[ψm, {u, 2}] + D ((((E^(-a u)) - 1))^2) ψm))/ψ)m /. {b -> 
              2  d - 1} /. {d -> ((Sqrt[2  D m])/(a ℏ))}]

The first term is \(TraditionalForm\`\(1\/2\) ℏ\ ω = ℏ\ a \@\( D\/\(2  m\)\)\) and the zero-point energy with the harmonic oscillator approximation.  The second term -(((ℏ^2) (a^2))/(8  m)) is the anharmonic correction.

The normalization is

norm2m = 
      Integrate[ψ(m^2), {u, -∞, ∞}, 
        Assumptions -> a > 0 && d > 0]

Actually, this wave function could have been guessed if one pays a careful attention to the asymptotic behaviors.  For u -> ∞, the potential asymptotes to a constant D, and hence the wave function must damp exponentially (e^(-κ u)) with κ = (Sqrt[2  m | E | ])/ℏ.  For u -> -∞, the potential rises extremely steeply as D (e^(-2 a u)).  It suggests that the energy eigenvalue becomes quickly irrelevant, and the behavior of the wave function must be given purely by the rising behavior of the potential.  By dropping the energy eigenvalue and looking at the Schrödinger equation,

	-((ℏ^2)/(2  
                    m)) (((d^2) ψ)/(d (u^2))) + 
          D (e^(-2 a u)) ψ = 0,

and change the variable to y = (e^(-a u)), we find

	\(TraditionalForm\`\(-\(ℏ\^2\/\(2  
                    m\)\)\) \((\(\(d\^2\) ψ\)\/\(d\ y\^2\) + \(1\/y\) \(d\ ψ\)\/\(d\ y\))\) + D\ ψ = 0\).

The second term in the parentheses is neglible for y -> ∞.  Therefore the wave function has the behavior ψ ∝ 
        (e^(-(Sqrt[2 m D]) y/ℏ a)).  Combining the behavior on both ends, the wave function has precisely the exact form given above.

This is the lesson: the one-dimensional potential problem is so simple that there are many ways to study the behavior of the wave function.  On the other hand, the real-world problem involves many more degrees of freedom.  In many cases, the Hamiltonian itself must be guessed.

Back to the Morse potential.  The case of the homework problem corresponds to D = 50, a = 1, m = 1, ℏ = 1.  Therefore the groundstate energy is \(TraditionalForm\`E\_0 = \(\(1\/2\) ℏ\ a \@\( D\/\(2  m\)\) - \(\(ℏ\^2\) a\^2\)\/\(8  m\) = \(5 - 1\/8 = 39\/8\)\)\).

ψmf = 
      Simplify[(1/(Sqrt[norm2m])) ψm /. {b -> 
                  2  d - 1} /. {d -> ((Sqrt[2  D m])/(a ℏ))} /. {D -> 50, a -> 1, m -> 1, ℏ -> 1}]

Plot[ψmf, {u, -1, 1}, PlotPoints -> 50];

The variational linear times Gaussian wave function was

Plot[ψf, {x, -1, 1}, PlotPoints -> 50];

Show[%, %%];

Quite close.

2. Fine-structure of the hydrogen atom

Preliminaries

We use the 2  s wavefunction

(Subscript[R, 2  s]) = (((2  a))^(-3/2)) (2 - (r/a)) 
          (E^(-r/(2  a)));

The full wavefunction is

(Subscript[ψ, 2  s]) = (1/(Sqrt[4  π])) (Subscript[R, 2  s]);

Here, a = (Subscript[a, 0])/Z = (ℏ^2)/(Z (e^2) m).  Similarly,

(Subscript[R, 2  p]) = (((2  a))^(-3/2)) (r/((Sqrt[3]) a)) 
          (E^(-r/(2  a)));

As a preparation for calculating the spin-orbit interaction, 

	\(TraditionalForm\`\(L\& -> \)·\(S\& -> \) = \(\(1\/2\) \((\(J\& -> \)\^2 - \(L\& -> \)\^2 - \(S\& -> \)\^2)\) = \(ℏ\^2\/2\) \((j(j + 1) - l(l + 1) - s(s + 1))\)\)\).
	
For j = l + (1/2),

	\(TraditionalForm\`\(L\& -> \)·\(S\& -> \) = \(\(ℏ\^2\/2\) \((\((l + 1\/2)\) \((l + 3\/2)\) - l(l + 1) - 
                3\/4)\) = \(ℏ\^2\/2\) l\)\),
	
while for j = l - (1/2),

	\(TraditionalForm\`\(L\& -> \)·\(S\& -> \) = \(\(ℏ\^2\/2\) \((\((l - 1\/2)\) \((l + 1\/2)\) - l(l + 1) - 
                3\/4)\) = \(-\(ℏ\^2\/2\)\) \((l + 1)\)\)\).
	
For the calculations of the relativistic corrections, we use the fact that

	\(TraditionalForm\`\(-\(1\/\(8 \( m\^3\) 
                  c\^2\)\)\) ⟨\((\(p\& -> \)\^2)\)\^2⟩ = \(-\(1\/\(8 \( m\^3\) 
                  c\^2\)\)\) ⟨\((\(-ℏ\^2\) Δ)\)\^2⟩\)
\(TraditionalForm\`\(\(=\)\(\(ℏ\^4\/\(8 \( m\^3\) 
                  c\^2\)\) \(∫\(d\^3\) 
                x\ \(ψ\^*\) ΔΔψ\) = \(-\(ℏ\^4\/\(8 \( m\^3\) 
                    c\^2\)\)\) \(∫\(d\^3\) \(\(x(Δψ)\)\^*\) \((Δψ)\)\)\)\)\)
\(TraditionalForm\`\(\(=\)\(\(-\(ℏ\^4\/\(8 \( m\^3\) 
                  c\^2\)\)\) \(∫\(d\^3\) \(x(
                d\^2\/\(d\ r\^2\) + \(2\/r\) 
                    d\/\(d\ r\) - \(l(l + 1)\)\/r\^2)\) \(\((R\ Y\_l\^m)\)\^*\) \((d\^2\/\(d\ r\^2\) + \(2\/r\) 
                    d\/\(d\ r\) - \(l(l + 1)\)\/r\^2)\) \((R\ Y\_l\^m)\)\)\)\)\)
=-((ℏ^4)/(8  (m^3) 
                  (c^2))) ∫(r^2) 
              d r (((d^2)/(d (r^2))) + (2/r) 
                    (d/(d r)) - ((l(l + 1))/(r^2))) R(
                r) (((d^2)/(d (r^2))) + (2/r) 
                    (d/(d r)) - ((l(l + 1))/(r^2))) R(
                  r).
	
The Darwin term for the Coulomb potential is proportional to 

	\(TraditionalForm\`Δ\ V\_c = \(Δ\ \(\(-Z\)\ e\^2\)\/r = \(4  π\ Z\ \(e\^2\) \(δ(\(x\& -> \))\) = 
            4  π\ \(ℏ\^2\/\(m\ a\)\) \(δ(\(x\& -> \))\)\)\)\).

2  s

First, the releativistic correction:

Integrate[-((ℏ^4)/(8  (m^3) 
                (c^2))) (((D[(Subscript[R, 2  s]), {r, 2}] + (2/r) 
                  D[(Subscript[R, 2  s]), r]))^2) (r^2), {r, 0, ∞}, 
      Assumptions -> a > 0]

Second, the spin-orbit interaction.  Because L = 0, it identically vanishes.

Third, the Darwin term.

((ℏ^2)/(8  (m^2) (c^2))) 
        4  π ((ℏ^2)/(m a)) ((Subscript[ψ, 2  s])^2) /. {r -> 
          0}

% + %%

2  (Subscript[p, 1/2])

First, the relativistic correction,

Integrate[-((ℏ^4)/(8  (m^3) 
                (c^2))) (((D[(Subscript[R, 2  p]), {r, 2}] + (2/r) 
                  D[(Subscript[R, 2  p]), r] - (2/(r^2)) (Subscript[R, 2  p])))^2) 
        (r^2), {r, 0, ∞}, Assumptions -> a > 0]

Second, the spin-orbit interaction.

Integrate[(g/(4  (m^2) 
              (c^2))) (1/(r^3)) ((ℏ^2)/(m a)) (-(ℏ^2)) ((Subscript[R, 2  p])^2) (r^2), {r, 0, ∞}, 
      Assumptions -> a > 0]

Third, the Darwin term.  Because the wavefunction vanishes at the origin, it is identically zero.

% + %% /. {g -> 2}

2  (Subscript[p, 3/2])

First, the relativistic correction.  We use the fact that \(TraditionalForm\`\(\(p\& -> \)\^2\) ψ = \(-ℏ\^2\) \(\(Y\_1\^m\)(d\^2\/\(d\ r\^2\) + \(2\/r\) d\/\(d\ r\) - 2\/r\^2)\) R\),

Integrate[-((ℏ^4)/(8  (m^3) 
                (c^2))) (((D[(Subscript[R, 2  p]), {r, 2}] + (2/r) 
                  D[(Subscript[R, 2  p]), r] - (2/(r^2)) (Subscript[R, 2  p])))^2) 
        (r^2), {r, 0, ∞}, Assumptions -> a > 0]

Second, the spin-orbit interaction.

Integrate[(g/(4  (m^2) 
              (c^2))) (1/(r^3)) ((ℏ^2)/(m a)) ((ℏ^2)/2) ((Subscript[R, 2  p])^2) (r^2), {r, 0, ∞}, 
      Assumptions -> a > 0]

Third, the Darwin term.  Because the wave function vanishes at the origin, it is identically zero.

% + %% /. {g -> 2}

Summary

Therefore, the 2  s and 2  (Subscript[p, 1/2]) states are still degenerate, and have the energy

	\(TraditionalForm\`\(-\(e\^2\/\(8  
                  a\)\)\) - \(5  ℏ\^4\)\/\(128 \( a\^4\) \(m\^3\) c\^2\) = \(-\(e\^2\/a\)\) \((1\/8 + \(5\/128\) α\^2)\)\),
	
while the 2  (Subscript[p, 3/2])states have the energy

	\(TraditionalForm\`\(-\(e\^2\/a\)\) \((1\/8 + \(1\/128\) α\^2)\)\),
	
where a = ((e^2)/(ℏ c)).  Numerically (we took g = 2 above for comparison between 2  s and 2  p), the energy shifts in eV are

constants = {e -> 1.60217653*(10^(-19))*2.997925*(10^9), 
          a -> 5.291772108*(10^(-9)), α -> 
            1/137.03599911};
ergtoev = 6.241506*(10^11);
(((-(e^2))/a) ((5/128) (α^2)) /. constants)*
      ergtoev
(((-(e^2))/a) ((1/128) (α^2)) /. constants)*
      ergtoev
% - %%

The fine structure is only 4.5×(10^(-5))eV in size!  Equivalently, 2  (Subscript[s, 1/2]) -> 2  (Subscript[p, 3/2])is a microwave transition of 2.76 cm.  This is why the physical constants and unit conversions used above are written to such high precision.  The values are taken from The NIST Reference on Constants, Units and Uncertainty (http://physics.nist.gov/cuu/index.html).

You should also look at the lecture note for the calculation of energy shifts for arbitrary n, l, j.  The states with the same j for a given n remain degenerate.  This degeneracy is lifted once the "Lamb shifts" are considered.

3. [optional] Harmonic oscillator

Exact result

The ground-state wavefunction for the unperturbed Hamiltonian is

(Subscript[ψ, 0]) = (((((m ω)/(π ℏ))))^(1/4)) 
          (E^(-m ω (x^2)/(2  ℏ)));

The exact result for the perturbed Hamiltonian is just

(Subscript[ψ, e]) = (Subscript[ψ, 0]) /. ω -> ω (Sqrt[ 1 + ϵ])

Simplify[Series[(Subscript[ψ, e]), {ϵ, 0, 2}]]

Perturbative result

We need to know the matrix elements. Using

	\(TraditionalForm\`V = \(\(1\/2\) ϵ\ m\ \(ω\^2\) 
            x\^2 = ϵ\ \(\(ℏ\ ω\)\/4\) \((a + a\^†)\)\^2\)\),
	
	\(TraditionalForm\`V\_00 = \(ϵ\ \(\(ℏ\ ω\)\/4\) ⟨0 | \((a + a\^†)\)\^2 | 
                0⟩ = ϵ\ \(ℏ\ ω\)\/4\)\),
	\(TraditionalForm\`V\_20 = \(ϵ\ \(\(ℏ\ ω\)\/4\) ⟨0 | \((a + a\^†)\)\^2 | 
                2⟩ = ϵ\ \(\(ℏ\ ω\)\/4\) \@2\)\),
	\(TraditionalForm\`V\_22 = \(ϵ\ \(\(ℏ\ ω\)\/4\) ⟨2 | \((a + a\^†)\)\^2 | 
                2⟩ = ϵ\ \(\(ℏ\ ω\)\/4\) \((3 + 2)\)\)\),
	\(TraditionalForm\`V\_42 = \(ϵ\ \(\(ℏ\ ω\)\/4\) ⟨2 | \((a + a\^†)\)\^2 | 
                4⟩ = ϵ\ \(\(ℏ\ ω\)\/4\) \@12\)\).
	
We use Eq. (5.1.44) in Sakurai:

	\(\(\(\(\(\(\(\(\(\(|\)\(0\)\)\(⟩\)\) = \(\(|\)\(0\^\((0)\)\)\)\)\(⟩\)\)\(+\)\) | 
                                        2\^\((0)\)\)\(⟩\)\) V\_20\/\(E\_0\^\((0)\) - E\_2\^\((0)\)\)\)\(+\)\)|(4^((0)))⟩(((Subscript[V, 42]) 
                                (Subscript[V, 20]))/((((Subscript[E, 0])^((0))) - 
                                    ((Subscript[E, 2])^((0)))) (((Subscript[E, 0])^((0))) - 
                                    ((Subscript[E, 4])^((0))))))
+|(2^((0)))⟩(((Subscript[V, 22]) 
                          (Subscript[V, 20]))/((((Subscript[E, 0])^((0))) - 
                              ((Subscript[E, 2])^((0)))) (((Subscript[E, 0])^((0))) - 
                              ((Subscript[E, 2])^((0))))))-|(2^((0)))⟩(((Subscript[V, 00]) 
                  (Subscript[V, 20]))/(((((Subscript[E, 0])^((0))) - ((Subscript[E, 2])^((0)))))^2))

	
	\(\(\(\(\(\(\(\(=\)\(\(|\)\(0\^\((0)\)\)\)\)\(⟩\)\)\(-\)\) | 
                                        2\^\((0)\)\)\(⟩\)\) ϵ\/\(4 \@ 2\)\)\(+\)\)|(4^((0)))⟩\(\(\(ϵ\^2\) \@6\)\/64\)+|(2^((0)))⟩(((ϵ^2) 
                        5)/(32 (Sqrt[ 2])))-|(2^((0)))⟩\(ϵ\^2\/\^\(32 \@ 2\)\)
	
	\(\(\(\(\(\(\(\(=\)\(\(|\)\(0\^\((0)\)\)\)\)\(⟩\)\)\(-\)\) | 
                                  2\^\((0)\)\)\(⟩\)\) ϵ\/\(4 \@ 2\)\)\(+\)\)|(4^((0)))⟩\(\(\(ϵ\^2\) \@6\)\/64\)+|(2^((0)))⟩((ϵ^2)/(8 (Sqrt[ 2]))).
	
This state is not normalized, and it needs to be renormalized by

	(Z^(-1))=⟨0 | 0⟩=1+((ϵ^2)/32)+O((ϵ^4)).
	
On the other hand, using (A.4.3) and (A.4.5) in Sakurai,

(Subscript[ψ, 2]) = (1/(Sqrt[2])) (((2  m ω (x^2))/ℏ) - 
              1) (((((m ω)/(π ℏ))))^(1/4)) 
          (E^(-m ω (x^2)/(2  ℏ)));

(Subscript[ψ, 4]) = (1/(2 (Sqrt[ 6]))) (4 (((((m ω (x^2))/ℏ)))^2) - 12 (( m ω (x^2))/ℏ) + 
              3) (((((m ω)/(π ℏ))))^(1/4)) 
          (E^(-m ω (x^2)/(2  ℏ)));

So the result from perturbation theory is

(Subscript[ψ, p]) = (1 - ((ϵ^2)/64)) ((Subscript[ψ, 0]) - (ϵ/(4 (Sqrt[ 2]))) (Subscript[ψ, 2]) + (((ϵ^2) (Sqrt[6]))/64) (Subscript[ψ, 4]) + ((ϵ^2)/(8 (Sqrt[ 2]))) (Subscript[ψ, 2]));

Simplify[Series[(Subscript[ψ, p]), {ϵ, 0, 2}]]

Its expansion agrees with that of the exact result.

4. Magnetic Field

We use the MKSA system.  The Hamiltonian of the electron in the magnetic field is written in the symmetric gauge where (Overscript[A,  ->])  = (B/2) (-y, x, 0),

	\(TraditionalForm\`H = \((\(p\& -> \) - e\ \(A\& -> \))\)\^2\/\(2  m\) - \(1\/\(4  π\ ϵ\_0\)\) \(Z\ e\^2\)\/r - 
          g\ \(e\/\(\(2\)  \(m\)\(\ \)\)\) \(s\& -> \)·\(B\& -> \)\)
\(TraditionalForm\`\(\(=\)\(\(1\/\(2  
                  m\)\) \((\((p\_x + e\ B\ y/2)\)\^2 + \((p\_y - e\ B\ x/2)\)\^2 + p\_z\^2)\) - \(1\/\(4  π\ ϵ\_0\)\) \(Z\ e\^2\)\/r - 
          g\ \(\(e\ B\)\/\(2  m\ c\)\) s\_z\)\)\)
\(TraditionalForm\`\(\(=\)\(\(1\/\(2  m\)\) \((\(p\& -> \)\^2 + 
                e\ \(B(\(p\_x\) y - \(p\_y\) x)\) + \(\(\(e\^2\) 
                        B\^2\)\/4\) \((x\^2 + 
                      y\^2)\))\) - \(1\/\(4  π\ ϵ\_0\)\) \(Z\ e\^2\)\/r - g\ \(\(e\ ℏ\ B\)\/\(2  m\ c\)\) m\_s\)\)\)
\(TraditionalForm\`\(\(=\)\(\(1\/\(2  m\)\) \((\(p\& -> \)\^2 - 
                e\ B\ L\_z + \(\(\(e\^2\) B\^2\)\/4\) \((x\^2 + 
                      y\^2)\))\) - \(1\/\(4  π\ ϵ\_0\)\) \(Z\ e\^2\)\/r - g\ \(\(e\ ℏ\ B\)\/\(2  m\ c\)\) \(\(m\_s\)\(.\)\)\)\)\)
	
The last term is already diagonalized for (Subscript[m, s]) = ±(1/2).

At O(B), the correction from the second term in the parenthesis vanishes because (Subscript[L, z]) = 0.  At O((B^2)), the last term in the parentheses contributes.

Because of the isotropy of the wave function, we can replace (x^2) + (y^2) -> (2/3) (r^2); computing:

(Subscript[R, 1  s]) = (a^(-3/2)) 2 
           (E^(-r/a));
Integrate[(((e^2) (B^2))/(8  
              m)) (2/3) (r^2) ((Subscript[R, 1  s])^2) (r^2), {r, 
        0, ∞}, Assumptions -> a > 0]

Another possible contribution at O((B^2)) coming from the second-order in the O(B) operator again vanishes because (Subscript[L, z]) = 0.  Therefore, the correction up to the second order,

	\(TraditionalForm\`Δ\ E = \(-g\)\ \(\(e\ ℏ\ B\)\/\(2\ m\)\) m\_s + \(a\^2\) \(\(e\^2\) B\^2\)\/\(4  m\)\).
	
It is useful to note that the second term leads to diamagnetism, a very weak repulsion from a magnet.  In addition, its effect on the index of refraction is quite small, suppressed by (α^2) compared to the effect of polarizability.  We see this by considering the modified permeability, where n is the number density of hydrogen atoms:

	\(TraditionalForm\`1\/μ\_0 + 
          n\ \(a\^2\) e\^2\/\(4  m\) = \(1\/μ\_0 + 
            n\ \(a\^2\) 
              a\ \(\(m\ e\^2\)\/\(4  π\ \(ϵ\_0\) ℏ\^2\)\) 
              e\^2\/\(4  m\) = 
          1\/μ\_0 + 
            n\ a\^3\ \(\(\ \)\(e\^4\)\)\/\(16  π\ \(ϵ\_0\) ℏ\^2\)\)\)
\(TraditionalForm\`\(\(=\)\(\(1\/μ\_0\) \((1 + 
                n\ a\^3\ \(\(\ \)\(\(e\^4\) μ\_0\)\)\/\(16  π\ \(ϵ\_0\) ℏ\^2\))\) = \(\(1\/μ\_0\) \((1 + 
                  n\ a\^3\ \(\(\ \)\(e\^4\)\)\/\(16  π\ \(ϵ\_0\^2\) \(ℏ\^2\) c\^2\))\) = \(1\/μ\_0\) \(\((1 + 
                    n\ a\^3\ π\ α\^2)\)\(\ \)\(.\)\)\)\)\)\)

5. Index of Refraction

Sakurai quotes the polarizability of the hydrogen atom α = (9/2) ((Subscript[a, 0])^3) (5.1.73), where (Subscript[a, 0]) is the Bohr radius.  Note that this result is in the Gaussian unit.  In the MKSA unit, the only modification is the energy levels from ((e^2)/(2 (Subscript[a, 0]))) (1/(n^2)) to \(TraditionalForm\`\(1\/\(4  π\ ϵ\_0\)\) \(e\^2\/\(2\ a\_0\)\) 1\/n\^2\), and hence the polarizability becomes α = 
        4  π (Subscript[ϵ, 0]) (9/2) ((Subscript[a, 0])^3).  The energy density of the electric field is ((Subscript[ϵ, 0])/2) (E^2), while the contribution of the hydrogen atoms is n (1/2) α (E^2), where n is the number density.  The total is then 

	(1/2) (E^2)((Subscript[ϵ, 0]) + 
              n α) = (1/2) (E^2) (Subscript[ϵ, 0])(
            1 + 4  π (9/2) ((Subscript[a, 0])^3) n).  
	
Therefore the permittivity of free space is changed to that of the gas, ϵ = (Subscript[ϵ, 0])(
          1 + 4  π (9/2) ((Subscript[a, 0])^3) n).  

At 0°C, 1atm, the number density of hydrogen molecules is \(TraditionalForm\`n = \(N\/V = \(p\/\(k\ T\) = \(101.3  kPa\)\/\(1.381\ \(10\^\(-23\)\) \(JK\^\(-1\)\) 273.15  K\)\)\)\), while (Subscript[a, 0]) = 0.529 (10^(-10))m.  We calculate, remembering to multiply the atomic polarizability above by 2:

n = 
        p/(k*T) /. {p -> 101300, k -> 1.381*(10^(-23)), 
            T -> 273.15};
eps = 
      4  π (9/2) ((Subscript[a, 0])^3)*2*n /. 
        (Subscript[a, 0]) -> 0.52910*(10^(-10))

Therefore, ϵ = (Subscript[ϵ, 0])* (1+0.000450).

The magnetic permeability also changes, but this correction is smaller by another factor of (α^2) (fine structure constant squared) using the calculation in the previous problem:

	\(TraditionalForm\`1\/μ\_0 + 
          n\ \(a\^2\) e\^2\/\(4  m\) = \(1\/μ\_0 + 
            n\ \(a\^2\) 
              a\ \(\(m\ e\^2\)\/\(4  π\ \(ϵ\_0\) ℏ\^2\)\) 
              e\^2\/\(4  m\) = 
          1\/μ\_0 + 
            n\ a\^3\ \(\(\ \)\(e\^4\)\)\/\(16  π\ \(ϵ\_0\) ℏ\^2\)\)\)
\(TraditionalForm\`\(\(=\)\(\(1\/μ\_0\) \((1 + 
                n\ a\^3\ \(\(\ \)\(\(e\^4\) μ\_0\)\)\/\(16  π\ \(ϵ\_0\) ℏ\^2\))\) = \(\(1\/μ\_0\) \((1 + 
                  n\ a\^3\ \(\(\ \)\(e\^4\)\)\/\(16  π\ \(ϵ\_0\^2\) \(ℏ\^2\) c\^2\))\) = \(1\/μ\_0\) \((1 + 
                  n\ a\^3\ π\ α\^2)\)\)\)\)\)

Therefore we ignore this correction.

The speed of light is given by (c^2)/(n^2) = 1/(ϵ μ), and hence in our case n = (Sqrt[ϵ/(Subscript[ϵ, 0])]):

(Sqrt[1 + eps])

The measurement shows n = 1.000140at λ = 590  mm.  The correction to the index of refraction is obtained with the correct order of magnitude, while not right on.  There are two reasons.  One is that we dealt with hydrogen atom, not molecule.  The second is that our calculation is for the spatially constant electric field, i.e. λ = ∞.

N.B.:  Sakurai's model for the polarization of the medium is even simpler than the usual treatment (cf. Jackson).  The polarization is not just \(TraditionalForm\`\(P\& -> \) = α\ n\ \(E\& -> \)\) where α is the polarizability from above and (Overscript[E,  ->])  is the externally applied field, but 

	\(TraditionalForm\`\(P\& -> \) = \(α\ n\ \(E\& -> \)\_total = α\ n\ \((\(E\& -> \) + \(\(4  π\)\/3\) \(P\& -> \))\)\)\).
	
In other words, the total electric field applied at a given location in the medium is equal to the sum of the externally applied field and the mean polarization of the medium.  Consequently, one finds

	\(TraditionalForm\`ϵ\/ϵ\_0 = 
        1 + \(n\ α\)\/\(1 - \(\(4  π\)\/3\) n\ α\) ≈ 
          1 + n\ α + \(\(4  π\)\/3\) \((n\ α)\)\^2 + 
            O\ [n\ α]\^3\).
	
In a denser medium where n α ~ 0.01 or greater this second order correction has a noticeable effect on the index of refraction.

(By the way, inverting the leftmost relation above gives the Clausius-Mossotti equation.)

6. [optional] Hyperfine splitting

See lecture notes on fine and hyperfine splitting.