1. n (p^2) states

The purpose of this problem is to understand quantiatively the impact of the Coulomb repulsion among electrons in multi-electron atoms and the Hund's rule that picked a particular \(TraditionalForm\`\(\^\(2  S + 1\)\) L\_J\) state as the ground state. We are not concerned with the spin-orbit interactions here. The n (p^2) configuration allows us to deal with the Coulomb repulsion as if we are dealing with only two electrons (namely we assume that the closed shell is very rigid), and we can work out the energy splittings among \(TraditionalForm\`\(\^1\) S\), \(TraditionalForm\`\(\^1\) D\), \(TraditionalForm\`\(\^3\) P\) states explicitly.  Even though the radial integrals depend explicitly on the radial wave function and hence on details that rely on something like Hartree-Fock, the ratio of the energy splittings is determined purely based on Clebsch-Gordan coefficients thanks to the power of spherical symmetry of the problem.  At the same time, the comparison to data provides useful cross checks and shows the limitation of the anlaylsis.

(a) Calculation of the Coulomb repulsion

We evaluate the Coulomb repulsion for three angular momentum states.

The states

The two-body wave function of electrons in the n (p^2) configuration is given by either 
\(TraditionalForm\`\(\(ψ\_orbital\%symmetric\)(\(x\& -> \)\_1, \(x\& -> \)\_2)\) \(\(ψ\_spin\%\(anti - symmetric\)\)(
          m\_\(s\_1\), m\_\(s\_2\))\)\)
or
\(TraditionalForm\`\(\(ψ\_orbital\%\(anti - 
                symmetric\)\)(\(x\& -> \)\_1, \(x\& -> \)\_2)\) \(\(ψ\_spin\%symmetric\)(m\_\(s\_1\), m\_\(s\_2\))\)\)
The spin part is simple because the only anti-symmetric wave function is
\(TraditionalForm\`\(\(\(\(\(\(1\/\@2\) \((\(\(|\)\(1\/2\)\), \(-1\)\/2⟩\)\)\(-\)\) | \(-1\)\/2, 
          1\/2\)\(⟩\)\)\()\)\)\)
for S = 0, while there are three possible symmetric wave functions
\(TraditionalForm\`\(\(\(\(\(\(\(\(\(|\)\(1\/2\)\), 
              1\/2\)\(⟩\)\), \ \(\(\(1\/\@2\) \((\(\(|\)\(1\/2\)\), \(-1\)\/2⟩\)\)\(+\)\) | \(-1\)\/2, 
            1\/2\)\(⟩\)\)\()\)\), \ \(\(|\)\(\(-1\)\/2\)\), \(-1\)\/2\)\(⟩\)\)\)
for S = 1.  For the evaluation of the Coulomb repulsion between two electrons, the spins play no role except for dictating the symmetry or anti-symmetry of the orbital wave function.

Because both electrons occupy the same orbital, the radial wave function is the product of the same function for two electron positions.  The only non-trivial part is the spherical harmonics so that a specific total L is chosen.  The form is
\(TraditionalForm\`\(ψ\_orbital\)(\(x\& -> \)\_1, \(x\& -> \)\_2) = \(R(r\_1)\) \(R(
            r\_2)\) \(∑\_\(m\_1, m\_2\)\(\(Y\_l\%\(m\_1\)\)(Ω\_1)\) \(\(Y\_l\%\(m\_2\)\)(Ω\_2)\) ⟨l\ l\ m\_1\ m\_2 | L\ M\_L⟩\)\)
The total orbital angular momentum is then specified, L and (Subscript[M, L]). The expression is summed over (Subscript[m, 1]) and (Subscript[m, 2]), but it is actually only one sum because (Subscript[M, L]) = (Subscript[m, 1]) + (Subscript[m, 2]).  But it doesn't hurt to keep both of them because wrong combinations automatically give zero contributions.  The important point is that the orbital wave function is either symmetric or anti-symmetric automatically due to the symmetry property of the Clebsch-Gordan coefficients.  As I mentioned in class, when you add the same integer angular momenta, l and l in our case, an even total angular momentum (L even) gives a symmetric wave function, while an odd total angular momentum (L odd) gives an anti-symmetric one.  I don't find this statement in Sakurai, but can be found in many other places.  For example at MathWorld, http://mathworld.wolfram.com/Clebsch-GordanCoefficient.html, Eq. (10) makes it clear.  Therefore, once the wave function is expressed as a summation over (Subscript[m, 1]) and (Subscript[m, 2]) as above, the (anti-)symmetry is already taken care of.  This is a special simplification when we add the same l, and doesn't apply when two electrons are in different orbitals.

Explicitly, the \(TraditionalForm\`\(\(\^1\)\(D\)\)\) state |L = 2, (Subscript[M, L]) = 2, S = 0, 
        (Subscript[M, S]) = 0⟩ is very simple because there is only one term in the orbital wave function, 
\(TraditionalForm\`\(\(\(ψ(\(x\& -> \)\_1, \(x\& -> \)\_2) = \(\(\(R(r\_1)\) \(R(
                    r\_2)\) \(\(Y\_1\%1\)(Ω\_1)\) \(\(Y\_1\%1\)(Ω\_2)\) \(1\/\@2\) \((\(\(|\)\(1\/2\)\), \(-1\)\/2⟩\)\)\(-\)\) | \(-1\)\/2, 
          1\/2\)\(⟩\)\)\()\)\)\)
Note the symmetry in the orbital wave function.

The 1  S state is more complicated because it involves three terms
\(TraditionalForm\`\(\(\(ψ(\(x\& -> \)\_1, \(x\& -> \)\_2) = \(R(r\_1)\) \(R(
                r\_2)\) \(1\/\@3\) \((\(\(Y\_1\%1\)(Ω\_1)\) \(\(Y\_1\%\(-1\)\)(Ω\_2)\) - \(\(Y\_1\%0\)(Ω\_1)\) \(\(Y\_1\%0\)(Ω\_2)\) + \(\(Y\_1\%\(-1\)\)(Ω\_1)\) \(\(Y\_1\%1\)(Ω\_2)\))\)

          \(\(\(1\/\@2\) \((\(\(|\)\(1\/2\)\), \(-1\)\/2⟩\)\)\(-\)\) | \(-1\)\/2, 
          1\/2\)\(⟩\)\)\()\)\)\)

The \(TraditionalForm\`\(\(\^3\)\(P\)\)\) state |L = 1, (Subscript[M, L]) = 0, S = 1, 
        (Subscript[M, S]) = 1⟩ as an example is
\(TraditionalForm\`\(\(ψ(\(x\& -> \)\_1, \(x\& -> \)\_2) = \(R(r\_1)\) \(R(
                r\_2)\) \(1\/\@2\) \((\(\(Y\_1\%1\)(Ω\_1)\) \(\(Y\_1\%\(-1\)\)(Ω\_2)\) - \(\(Y\_1\%\(-1\)\)(Ω\_1)\) \(\(Y\_1\%1\)(Ω\_2)\))\) | 1\/2, 
        1\/2\)\(⟩\)\)\)

working it out

The Coulomb potential is expanded in the spherical harmonics in the usual way,
\(TraditionalForm\`e\^2\/r\_12 = \(e\^2\) \(∑\_\(λ, μ\)\(\(r\_ < \)\^λ\/\(r\_ > \)\^\(λ + 1\)\) \(\(4  π\)\/\(2  λ + 
                    1\)\) \(\(\(Y\_λ\%μ\)(Ω\_1)\)\^*\) \(\(Y\_λ\%μ\)(Ω\_2)\)\)\).
Here, I avoided to use l', m' because they are hard to distinguish, and instead used λ, μ.

Putting everything together, the expectation value is
\(TraditionalForm\`⟨n\ p\^2, 
          L\ M\_L | e\^2\/r\_12 | n\ p\^2, 
          L\ M\_L⟩ = ∫d\ \(\(x\& -> \)\_1\) \(∫d\ \(\(x\& -> \)\_2\) \(ψ\_spin\^*\) \(ψ\_orbital\^*\) \(e\^2\/r\_12\) \(ψ\_orbital\) ψ\_spin\)\)
The spin part of the wave function gives just its normalization, namely 1.   Now writing it out explicitly, 
\(TraditionalForm\`∫d\ \(\(x\& -> \)\_1\) \(∫d\ \(\(\(\(x\& -> \)\_2\)(\(R(r\_1)\) \(R(
                        r\_2)\) \(∑\_\(\(m'\)\_1, \(m'\)\_2\)\(\(Y\_l\%\(\(m'\)\_1\)\)(Ω\_1)\) \(\(Y\_l\%\(\(m'\)\_2\)\)(Ω\_2)\) ⟨l\ l\ \(m'\)\_1\ \(m'\)\_2 | 
                              L\ M\_L⟩\))\)\^*\) \((\(e\^2\) \(∑\_\(λ, μ\)\(\(r\_ < \)\^λ\/\(r\_ > \)\^\(λ + 1\)\) \(\(4  π\)\/\(2  λ + 
                            1\)\) \(\(\(Y\_λ\%μ\)(Ω\_1)\)\^*\) \(\(Y\_λ\%μ\)(Ω\_2)\)\))\) \((\(R(
                    r\_1)\) \(R(
                    r\_2)\) \(∑\_\(m\_1, m\_2\)\(\(Y\_l\%\(m\_1\)\)(Ω\_1)\) \(\(Y\_l\%\(m\_2\)\)(Ω\_2)\) ⟨l\ l\ m\_1\ m\_2 | L\ M\_L⟩\))\)\)\)
Note that in the "bra" or complex-conjugated wave function, the summation is taken over (Subscript[(m'), 1]) and (Subscript[(m'), 2])to be distinguished from the summation over (Subscript[m, 1])and (Subscript[m, 2]) in the "ket".  This is very important.

Now we switch to the spherical coordinates and use the identity Eq. (3.7.73) in Sakurai,
\(∫d\ Ω\ \(Y\_l\%m\^*\) \(Y\_\(l\_1\)\%\(m\_1\)\) Y\_\(l\_2\)\%\(m\_2\)\)=(Sqrt[(((2  (Subscript[l, 1]) + 1) (2  (Subscript[l, 2]) + 1))/(4  π(
                        2  l + 1)))])⟨(Subscript[l, 1]) (Subscript[l, 2]) 00|l 0⟩\(⟨\(l\_1\) \(l\_2\) \(m\_1\) m\_2 | l\ m⟩\).
We also use the fact that the Clebsch-Gordan coefficients are all real.  We obtain
⟨n (p^2), L (Subscript[M, L]) | ((e^2)/(Subscript[r, 12])) | n (p^2), 
          L (Subscript[M, L])⟩
\(\(=\)\(∑\_\(λ, μ, m\_1, m\_2, \(m'\)\_1, \(m'\)\_2\)∫\(r\_1\^2\) 
                d\ \(r\_1\) \(∫\(r\_2\^2\) 
                    d\ \(r\_2\) \(\(\(e\^2\) \(r\_ < \)\^λ\)\/\(r\_ > \)\^\(λ + 1\)\) \(\(4  π\)\/\(2  λ + 1\)\) \(\(R(
                          r\_1)\)\^2\) \(\(R(
                          r\_2)\)\^2\) ⟨l\ l\ \(m'\)\_1\ \(m'\)\_2 | L\ M\_L⟩ ⟨l\ l\ m\_1\ m\_2 | 
                        L\ M\_L⟩\)\)\)
(Sqrt[(((2  λ + 1))/(4  π))])⟨l λ 0 0|l 0⟩\(⟨l\ λ\ \(\(m'\)\_1\) μ | 
              l\ m\_1⟩\)(Sqrt[(((2  λ + 
                      1))/(4  π))])⟨l λ 0 0|l 0⟩\(⟨l\ λ\ \(\(m'\)\_2\) μ | l\ m\_2⟩\)
The radial integral depends only on λ, and we call it
\(TraditionalForm\`\(F\^λ\)(n\ p, 
          n\ p) = ∫\(r\_1\^2\) 
            d\ \(r\_1\) \(∫\(r\_2\^2\) 
                d\ \(r\_2\) \(\(\(e\^2\) \(r\_ < \)\^λ\)\/\(r\_ > \)\^\(λ + 1\)\) \(\(R(r\_1)\)\^2\) \(R(r\_2)\)\^2\)\)
Then the result is
⟨n (p^2), L (Subscript[M, L]) | ((e^2)/(Subscript[r, 12])) | n (p^2), 
          L (Subscript[M, L])⟩
=(∑_λ)(F^λ)(n p, n p)(∑_(μ, (Subscript[m, 1]), (Subscript[m, 2]), (Subscript[(m'), 1]), (Subscript[(m'), 2])))⟨l l (Subscript[(m'), 1]) (Subscript[(m'), 2]) | 
                      L (Subscript[M, L])⟩\(⟨l\ l\ m\_1\ m\_2 | L\ M\_L⟩\)((⟨l λ 0 0|l 0⟩)^2)\(⟨l\ λ\ \(\(m'\)\_1\) μ | l\ m\_1⟩\)\(⟨l\ λ\ \(m\_2\) μ | l\ \(m'\)\_2⟩\)

The summation is over μ, (Subscript[m, 1]), (Subscript[m, 2]), (Subscript[(m'), 1]), (Subscript[(m'), 2]), but there are relations (Subscript[M, L]) = (Subscript[m, 1]) + (Subscript[m, 2]) = (Subscript[(m'), 1]) + (Subscript[(m'), 2]), μ = (Subscript[m, 1]) - (Subscript[(m'), 1]) = (Subscript[(m'), 2]) - (Subscript[m, 2]) and there is actually only one summation, while there is an overall summation over λ.  But it doesn't hurt to sum over all of them because the wrong combinations automatically give zero.

The general expression would be

Sum[
      ClebschGordan[{l, (Subscript[(m'), 1])}, {l, (Subscript[(m'), 2])}, {L, (Subscript[M, L])}] 
        ClebschGordan[{l, (Subscript[m, 1])}, {l, (Subscript[m, 2])}, {L, 
            (Subscript[M, L])}] ClebschGordan[{l, 0}, {λ, 0}, {l, 0}]^2 
        ClebschGordan[{l, (Subscript[(m'), 1])}, {λ, μ}, {l, (Subscript[m, 1])}] 
        ClebschGordan[{l, 
            (Subscript[m, 2])}, {λ, μ}, {l, (Subscript[(m'), 2])}], {(Subscript[m, 1]), -l, 
        l}, {(Subscript[m, 2]), -l, l}, {(Subscript[(m'), 1]), -l, l}, {(Subscript[(m'), 2]), -l, 
        l}, {μ, -λ, λ}]

If you substitute (Subscript[M, L]) = 0 etc into the above expression, Mathematica complains.  You have to write "0" in there explicitly.  I also take the sum over all of the variables, and most of them are unphysical because they don't satisfy the requirements like (Subscript[m, 1]) + (Subscript[m, 2]) = (Subscript[m, 3]).  Mathematica returns zero in such cases and there is no harm to sum over all of them, together with generous warning signs you are forced to stare at.

apply it to our problem

From the spherical symmetry of the problem, the result should be independent of the choice of (Subscript[M, L]).  We take (Subscript[M, L]) = 0 for the purpose of actual evaluations.  Because we have l = 1, it is

Sum[
      ClebschGordan[{1, (Subscript[(m'), 1])}, {1, (Subscript[(m'), 2])}, {L, 0}] 
        ClebschGordan[{1, (Subscript[m, 1])}, {1, (Subscript[m, 2])}, {L, 
            0}] ClebschGordan[{1, 0}, {λ, 0}, {1, 0}]^2 
        ClebschGordan[{1, (Subscript[(m'), 1])}, {λ, μ}, {1, (Subscript[m, 1])}] 
        ClebschGordan[{1, 
            (Subscript[m, 2])}, {λ, μ}, {1, (Subscript[(m'), 2])}], {(Subscript[m, 1]), -1, 
        1}, {(Subscript[m, 2]), -1, 1}, {(Subscript[(m'), 1]), -1, 1}, {(Subscript[(m'), 2]), -1, 
        1}, {μ, -λ, λ}]

\(TraditionalForm\`\(\(\^1\)\(S\)\)\)

The coefficient of (F^0) is

Sum[
      ClebschGordan[{1, (Subscript[(m'), 1])}, {1, (Subscript[(m'), 2])}, {0, 0}] 
        ClebschGordan[{1, (Subscript[m, 1])}, {1, (Subscript[m, 2])}, {0, 
            0}] ClebschGordan[{1, 0}, {0, 0}, {1, 0}]^2 
        ClebschGordan[{1, (Subscript[(m'), 1])}, {0, μ}, {1, (Subscript[m, 1])}] 
        ClebschGordan[{1, (Subscript[m, 2])}, {0, μ}, {1, (Subscript[(m'), 2])}], {(Subscript[m, 1]), -1, 
        1}, {(Subscript[m, 2]), -1, 1}, {(Subscript[(m'), 1]), -1, 1}, {(Subscript[(m'), 2]), -1, 
        1}, {μ, 0, 0}]

and the coefficient of (F^2) is

Sum[
      ClebschGordan[{1, (Subscript[(m'), 1])}, {1, (Subscript[(m'), 2])}, {0, 0}] 
        ClebschGordan[{1, (Subscript[m, 1])}, {1, (Subscript[m, 2])}, {0, 
            0}] ClebschGordan[{1, 0}, {2, 0}, {1, 0}]^2 
        ClebschGordan[{1, (Subscript[(m'), 1])}, {2, μ}, {1, (Subscript[m, 1])}] 
        ClebschGordan[{1, (Subscript[m, 2])}, {2, μ}, {1, (Subscript[(m'), 2])}], {(Subscript[m, 1]), -1, 
        1}, {(Subscript[m, 2]), -1, 1}, {(Subscript[(m'), 1]), -1, 1}, {(Subscript[(m'), 2]), -1, 
        1}, {μ, -2, 2}]

\(TraditionalForm\`\(\(\^1\)\(D\)\)\)

The coefficient of (F^0)is

Sum[
      ClebschGordan[{1, (Subscript[(m'), 1])}, {1, (Subscript[(m'), 2])}, {2, 0}] 
        ClebschGordan[{1, (Subscript[m, 1])}, {1, (Subscript[m, 2])}, {2, 
            0}] ClebschGordan[{1, 0}, {0, 0}, {1, 0}]^2 
        ClebschGordan[{1, (Subscript[(m'), 1])}, {0, μ}, {1, (Subscript[m, 1])}] 
        ClebschGordan[{1, (Subscript[m, 2])}, {0, μ}, {1, (Subscript[(m'), 2])}], {(Subscript[m, 1]), -1, 
        1}, {(Subscript[m, 2]), -1, 1}, {(Subscript[(m'), 1]), -1, 1}, {(Subscript[(m'), 2]), -1, 
        1}, {μ, -0, 0}]

and the coefficient of (F^2)is

Sum[
      ClebschGordan[{1, (Subscript[(m'), 1])}, {1, (Subscript[(m'), 2])}, {2, 0}] 
        ClebschGordan[{1, (Subscript[m, 1])}, {1, (Subscript[m, 2])}, {2, 
            0}] ClebschGordan[{1, 0}, {2, 0}, {1, 0}]^2 
        ClebschGordan[{1, (Subscript[(m'), 1])}, {2, μ}, {1, (Subscript[m, 1])}] 
        ClebschGordan[{1, (Subscript[m, 2])}, {2, μ}, {1, (Subscript[(m'), 2])}], {(Subscript[m, 1]), -1, 
        1}, {(Subscript[m, 2]), -1, 1}, {(Subscript[(m'), 1]), -1, 1}, {(Subscript[(m'), 2]), -1, 
        1}, {μ, -2, 2}]

\(TraditionalForm\`\(\(\^3\)\(P\)\)\)

The coefficient of (F^0)is

Sum[
      ClebschGordan[{1, (Subscript[(m'), 1])}, {1, (Subscript[(m'), 2])}, {1, 0}] 
        ClebschGordan[{1, (Subscript[m, 1])}, {1, (Subscript[m, 2])}, {1, 
            0}] ClebschGordan[{1, 0}, {0, 0}, {1, 0}]^2 
        ClebschGordan[{1, (Subscript[(m'), 1])}, {0, μ}, {1, (Subscript[m, 1])}] 
        ClebschGordan[{1, (Subscript[m, 2])}, {0, μ}, {1, (Subscript[(m'), 2])}], {(Subscript[m, 1]), -1, 
        1}, {(Subscript[m, 2]), -1, 1}, {(Subscript[(m'), 1]), -1, 1}, {(Subscript[(m'), 2]), -1, 
        1}, {μ, -0, 0}]

and the coefficient of (F^2)is

Sum[
      ClebschGordan[{1, (Subscript[(m'), 1])}, {1, (Subscript[(m'), 2])}, {1, 0}] 
        ClebschGordan[{1, (Subscript[m, 1])}, {1, (Subscript[m, 2])}, {1, 
            0}] ClebschGordan[{1, 0}, {2, 0}, {1, 0}]^2 
        ClebschGordan[{1, (Subscript[(m'), 1])}, {2, μ}, {1, (Subscript[m, 1])}] 
        ClebschGordan[{1, (Subscript[m, 2])}, {2, μ}, {1, (Subscript[(m'), 2])}], {(Subscript[m, 1]), -1, 
        1}, {(Subscript[m, 2]), -1, 1}, {(Subscript[(m'), 1]), -1, 1}, {(Subscript[(m'), 2]), -1, 
        1}, {μ, -2, 2}]

Putting together

Putting them together, we find
\(TraditionalForm\`⟨\(n\ p\^2\) \(\(\(,\)\(\ \)\)\^1\) S | 
            e\^2\/r\_12 | \(n\ p\^2\) \(\(\(,\)\(\ \)\)\^1\) 
              S⟩ = F\^0 + \(2\/5\) F\^2\)
\(TraditionalForm\`⟨\(n\ p\^2\) \(\(\(,\)\(\ \)\)\^1\) D | 
            e\^2\/r\_12 | \(n\ p\^2\) \(\(\(,\)\(\ \)\)\^1\) 
              D⟩ = F\^0 + \(1\/25\) F\^2\)
\(TraditionalForm\`⟨\(n\ p\^2\) \(\(\(,\)\(\ \)\)\^3\) P | 
            e\^2\/r\_12 | \(n\ p\^2\) \(\(\(,\)\(\ \)\)\^3\) 
              P⟩ = F\^0 - \(1\/5\) F\^2\)

side note

Just in case you are dismayed by the complicated combination of Clebsch-Gordan coefficients above, it can be recast to a more compact form using Wigner's 3j- and 6j-symbols.  Messiah has a good compilation of useful formulae.  For our problem, 
⟨n (l^2) ((, )^(2  S + 
                      1)) L | 
              ((e^2)/(Subscript[r, 12])) | n (l^2) ((, )^(2  S + 1)) 
                L⟩=(∑_λ)(F^λ)(n l,  
                n l)(((-1))^L)(((2  l + 1))^2)(((({{l, l, λ}, {0, 0, 0}})))^2){({{l, l, L}, {l, l, λ}})}
This fact also makes it clear why the result doesn't depend on (Subscript[M, L]).

Matt Cargo lectured on the spin network which gives this expression by inspection *right away*.  Ask him for references.

The expression can be evaluated readily by Mathematica

(((-1))^L) (((2  l + 1))^2) ThreeJSymbol[{l, 0}, {l, 
              0}, {λ, 0}]^2 
      SixJSymbol[{l, l, L}, {l, l, λ}]

(((-1))^0) (((2*1 + 1))^2) ThreeJSymbol[{1, 0}, {1, 
              0}, {0, 0}]^2 SixJSymbol[{1, 1, 0}, {1, 1, 0}]

(((-1))^0) (((2*1 + 1))^2) ThreeJSymbol[{1, 0}, {1, 
              0}, {2, 0}]^2 SixJSymbol[{1, 1, 0}, {1, 1, 2}]

(((-1))^2) (((2*1 + 1))^2) ThreeJSymbol[{1, 0}, {1, 
              0}, {0, 0}]^2 SixJSymbol[{1, 1, 2}, {1, 1, 0}]

(((-1))^2) (((2*1 + 1))^2) ThreeJSymbol[{1, 0}, {1, 
              0}, {2, 0}]^2 SixJSymbol[{1, 1, 2}, {1, 1, 2}]

(((-1))^1) (((2*1 + 1))^2) ThreeJSymbol[{1, 0}, {1, 
              0}, {0, 0}]^2 SixJSymbol[{1, 1, 1}, {1, 1, 0}]

(((-1))^1) (((2*1 + 1))^2) ThreeJSymbol[{1, 0}, {1, 
              0}, {2, 0}]^2 SixJSymbol[{1, 1, 1}, {1, 1, 2}]

(b) comparison to data

Group IV

In this group, there are two electrons in a p orbital.  We take the average of levels split by the spin-orbit coupling weighted by the number of states 2  J + 1.  The predicted ratio of energy splittings is
\(TraditionalForm\`\(\^1 S\( - \^1\)D\)\/\(\^1 D\( - \^3\)P\) = \(\(2\/5 - 1\/25\)\/\(1\/25 + 1\/5\) = 3\/2\)\).
On the other hand, the data says

C:	1.13
Si:	1.48
Ge:	1.50
Sn:	1.39
Pb:	1.62

The data agree very well with the prediction for Si and Ge.  It fails for Sn and Pb, namely for very large atoms, because the spin-orbit coupling becomes important as ∝(((Z α))^2)so that the single-particle wave function is no longer described well by separate l and s quantum numbers.  Therefore the LS coupling scheme ceases to be a good approximation.  On the other hand, the prediction fails for C, too.  This is because two electrons we consider have a large enough impact on the remaining four electrons.  Even though the four electrons form a closed shell, it is so small that it gets disturbed by the high cost of excitation due to the two outer-shell electrons.  Formally this effect is at the second-order in perturbation theory while we considered only the first-order effect (namely just the expectation values).

Group VI

In this group, there are two holes in a p orbital.  The data says
  
O: 	1.14
S: 	1.43
Se:	1.49
Te:	1.49
Po:	no data

Again the agreement is very good for the intermediate size atoms, while it fails for the smallest one, O.  The reason is expected to be the same as the Group IV.

What the pros do

It is instructive to look at theretical predictions (postdictions) by pros who use multi-configuration Hartree-Fock.  Such calculations can deal with both the impact of outer shell electrons on the closed shells and/or large relativistic effects for large Z.  I copy results by Charlotte Froese Fischer and co. from http://atoms.vuse.vanderbilt.edu/.  For C, where the impact of the two 2  p electrons on the inner shell appear important, the comparison to the data is
level		data		MCHF
\(TraditionalForm\`\(\(\^3\)\(P\_0\)\)\)		0		0
\(TraditionalForm\`\(\(\^3\)\(P\_1\)\)\)		16.41671	16.33
\(TraditionalForm\`\(\(\^3\)\(P\_2\)\)\)		43.41350	43.03
\(TraditionalForm\`\(\(\^1\)\(D\_2\)\)\)		10192.66	10268.23
\(TraditionalForm\`\(\(\^1\)\(S\_0\)\)\)		21648.02	21818.60
Impressive.

2. Thomas-Fermi energy levels

The purpose of this problem is to understand qualitatively the single-particle energy levels in multi-electron atoms.  Due to the screening of the nuclear charge by the inner-shell electrons, the degeneracy among the orbitals of the same principal quantum number n is lifted.  This explains the empirical fact of the "Aufbau principle", namely how to build atoms filling the orbitals in the order of 1  s, 2  s, 2  p, 3  s, 3  p, 4  s, 3  d, 4  p, 5  s, 4  d, 5  f, 5  p, 6  s, etc

For a useful comparison, the excited levels of Al provides the energy levels above 3  p, as single-electron levels in the potential provided by the 1  (s^2) 2  (s^2) 2  (p^6) 3  (s^2) closed shell.

Thomas-Fermi for \(TraditionalForm\`\(Al\^+\)\)

I have given the following note to you, and you didn't have to work this out.

Following the discusions on the positive ions in the lecture notes, we need to find χ'(0) so that (Subscript[x, 0]) χ' ((Subscript[x, 0])) = -(Z - N)/
            Z = -1/13 for \(TraditionalForm\`\(Al\^+\)\).  After some trial-and-error, I find χ' (0) = -1.5881063075 works:

chisol = 
      NDSolve[{(1 + (4/3) (x^(3/2))) y''[x] + 
              4  (x^(1/2)) y'[x] + 
              If[x == 0, 
                0, y[x]/(x^(1/2)) (1 - (((1 + (4/3) 
                                  (x^(3/2))))^(3/2)) 
                        Max[y[x], 0]^(1/2))] == 0, 
          y[0] == 1, y'[0] == -1.5881063075}, y, {x, 0, 300}]; χ[
        x_] := (1 + 4/3*(x^(3/2))) y[x]; 
    Plot[χ[x] /. chisol, {x, 0, 300}]

Plot[χ[x] /. chisol, {x, 0, 30}]

Look for where χ((Subscript[x, 0])) = 0,

(Subscript[x, 0]) = x /. FindRoot[χ[x] /. chisol, {x, 0, 100}]

Now make sure (Subscript[x, 0]) χ' ((Subscript[x, 0])) = -1/Z,

Z x χ'[x] /. 
          chisol /. {x -> (Subscript[x, 0])} /. {Z -> 13}

Using the definition e Φ(
            r) = ((Z (e^2))/r) χ(r) and variables r =  (Z^(-1/3)) b x, b = (1/2) (((((3  π)/4)))^(2/
                3)) ((ℏ^2)/(m (e^2))), we define (in the atomic unit e = m = ℏ = 1)

Φ[
        r_] := -(Z/r) χ[((((Z^(-1/
                                3)) (1/2) (((((3  π)/4)))^(2/3))))^(-1)) 
                  r] /. {(Subscript[r, 0]) -> (Z^(-1/
                        3)) (1/2) (((((3  π)/4)))^(2/3)) 
                  (Subscript[x, 0])} /. chisol[[1]] /. {Z -> 13}

Plot[r Φ[r], {r, 0, 20}]

The Fermi energy is just (Subscript[ϵ, F]) = ((-(e^2))/(Subscript[r, 0])), where (Subscript[r, 0]) = (Z^(-1/3)) 
            b (Subscript[x, 0]) = (Z^(-1/
                  3)) (1/2) (((((3  π)/4)))^(2/3)) 
            (Subscript[x, 0]) in the atomic unit.

ϕ[
        r_] := -(Z/r) χ[((((Z^(-1/
                                  3)) (1/2) (((((3  π)/4)))^(2/3))))^(-1)) r] - 
              (1/(Subscript[r, 0])) /. {(Subscript[r, 0]) -> (Z^(-1/
                        3)) (1/2) (((((3  π)/4)))^(2/3)) 
                  (Subscript[x, 0])} /. chisol[[1]] /. {Z -> 13}

(Z^(-1/3)) (1/2) (((((3  π)/4)))^(2/3)) 
        (Subscript[x, 0]) /. {Z -> 13}

Plot[r ϕ[r], {r, 0, 10}]

Indeed, ϕ(r) becomes (-1) ((e^2)/r) due to the overall positive charge beyond r > (Subscript[r, 0]) where the electron number density vanishes.  The coefficient is consistent with -1 within the tiny numerical errors.  To make this point clear, let us plot it on bigger scales,

Plot[r ϕ[r], {r, 0, 10}, PlotRange -> {-14, 0}, 
      PlotStyle -> RGBColor[0, 1, 0]]

r ϕ(r) basically is the total charge contained in the radius r.  At the origin, it is 13. Away from the origin, it asymptotes to 1, which is the charge of \(TraditionalForm\`\(Al\^+\)\).  

Therefore, for orbitals close to the origin (like 1  s), we expect it the energy is close to the hydrogenic spectrum with Z = 13, while for high orbitals far away from the origin (like 5  g), the energy is close to that with Z = 1.  This qualitative understanding helps to start guessing the energy eigenvalues.

l=0

no nodes (1s)

With no nodes (ground state).  For the hydrogen-like spectrum, we expect ε = -(13^2) = -169.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                0} /. {ε -> -103.77652277489802}, ρ, {r, 0, 5}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 5}, 
      PlotRange -> {-1, 1}]

one node (2s)

Now for the higher levels.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                0} /. {ε -> -7.377516185}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-0.2, 0.2}]

two nodes (3s)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                0} /. {ε -> -0.55563}, ρ, {r, 
          0, 20}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 20}, 
      PlotRange -> {-0.2, 0.2}]

three nodes (4s)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                0} /. {ε -> -0.1746}, ρ, {r, 0, 
          30}]; Plot[Evaluate[ρ[r] /. sol], {r, 0, 30}, 
      PlotRange -> {-0.1, 0.1}]

four nodes (5s)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                0} /. {ε -> -0.086365843}, ρ, {r, 0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.03, 0.03}]

five nodes (6s)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                0} /. {ε -> -0.0515225}, ρ, {r, 
          0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.02, 0.02}]

l=1

no nodes (2p)

With no nodes (ground state), we expect ε = ((-169)/4) = -42.25.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                1} /. {ε -> -4.6442405}, ρ, {r, 
          0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-0.01, 0.1}]

one node (3p)

Now for the higher levels.  With one node, we expect ε = ((-1)/9),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                1} /. {ε -> -0.27415}, ρ, {r, 
          0, 30}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 30}, 
      PlotRange -> {-0.015, 0.015}]

two nodes (4p)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                1} /. {ε -> -0.116966184}, ρ, {r, 0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.015, 0.015}]

three nodes (5p)

With two nodes, we expect ε = ((-1)/25),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                1} /. {ε -> -0.0648106}, ρ, {r, 
          0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.001, 0.001}]

four nodes (6p)

With two nodes, we expect ε = ((-1)/36),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                1} /. {ε -> -0.04115}, ρ, {r, 
          0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.0003, 0.001}]

l=2

no nodes (3d)

With no nodes (ground state), we expect ε = ((-169)/9) = -18.78

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                2} /. {ε -> -0.112526525}, ρ, {r, 0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.001, 0.01}]

one node (4d)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                2} /. {ε -> -0.0633465}, ρ, {r, 
          0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.0004, 0.0003}]

two nodes (5d)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                2} /. {ε -> -0.0405}, ρ, {r, 0, 
          100}]; Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.0004, 0.0003}]

l=3

no nodes (4f)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                3} /. {ε -> -0.062517}, ρ, {r, 
          0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.001, 0.01}]

one node (5f)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                3} /. {ε -> -0.04002}, ρ, {r, 
          0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.0003, 0.0003}]

l=4

no nodes (5g)

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[r] - 
                      2  ϕ[r] ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-13)/(l + 1))} /. {l -> 
                4} /. {ε -> -0.04002}, ρ, {r, 
          0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.001, 0.01}]

Coparison to data

Values of ε calculated above are 2  E in the atomic unit (Hartree).  Dividing them by two:

\({\(-103.77652277489802\), 
\(-7.377516185\), 
\(-4.6442405\), 
\(-0.55563\), 
\(-0.27415\), 
\(-0.1746\), 
\(-0.116966184\), 
\(-0.112526525\), 
\(-0.086365843\), 
\(-0.0648106\), 
\(-0.0633465\), 
\(-0.062517\), 
\(-0.0515225\), 
\(-0.04115\), 
\(-0.0405\), 
\(-0.04002\), 
\(-0.04002\)}/2\)

ordering is
1s	-51.8883
2s	-3.68876
2p	-2.32212
3s	-0.277815
3p	-0.137075
4s	-0.0873
4p	-0.0584831
3d	-0.0562633
5s	-0.0431829
5p	-0.0324053
4d	-0.0316733
4f	-0.0312585
6s	-0.0257613
6p	-0.020575
5d	-0.02025
5f	-0.02001
5g	-0.02001

On the other hand, the data, after taking weighted average over the spin-orbit splittings and converting (cm^(-1)) to Hartree, is
3p	-0.21975
4s	-0.10454
3d	-0.087443
4p	-0.072221
5s	-0.048273
4d	-0.042608
5p	-0.036482
4f	-0.031725
6s	-0.027964
4d	-0.027545
6p	-0.022527
5f	-0.020275
5g	-0.020071
5d	-0.018739

Note that "4d" is repeated twice.  This is because 3  s 3  (p^2)(D) configuration gives exactly the same quantum number \(TraditionalForm\`\(\^2\) D\) as the 3  (s^2) n d configurations, and hence effectively inserts another state into the 3  (s^2) n d sequence.  V. Kaufman and W.C. Martin, J. Phys. Chem. Ref. Data, 20, 775 (1991) decided to use "4d" label twice so that at very high n d they line up well with n f and n g orbitals.  This issue also explains why the n d states show poor agreement beween the calculation and data.  Apparently this new state is somewhere between our 3  d and 4  d.  Due to the no-level-crossing theorem, the states repel by mixing, and the 3  d state is pushed down, while the normal 4  d state is pushed up, and the new "4d" state is inserted between them.  The mixing is less important as you go to higher states indeed, and the 5d state shows reasonable agreement.  Nonetheless, the data shows the 5d state is above 5f and 5g as a consequence of this mixing.

The graphical comparison between this calculation and the data is shown at http://hitoshi.berkeley.edu/221B/Aldata.pdf

Even though the details are different, the Thomas-Fermi calculation shows important qualitative features that are consistent with data:
(1) higher l orbitals have higher energies for the same principal quantum number n
(2) 4  p and 3  d are very close
(3) 5  p, 4  d, and 4  f are very close
(4) 6  p, 5  d, 5  f, and 5  g are very close

The quantiative agreement is something like 20% level,which I find remarkble for a crude model as Thomas-Fermi and Z = 13 which is certainly not very close to infinity.

Ask the Pros

Again consulting http://atoms.vuse.vanderbilt.edu/, here is the comparison in (cm^(-1))

configuration	state	MCHF		data
3  (s^2) 3  p	\(TraditionalForm\`\(\(\ \^2\)\(P\_\(1/2\)\%o\)\)\) 	0.0000	
		\(TraditionalForm\`\(\(\ \^2\)\(P\_\(3/2\)\%o\)\)\)	90.10  		112.061

3  (s^2) 4  s	\(TraditionalForm\`\(\(\ \^2\)\(S\_\(1/2\)\)\)\)	25582.85	25347.756

3  s 3  (p^2)	\(TraditionalForm\`\(\(\ \^4\)\(P\_\(1/2\)\)\)\)	28 080.16	29020.41
		\(TraditionalForm\`\(\(\ \^4\)\(P\_\(3/2\)\)\)\)	28 117.91	29066.96
		\(TraditionalForm\`\(\(\ \^4\)\(P\_\(5/2\)\)\)\)	28 178.99	29142.78
		
3  (s^2) 3  d	\(TraditionalForm\`\(\(\^2\)\(D\_\(3/2\)\)\)\)	32184.61	32435.453
		\(TraditionalForm\`\(\(\^2\)\(D\_\(5/2\)\)\)\)	32186.17	32436.796

3  (s^2) 4  p	 \(TraditionalForm\`\(\(\ \^2\)\(P\_\(1/2\)\%o\)\)\)	33047.68	32949.807
                         \(TraditionalForm\`\(\(\ \^2\)\(P\_\(3/2\)\%o\)\)\) 33060.38	32965.639

The list doesn't go up to as high levels as ours do. But it is impressive.

3. Thomson's Plum Pudding Model

This problem has an obvious historical importance.  The Rutherford's experiment discriminated the "Rutherford model of atoms" against the "Thomson's Plum Pudding Model of atoms."  Therefore we need to see what the predicted difference between two models is for the Rutherford's experiment.

First the Born approximation.  Using the form factor defined in the lecture notes (Eq. (16)),
\(TraditionalForm\`F(
          q) = \(1\/Z\) \(∫d\ \(x\& -> \) \(\(ρ\_N\)(\(x\& -> \))\)\ e\^\(i\ \(q\& -> \)·\(x\& -> \)\)\)\).  For this problem, \(TraditionalForm\`\(ρ\_N\)(\(x\& -> \)) = \(Z\/\((2  π\ σ)\)\^\(3/2\)\) e\^\(\(-\(x\& -> \)\^2\)/2  σ\^2\)\).  Therefore,
\(TraditionalForm\`F(
          q) = \(1\/Z\) \(∫d\ \(x\& -> \) \(Z\/\((2  π\ σ)\)\^\(3/2\)\) 
              e\^\(\(-\(x\& -> \)\^2\)/2  σ\^2\)\ e\^\(i\ \(q\& -> \)·\(x\& -> \)\)\)\)
\(TraditionalForm\`\(\(=\)\(∫d\ \(x\& -> \) \(1\/\((2  π\ σ)\)\^\(3/2\)\) 
            e\^\(\(-\((\(x\& -> \) - i\ \(σ\^2\) \(q\& -> \))\)\^2\)/2  σ\^2 - \(σ\^2\) \(q\& -> \)\^2/2\)\)\)\)
=(e^(-(σ^2) (((Overscript[q,  ->]) )^2)/
                2))

Using Eqs. (17) and (12), and (q^2) = 2  (k^2)(1 - cos θ) = 
          4  (k^2) (sin^2)(θ/2),
((d σ)/(d Ω))=\(\(\(\((2  
                            m)\)\^2\) \((Z\ Z'\ e\^2)\)\^2\)\/\(16 \(\((ℏ\ k)\)\^2\) \(\(sin\^4\)(θ/2)\)\)\)\(e\^\(\(-4\) \(σ\^2\) \(k\^2\) \(\(sin\^2\)(θ/2)\)\)\).

Now we can plot both the point-like case and the plum pudding case.  We take E = 10  MeV = 1.60 (10^(-5)) dyn for the α-particle, and σ≈0.528 Å.  The mass of the α-particle is 4.002602u, where 1  u = 1.6605×(10^(-24)) g. The momentum of the α-particle is p = (Sqrt[2  (Subscript[m, α]) E]) = 
          1.46×(10^(-14)) g cm/sec.

(Sqrt[2*4.002602*1.6605 (10^(-24))*1.60 (10^(-5)) ])

(((((2  m))^2) (((Z Z' (e^2)))^2))/(16 (((ℏ k))^4) Sin[θ/2]^4)) /. {k -> 1.46 (10^(-14))/ℏ, 
        Z -> 79, Z' -> 2, 
        e -> (Sqrt[(1/137) 197 (10^6) (10^(-13)) 1.6 (10^(-12))]), m -> 9.11 (10^(-28))}

<< Graphics`Graphics`

Here is the differential cross section by the point-like nucleus as a function of the scattering angle θ,

LogPlot[
      6.034076289556489`*^-33 Csc[(θ/2)]^4, {θ, 
        0, π}]

Here is the form factor

(E^(-4 (σ^2) (k^2) Sin[θ/2]^2)) /. {σ -> 0.528 (10^(-8)), 
          k -> 1.46 (10^(-14))/ℏ} /. {ℏ -> 
          1.054 (10^(-27))}

Multiplying the form factor on the point-like differential cross section

LogPlot[
      6.034076289556489`*^-33 Csc[(θ/2)]^4 (E^(-2.1396972386751106`*^10 Sin[(θ/2)]^2)), {θ, 0, (π/16)}]

Here is the comparison of differential cross sections between the Rutherford's and Thomson's models.

Show[%%%, %]

Namely, the Gaussian form factor drops off so quickly that there is basically no cross section at all at any large angles, not to mention backwards (θ=π).  This comparison makes it clear why people had to abandon Thomson's model after the experiment saw the backward scattered α-particle.