HW #7

1. Fermi energy

This is a problem meant to remind ourselves Fermi-degenerate gas which is assumed in the Thomas-Fermi model of atoms, not a realistic condensed-matter problem which is beyond the scope of this course.  In any case, here we go.

The number density in the Fermi-degenerate gas is
\(TraditionalForm\`n = \(g \(∫\^p\_F\(\( d\^3\) p\)\/\((2  π\ ℏ)\)\^3\) = 
          g  1\/\((2  π\ ℏ)\)\^3\ \(\(4  π\)\/3\) 
            p\_F\^3\)\)
where g = 2 for the spin degrees of freedom for electrons, and \(TraditionalForm\`n = \(a\^\(-3\) = \((3.61  Å)\)\^\(-3\)\)\) for copper.  We then find
\(TraditionalForm\`p\_F = \(\(2  π\ ℏ\)\/a\) \((3\/\(4`  π\ g\))\)\^\(1/3\)\), and hence the Fermi energy
\(TraditionalForm\`E\_F = \(p\_F\^2\/\(2  m\) = \(1\/\(2  
                  m\)\) \(\((\(2  π\ ℏ\)\/a)\)\^2\) \((3\/\(4`  π\ g\))\)\^\(2/3\)\)\)

(1/(1.60 (10^(-19)))) (1/(2 
               (Subscript[m, e]))) (((((2  π ℏ)/a)))^2) ((((3/(4  π g))))^(2/3)) /. {ℏ -> 1.054 (10^(-34)), 
        (Subscript[m, e]) -> 9.11 (10^(-31)), a -> 3.61 (10^(-10)), 
        g -> 2}

(Subscript[E, F]) = 2.80 eV.  This is much higher than the room temperature, k T~(1/40) eV.  Therefore, the room temperature is very *cold* for electrons in a metal.

I'm not asking for more in this problem, but to be more realistic, we need to know that the copper is not a simple cubic lattice but it is a face-centered cubic, which makes the density four times larger.  Then the four momentum larger by (4^(1/3)) and the Fermi energy by (4^(2/3)).

(4^(2/3)) (1/(1.60 (10^(-19)))) (1/(2 
               (Subscript[m, e]))) (((((2  π ℏ)/a)))^2) ((((3/(4  π g))))^(2/3)) /. {ℏ -> 1.054 (10^(-34)), 
        (Subscript[m, e]) -> 9.11 (10^(-31)), a -> 3.61 (10^(-10)), 
        g -> 2}

The experimental value is 7.00 eV, slightly lower (see, e.g., http://hyperphysics.phy-astr.gsu.edu/hbase/tables/fermi.html).

At this level of accuracy, one has to consider the fact that the electrons are not entirely free, but interact with the ions and hence the lattice.  This is the electron-phonon interaction.  It is know to result in an "electronic mass enhancement", which reduces the Fermi energy correspondingly.  For example, a paper F.S. Khan, P.B. Allen, W.H. Butler, and F.J. Pinski, Phys. Rev. B26, 1538-1548 (1982) calculates the effect.

2. The Thomas-Fermi Model.

(a) Solving for χ(x)

Following the technical remark in the lecture notes we substitute χ(x) = (1 + (4/3) (x^(3/2))) y(x)and solve for y(x). The initial condition for χ ' (0) is found by playing around with its value until we get the desired boundary condition χ(x -> ∞) = 0. The solution is extremely sensitive to χ '(0), so we need to do a trial and error to figure out the correct value (unfortunately, the value may be version and platform dependent due to numerical inaccuracies).

sol = 
      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)) y[x]^(1/2))] == 
            0, y[0] == 1, y'[0] == -1.5880709987255}, 
        y, {x, 0, 300}]; χ[x_] := (1 + 4/3*(x^(3/2))) y[x]; 
    Plot[Evaluate[χ[x] /. sol], {x, 0, 300}, 
      PlotRange -> {-0.01, 0.01}]

Plot[Evaluate[χ[x] /. sol], {x, 0, 100}, 
      AxesLabel -> {x, χ}, PlotRange -> {0, 1}]

Plot[Evaluate[χ[x] /. sol], {x, 0, 10}, 
      AxesLabel -> {x, χ}, PlotRange -> {0, 1}]

(b)

If we neglect inter-electron interaction we know that to leading order (in α) the energy depends only on the quantum number n and not on l. How will the effect of screening change this? 
Zeχ(r) is the effective charge an electron at a distance r from the nucleus see's through the screening of the electrons closer to the nucleus. An electron that resides in an orbital that is further away from the nucleus will feel less of the nuclear charge (by the shape of χ) and will therefore be less tightly bound, or with a higher energy. Clearly electrons with higher l spend most of their time further away from the nucleus (semi-classically). In other words, l = 0 corresponds to radial motion, probing closer to the nucleus, and l = n corresponds to the most circular motion, staying further away form the nucleus. Intermediate values of l correspond to an elliptical orbit. Thus, electrons with higher l will have higher energy.

(c) The Total Binding Energy

The total energy of the atom is the sum of the kinetic and the potential energies.  

We cannot simply compute -(1/2) ∫e ρ(
              r) ϕ(r) , with ϕ(r) the Thomas-Fermi  potential, because this would only count half the electron-nuclear interaction energy since -e ρ does not include the nuclear charge density .  If we tried to remedy this by calculating (1/2) ∫ (-e ρ(r)  + 
                Z e (δ^3)((Overscript[r, ⇀]))) ϕ(
              r), we would have included too much--our  result would diverge because of the self-energy of the nuclear point charge.  Said once more, the potential energy is almost (1/2) ∫ (-e ρ (r) + 
                Z e (δ^3)((Overscript[r, ⇀]))) ϕ(
              r) = \(TraditionalForm\`1\/2\ \(∫\ \((\(-e\)\ ρ\ \((r)\) + 
                Z\ e\ \(\(δ\^3\)(
                    r\&⇀)\))\)\ \(\((\(-e\)\ \(ρ(r')\)\  + 
                    Z\ e\ \(\(δ\^3\)(
                        r\&⇀')\))\)\/\(\(|\)\(\ \)\(r\&⇀\  - \ 
                    r\&⇀'\)\(\ \)\(|\)\)\) \(d\^3\) \(r\&⇀\) \(d\^3\) r\&⇀'\)\), but we must subtract off the nucleus-nucleus self-energy  piece which is the product of the two δ-functions.  Left over are exactly the electron-electron  and electron-nucleus pieces mentioned below.

Kinetic energy

To calculate the kinetic energy we assume (as in the notes) the electrons are distributed evenly in phase space. The number density in phase space is
n(r, p)={({{(2/(((2π))^3)), for (((p^2) )/(2  m)) < 
                    eΦ}, {0, \(for\ \(\(p\^2\)\(\ \)\)\/\(2  m\) > \(\(eΦ\)\(.\)\)\)}})
So the kinetic energy is 
(Subscript[E, K])=∫(d^3) p (d^3) 
                r ((p^2)/(2  m)) n(p, r)=(2/(((2π))^3))(∫_0^∞)dr (r^2)(∫_0^((Subscript[p, f])))dp (p^2)((p^2)/(2  m))
=(1/((((2π))^3)m))∫dr(1/5)(((2me Φ))^(5/
                      2)).
Plugging in Φ = (Ze/r) χ(r) and changing the integration variable (Overscript[r, ⇀])  =  
        (Overscript[x, ⇀]) (Z^(-1/3)) b,as defined in the notes with b = (1/2) (((((3  π)/4)))^(2/
                3)) ((ℏ^2)/(m (e^2))) we get 
(Subscript[E, K])=(3/5)(((Z^(7/3)) (e^2))/b)∫((χ^(5/2))/(x^(1/2)))dx.

Electron-Nucleus potential energy

The potential energy is
\(TraditionalForm\`E\_\(e\ A\) = \(-\(∫d\ \(r\& -> \) \(ρ(r)\)\ \(Z\ e\^2\)\/r\)\)\)
Changing the variable  r = ((Z^(-1/3)) b) x, and  ρ(
          r) = ((Z^2)/(4  π (b^3))) ((((χ/x)))^(3/2)),
\(TraditionalForm\`E\_\(e\ A\) = \(-\(\(\(Z\^\(7/3\)\) 
                  e\^2\)\/b\)\) \(∫d\ x\ \(χ(x)\)\^\(3/2\)\/x\^\(\(1/2\)\(\ \)\)\)\).

Electron-Electron potential energy

The potential energy is
\(TraditionalForm\`E\_\(e\ e\) = \(1\/2\) \(∫d\ \(\(r\& -> \)\_1\) \(∫d\ \(\(r\& -> \)\_2\) \(ρ(
                    r\_1)\) \(e\^2\/\(\(|\)\(\(r\& -> \)\_1 - \(r\& -> \)\_2\)\(|\)\)\) \(ρ(r\_2)\)\)\)\).
The factor of a half is needed to do the double counting of the Coulomb repulsion.  It is the continuum generalization of \(∑\_\(i < j\)e\^2\/r\_\(i\ j\) = \(1\/2\) ∑\_\(i, j\)\)((e^2)/(Subscript[r, i j])), where the term i = j cannot be avoided for a continuous medium.

Expanding the Coulomb potential in spherical harmonics, only the l = 0 term contributes because ρ(r) is spherically symmetric.  Therefore,
(Subscript[E, e e])
\(TraditionalForm\`\(\(=\)\(\(1\/2\) \(∫\(r\_1\^2\) 
              d\ \(r\_1\) 
              d\ \(Ω\_1\) \(∫\(r\_2\^2\) d\ \(r\^2\) 
                  d\ \(Ω\_2\) \(ρ(r\_1)\) \(ρ(
                    r\_2)\) \(e\^2\/\(r\_ > \)\) 
                  4  π\ \(\(Y\_0\%0\)(Ω\_1)\) \(\(Y\_0\%0\)(Ω\_2)\)\)\)\)\)\)
\(TraditionalForm\`\(\(=\)\(\(e\^2\/2\) \(\((4  π)\)\^2\) \(∫d\ \(r\_1\) \(r\_1\^2\) \(∫d\ \(r\_2\) \(r\_2\^2\) \(ρ(r\_1)\) \(ρ(r\_2)\) 1\/\(r\_ > \)\)\)\)\)\)
\(TraditionalForm\`\(\(=\)\(\(\(\(e\^2\)(
                4  π)\)\^2\) \(∫\_0\%∞ 
              d\ \(r\_1\) \(r\_1\^2\) \(∫\_0\%\(r\_1\)d\ \(r\_2\) \(r\_2\^2\) \(ρ(r\_1)\) \(ρ(r\_2)\) 1\/r\_1\)\)\)\)\)
Now changing the variables,
\(E\_\(e\ e\) = \(\(\(e\^2\)(
                  4  π)\)\^2\) \(\((\(Z\^\(\(-1\)/3\)\) 
                    b)\)\^5\) \(\((Z\^2\/\(4  π\ b\^3\))\)\^2\) \(∫\_0\%∞ 
                d\ \(x\_1\) \(x\_1\^2\) \(∫\_0\%\(x\_1\)d\ \(x\_2\) \(\(\(x\_2\^2\)(\(χ(x\_1)\)\/x\_1)\)\^\(3/
                          2\)\) \(\((\(χ(x\_2)\)\/x\_2)\)\^\(3/2\)\) 
                    1\/x\_1\)\)\)
=(Z^(7/3))((e^2)/b)(∫_0^∞)d (Subscript[x, 1])(((χ((Subscript[x, 1])))^(3/2))/((Subscript[x, 1])^(1/2)))(∫_0^((Subscript[x, 1])))d (Subscript[x, 2])((Subscript[x, 2])^(1/2))((χ(
                          (Subscript[x, 2])))^(3/2))
Using the differential equation χ'' = ((χ^(3/2))/(x^(1/2))), we can simplify it drastically,
(Subscript[E, e e]) = (Z^(7/
                3)) ((e^2)/b) (∫_0^∞) 
              d (Subscript[x, 1]) χ'' ((Subscript[x, 1])) (Subscript[((∫_0)^x), 1]) 
                  d (Subscript[x, 2]) (Subscript[x, 2]) χ'' ((Subscript[x, 2]))
Integrating it in parts in (Subscript[x, 1]),
\(TraditionalForm\`E\_\(e\ e\) = \(Z\^\(7/
                3\)\) \(e\^2\/b\) {\([χ' \((x\_1)\) \(∫\_0\^x\_1 d\ \(x\_2\) \(x\_2\) χ'' \((x\_2)\)\)]\)\_0\%∞ - ∫\_0\%∞ 
                  d\ \(x\_1\) χ' \((x\_1)\) \(x\_1\) χ'' \((x\_1)\)}\)
The first term vanishes at the boundaries, while the second term can be integrated in parts again,
\(TraditionalForm\`E\_\(e\ e\) = \(Z\^\(7/
                3\)\) \(e\^2\/b\) {\(-\([x\ \(1\/2\) χ' \((x)\)\^2]\)\_0\%∞\) + ∫\_0\%∞ 
                  d\ x\ \(1\/2\) χ' \((x)\)\^2}\)
\(TraditionalForm\`\(\(=\)\(\(1\/2\) \(Z\^\(7/
                3\)\) \(e\^2\/b\) \(∫\_0\%∞ χ' \((x)\)\^2\)\)\)\).

The total

We have found
(Subscript[E, K])=(3/5)(((Z^(7/3)) (e^2))/b)(∫_0^∞)((χ^(5/2))/(x^(1/2)))dx
(Subscript[E, e A]) = -(((Z^(7/3)) 
                  (e^2))/b) (∫_0^∞) 
              d x (((χ(x))^(3/2))/(x^(1/2 )))
\(TraditionalForm\`E\_\(e\ e\) = \(1\/2\) \(Z\^\(7/
                3\)\) \(e\^2\/b\) \(∫\_0\%∞ χ' \((x)\)\^2\)\)

Numerical integrals

We can let Mathematica do these integrals numerically.

\(\(3\/5\) 
      NIntegrate[
        Evaluate[χ[x]\^\(5/2\)\/x\^\(1/2\) /. sol], {x, 0, 
          300}]\)

\(\(-NIntegrate[
        Evaluate[χ[x]\^\(3/2\)\/x\^\(1/2\) /. sol], {x, 0, 
          300}]\)\)

(1/2) 
      NIntegrate[Evaluate[χ'[x]^2 /. sol], {x, 0, 300}]

% + %%

Therefore,
(Subscript[E, K]) = 0.680602 (Z^(7/3)) ((e^2)/b)
(Subscript[E, e A]) = -1.58807 (Z^(7/3)) ((e^2)/b)
(Subscript[E, e e]) = 0.226867 (Z^(7/3)) ((e^2)/b)
and the total potential energy is
(Subscript[E, e A]) + 
          (Subscript[E, e e]) = -1.3612 (Z^(7/2)) 
            ((e^2)/b) = -2 (Subscript[E, K])
Numerically, the total binding energy is 
|(Subscript[E, K]) + (Subscript[E, e A]) + 
            (Subscript[E, e e])| = (Subscript[E, K]) = 20.92  eV (Z^(7/2))

(In general, it is a bad idea to perform a numerical integration for an integrand that is singular at the boundary (e.g., (x^(-1/2))), but apparently Mathematica handles it very well.)

Analytic method

Actually, the integrals can be performed analytically using the fact that \(TraditionalForm\`χ'' = χ\^\(3/2\)\/x\^\(\(1/2\)\(\ \)\)\), χ(0) = 1, and χ(∞) = 0.  We start with the electron-nucleus potential energy,
\(TraditionalForm\`E\_\(e\ A\) = \(\(-\(\(\(Z\^\(7/3\)\) 
                    e\^2\)\/b\)\) \(∫\_0\%∞ 
                d\ x\ \(χ(x)\)\^\(3/2\)\/x\^\(\(1/2\)\(\ \)\)\) = \(-\(\(\(Z\^\(7/3\)\) e\^2\)\/b\)\) \(∫\_0\%∞ 
                d\ x\ χ'' \((x)\)\)\)\)
\(TraditionalForm\`\(\(=\)\(\(-\(\(\(\(Z\^\(7/3\)\) e\^2\)\/b\)[χ' \((x)\)]\)\_0\%∞\) = \(-\(\(\(Z\^\(7/3\)\) 
                    e\^2\)\/b\)\) \((\(-χ'\) \((0)\))\)\)\)\)
The other two terms require more technical manipulations.  We first consider the following integral that appears in the kinetic energy,
\(TraditionalForm\`∫\_0\%∞\( χ\^\(5/2\)\/x\^\(1/2\)\) 
            d\ x = \(∫\_0\%∞ χ\ \(χ\^\(3/2\)\/x\^\(1/2\)\) d\ x = ∫\_0\%∞ χ\ χ''\ d\ x\)\)
\(TraditionalForm\`\(\(=\)\(\([χ\ χ']\)\_0\%∞ - ∫\_0\%∞ χ' χ' 
                d\ x = \(-χ'\) \((0)\) - ∫\_0\%∞\( χ'\^2\) d\ x\)\)\)
On the other hand, there is another way to do the integration by parts,
\(TraditionalForm\`∫\_0\%∞\( χ\^\(5/2\)\/x\^\(1/2\)\) 
            d\ x = \([2 \( x\^\(1/2\)\) χ\^\(5/2\)]\)\_0\%∞ - ∫\_0\%∞ 
              2 \( x\^\(1/2\)\) \(5\/2\) \(χ\^\(3/2\)\) χ'\ d\ x\)
\(TraditionalForm\`\(\(=\)\(\(-5\) \(∫\_0\%∞\( 
                    x\^\(1/2\)\) \(χ\^\(3/
                      2\)\) χ'\ d\ x\) = \(-5\) \(∫\_0\%∞ x\ χ''\ χ'\ d\ x\)\)\)\)
\(TraditionalForm\`\(\(=\)\(\(-5\) \(∫\_0\%∞ 
                x\ \(d\/\(d\ x\)\) \((\(1\/2\) χ'\^2)\) 
                d\ x\) = \(-5\) {\([x\ \(1\/2\) χ'\^2]\)\_0\%∞ - ∫\_0\%∞\( 1\/2\) \(χ'\^2\) d\ x}\)\)\)
\(TraditionalForm\`\(\(=\)\(\(5\/2\) \(∫\_0\%∞\( χ'\^2\) d\ x\)\)\)\)
Comparing the two equations, we find
\(TraditionalForm\`∫\_0\%∞\( χ\^\(5/2\)\/x\^\(1/2\)\) d\ x = \(5\/7\) \((\(-χ'\) \((0)\))\)\)
\(TraditionalForm\`∫\_0\%∞\( χ'\^2\) 
            d\ x = \(2\/7\) \((\(-χ'\) \((0)\))\)\)
Therefore,
(Subscript[E, K]) = 
        (3/7) (Z^(7/3)) ((e^2)/b) (-χ' (0))
(Subscript[E, e A]) = - 
            (Z^(7/3)) ((e^2)/b) (-χ' (0))
(Subscript[E, e e]) = 
        (2/7) (Z^(7/3)) ((e^2)/b) (-χ' (0))
The above numerical results agree very well with this analyic result.  On the other hand, the analytic result clearly shows (Subscript[E, e A]) + (Subscript[E, e e]) = -2 (Subscript[E, K]), which can be understood using the virial theorem explained below.

Virial Theorem

This is not a part of the homework problem, but it is very useful to know.  In classical mechanics in a bound system, the kinetic and potential energies averaged over a long time-interval are related in a simple fashion.  

As an example, take a single particle with the Hamiltonian H = ((p^2)/(2  m)) + λ (x^(2  n)) where λ is an arbitrary constant.  The Hamilton equation of motion is
(d/(d t)) 
          x = ∂ (H/∂) p = (p/m),
(d/(d t)) 
          p = -∂ (H/∂) x = -2 
            n λ (x^(2  n - 1)).

We first show that the average of any total time derivative (d/(d t)) A(x, p) over a long interval vanishes, 
\(TraditionalForm\`⟨\(d\/\(d\ t\)\) 
            A⟩ = \(lim\_\(T -> ∞\)\(1\/T\) \(∫\_0\%T A\ d\ t\) = 
          lim\_\(T -> ∞\)\(1\/T\) \((A(T) - A(0))\)\)\).
Because of the assumption that the system is bound, the difference A(T) - A(0) must be finite, and hence the limit goes to zero.

Now we apply the statement to A(x, p) = x p.  Then,
\(TraditionalForm\`0 = \(⟨\(d\/\(d\ t\)\) \((x\ p)\)⟩ = \(⟨\(\(d\ x\)\/\(d\ t\)\) 
                  p⟩ + ⟨x\ \(d\ p\)\/\(d\ t\)⟩ = ⟨\(p\/m\) 
                  p⟩ + ⟨x(\(-2\) 
                    n\ λ\ x\^\(2  n - 1\))⟩\)\)\)
\(TraditionalForm\`\(\(\(=\)\(⟨p\^2\/m⟩ - 
            2  n ⟨λ\ x\^\(2  n\)⟩ = 
          2 ⟨K⟩ - 
            2  n ⟨V⟩\)\)\(,\)\)\)
where K = ((p^2)/(2  m)) is the kinetic energy and V = λ (x^(2  n)) the potential energy.  Therefore we find ⟨K⟩ = 
        n ⟨V⟩.

In quantum mechanics, the time-average is replaced by the expectation value.  If the state is stationary, \(TraditionalForm\`\(\(\(\(\(\(i\ ℏ\ d\/\(d\ t\) | ψ\)\(⟩\)\) = H | ψ\)\(⟩\)\) = 
          E | ψ\)\(⟩\)\)\), and hence
\(TraditionalForm\`i\ ℏ \( 
              d\/\(d\ t\)\) ⟨ψ | 
              A | ψ⟩ = \(⟨ψ | \([A, H]\) | ψ⟩ = \(⟨ψ | \((A\ E - E\ A)\) | ψ⟩ = 0\)\)\).  The rest of the discussion is the same, 
\(TraditionalForm\`0 = \(\(d\/\(d\ t\)\) ⟨ψ | 
                x\ p | ψ⟩ = \(⟨ψ | \(p\/m\) p - 
                  x\ ∂\ V\/\(d\ x\) | ψ⟩ = 
            2 ⟨K⟩ - 
              2  n ⟨V⟩\)\)\).

For the single-electron atom, the application is straight-forward. 
\(TraditionalForm\`H = \(p\& -> \)\^2\/\(2  m\) - \(Z\ e\^2\)\/r\)
\(TraditionalForm\`\(d\ \(x\& -> \)\)\/\(d\ t\) = \(p\& -> \)\/m\),
\(TraditionalForm\`\(d\ \(p\& -> \)\)\/\(d\ t\) = \(-\(\(Z\ e\^2\)\/r\^3\)\) \(x\& -> \)\).
Therefore, 
\(TraditionalForm\`0 = \(\(d\/\(d\ t\)\) ⟨\(x\& -> \)·\(p\& -> \)⟩ = \(⟨\(p\& -> \)\/m·\(p\& -> \) - \(x\& -> \)·\(Z\ e\^2\)\/r\^3 \( 
                      x\& -> \)⟩ = 
            2 ⟨K⟩ + ⟨V⟩\)\)\).
We find \(TraditionalForm\`⟨V⟩ = \(-2\) ⟨K⟩\), which can be explicitly verified for any bound-state solution of the hydrogen-like atom.

It can be further generalized to the multi-electron case. 
\(TraditionalForm\`H = ∑\_i\((\(p\& -> \)\_i\^2\/\(2  m\) - \(Z\ e\^2\)\/r\_i)\) + \(1\/2\) \(\(∑\^'\)\_\(i, j\)e\^2\/r\_\(i\ j\)\)\),
where \(TraditionalForm\`\(∑\^'\)\) excludes i = j.  Then,
\(TraditionalForm\`\(d\ \(x\& -> \)\_k\)\/\(d\ t\) = \(p\& -> \)\_k\/m\)
\(TraditionalForm\`\(d\ \(p\& -> \)\_k\)\/\(d\ t\) = \(-\(\(Z\ e\^2\)\/r\_k\^3\)\) \(x\& -> \)\_k + ∑\_\(i ≠ k\)\(e\^2\/r\_\(i\ k\)\^3\) \((\(x\& -> \)\_k - \(x\& -> \)\_i)\)\)
Note that k is not summed over in the second term, and a factor of a half disappeared, because
\(TraditionalForm\`\(∂\/∂\ \(x\& -> \)\_k\) \(\(∑\^'\)\_\(i, j\)e\^2\/r\_\(i\ j\)\) = \(∑\^'\)\_\(i, j\)\(e\^2\/r\_\(i\ j\)\^3\) \((\((\(x\& -> \)\_k - \(x\& -> \)\_j)\) δ\_\(i\ k\) + \((\(x\& -> \)\_k - \(x\& -> \)\_i)\) δ\_\(j\ k\))\)\)
\(TraditionalForm\`\(\(=\)\(2 \(∑\_\(i ≠ k\)\(e\^2\/r\_\(i\ k\)\^3\) \((\(x\& -> \)\_k - \(x\& -> \)\_i)\)\)\)\)\)
Therefore,
\(TraditionalForm\`∑\_k\( p\& -> \)\_k·\(d\ \(x\& -> \)\_k\)\/\(d\ t\) = 2  K\),
\(TraditionalForm\`∑\_k\( x\& -> \)\_k·\(d\ \(p\& -> \)\_k\)\/\(d\ t\) = \(∑\_k\((\(-\(\(Z\ e\^2\)\/r\_k\)\))\) + ∑\_k\( x\& -> \)\_k·\(∑\_\(i ≠ k\)\(e\^2\/r\_\(i\ k\)\^3\) \((\(x\& -> \)\_k - \(x\& -> \)\_i)\)\) = V\)\)
A trick was used in the second term. It is summed over both i and k, and hence we can add another term with i and k interchanged, 
\(TraditionalForm\`∑\_k\( x\& -> \)\_k·\(∑\_\(i ≠ k\)\(e\^2\/r\_\(i\ k\)\^3\) \((\(x\& -> \)\_k - \(x\& -> \)\_i)\)\)\)
\(TraditionalForm\`\(\(=\)\(\(∑\^'\)\_\(i, k\)\(e\^2\/r\_\(i\ k\)\^3\) \(x\& -> \)\_k·\((\(x\& -> \)\_k - \(x\& -> \)\_i)\)\)\)\)
\(TraditionalForm\`\(\(=\)\(\(1\/2\) \(\(∑\^'\)\_\(i, k\)\(e\^2\/r\_\(i\ k\)\^3\) \((\(x\& -> \)\_k·\((\(x\& -> \)\_k - \(x\& -> \)\_i)\) + \(x\& -> \)\_i·\((\(x\& -> \)\_i - \(x\& -> \)\_k)\))\)\)\)\)\)
\(TraditionalForm\`\(\(=\)\(\(1\/2\) \(\(∑\^'\)\_\(i, k\)\(e\^2\/r\_\(i\ k\)\^3\) \((\(x\& -> \)\_i - \(x\& -> \)\_k)\)·\((\(x\& -> \)\_i - \(x\& -> \)\_k)\)\)\)\)\)
\(TraditionalForm\`\(\(=\)\(\(1\/2\) \(\(∑\^'\)\_\(i, j\)e\^2\/r\_\(i\ j\)\)\)\)\)
Therefore,
\(TraditionalForm\`0 = \(∑\_k⟨\(p\& -> \)\_k·\(d\ \(x\& -> \)\_k\)\/\(d\ t\) + \(x\& -> \)\_k·\(d\ \(p\& -> \)\_k\)\/\(d\ t\)⟩ = 
          2 ⟨K⟩ + ⟨V⟩\)\),
the same result as in the single-electron case.

comparison to data

I could not find the experimental values for large atoms (if you do, let me know!), but found numbers up to Z ≤ 8 in J. F. Barnes and R. D. Cowan, Phys. Rev. 132, 236-240 (1963).  Here is the comparison:

<< Graphics`Graphics`

LogPlot[20.92 (Z^(7/3)), {Z, 1, 8}]

LogListPlot[{{1, 13.60}, {2, 78.98}, {3, 203.4}, {4, 399.0}, {5, 
          670.8}, {6, 1030}, {7, 1486}, {8, 2043}}, 
      PlotStyle -> PointSize[0.02]]

Show[%, %%]

Given how crude the model is, and the comparison is made for small Z where the Fermi gas approximation is supposed to be not good, I find this agreement remarkably good.

(d)

We want to solve the equation Z - 1=(∫_0^R)d\(r\&_\) ρ(r) and get a plot of R(Z). 
(∫_0^R)d(Overscript[r,  ->])  ρ(r)=4π(∫_0^X)(((b (Z^(-1/3))))^3)dx (x^2) ((Z^2)/(4  π (b^3)))(((((χ(x))/x)))^(3/2))=Z(∫_0^X)dx ((χ(x))^(3/
                        2))(Sqrt[x])
Using the differential equation,
(∫_0^R)d(Overscript[r,  ->])  ρ(r)=Z(∫_0^X)dx χ'' x=Z{(Subscript[[χ' x], 0]^X)-(∫_0^X) χ' d x}
=Z {χ' (X) X - χ(X) + χ(0)} = 
            Z {χ' (X) X - χ(X) + 1}
Setting the r.h.s equal to Z - 1 we get
χ' (X) X - χ(X) + 1 = 1 - (1/Z),
or
Z = (((χ(X) - χ' (X) X))^(-1)).

size = 
      ParametricPlot[
        Evaluate[{(((χ[x] - χ'[x] x))^(-1)), (((χ[
                            x] - χ'[x] x))^(1/3)) b x /. 
                b -> 0.5} /. sol], {x, 0, 30}, 
        AxesLabel -> {Z, R[(Overscript[A, o])]}, 
        PlotRange -> {{0, 100}, {0, 3}}]

This plot demonstrates the weak dependence of R on Z. As Z changes form 10 to 80 R changes only by ~20%.

comparison to data

There are many definitions of atomic size, van der Waals radii, covalent radii, ionic radii, metallic radii.  Here we use the covalent radii I got from http://www.ccdc.cam.ac.uk/products/csd/radii/

Rcovalent = {{1.0000, 0.23000}, {2.0000, 1.5000}, {3.0000, 
          0.68000}, {4.0000, 0.35000}, {5.0000, 0.83000}, {6.0000, 
          0.68000}, {7.0000, 0.68000}, {8.0000, 0.68000}, {9.0000, 
          0.64000}, {10.000, 1.5000}, {11.000, 0.97000}, {12.000, 
          1.1000}, {13.000, 1.3500}, {14.000, 1.2000}, {15.000, 
          1.0500}, {16.000, 1.0200}, {17.000, 0.99000}, {18.000, 
          1.5100}, {19.000, 1.3300}, {20.000, 0.99000}, {21.000, 
          1.4400}, {22.000, 1.4700}, {23.000, 1.3300}, {24.000, 
          1.3500}, {25.000, 1.3500}, {26.000, 1.3400}, {27.000, 
          1.3300}, {28.000, 1.5000}, {29.000, 1.5200}, {30.000, 
          1.4500}, {31.000, 1.2200}, {32.000, 1.1700}, {33.000, 
          1.2100}, {34.000, 1.2200}, {35.000, 1.2100}, {36.000, 
          1.5000}, {37.000, 1.4700}, {38.000, 1.1200}, {39.000, 
          1.7800}, {40.000, 1.5600}, {41.000, 1.4800}, {42.000, 
          1.4700}, {43.000, 1.3500}, {44.000, 1.4000}, {45.000, 
          1.4500}, {46.000, 1.5000}, {47.000, 1.5900}, {48.000, 
          1.6900}, {49.000, 1.6300}, {50.000, 1.4600}, {51.000, 
          1.4600}, {52.000, 1.4700}, {53.000, 1.4000}, {54.000, 
          1.5000}, {55.000, 1.6700}, {56.000, 1.3400}, {57.000, 
          1.8700}, {58.000, 1.8300}, {59.000, 1.8200}, {60.000, 
          1.8100}, {61.000, 1.8000}, {62.000, 1.8000}, {63.000, 
          1.9900}, {64.000, 1.7900}, {65.000, 1.7600}, {66.000, 
          1.7500}, {67.000, 1.7400}, {68.000, 1.7300}, {69.000, 
          1.7200}, {70.000, 1.9400}, {71.000, 1.7200}, {72.000, 
          1.5700}, {73.000, 1.4300}, {74.000, 1.3700}, {75.000, 
          1.3500}, {76.000, 1.3700}, {77.000, 1.3200}, {78.000, 
          1.5000}, {79.000, 1.5000}, {80.000, 1.7000}, {81.000, 
          1.5500}, {82.000, 1.5400}, {83.000, 1.5400}, {84.000, 
          1.6800}, {85.000, 1.2100}, {86.000, 1.5000}, {87.000, 
          1.5000}, {88.000, 1.9000}, {89.000, 1.8800}, {90.000, 
          1.7900}, {91.000, 1.6100}, {92.000, 1.5800}, {93.000, 
          1.5500}, {94.000, 1.5300}, {95.000, 1.5100}, {96.000, 
          0.99000}, {97.000, 1.5400}, {98.000, 1.8300}, {99.000, 
          1.5000}, {100.00, 1.5000}, {101.00, 1.5000}, {102.00, 
          1.5000}, {103.00, 1.5000}, {104.00, 1.5000}, {105.00, 
          1.5000}, {106.00, 1.5000}, {107.00, 1.5000}, {108.00, 
          1.5000}, {109.00, 1.5000}, {110.00, 1.5000}}

data = 
      ListPlot[Rcovalent, PlotStyle -> PointSize[0.01]]

Show[size, data]

They are not expected to be right on because of the arbitrariness in definitions.  However, the qualitative feature of the rise and near saturation is well reproduced. If I use a "fudge" factor of about 0.6 on the Thomas-Fermi calculations,

size2 = 
      ParametricPlot[
        Evaluate[{(((χ[x] - χ'[x] x))^(-1)), 
              0.6 (((χ[x] - χ'[x] x))^(1/3)) b x /. 
                b -> 0.5} /. sol], {x, 0, 30}, 
        AxesLabel -> {Z, R[(Overscript[A, o])]}, 
        PlotRange -> {{0, 100}, {0, 2}}]

Show[size2, data]

Not bad at all.

In summary, the Thomas-Fermi model describes the "bulk" properties of atoms reasonably well, while it cannot keep track of the periodicity in chemical properties because we never use the information of electronic configurations.

I see many plots that show the atomic size.  As far as I understand, these plots that show the periodicity beautifully are calculated, not measured.  Here is what I got from http://web.mit.edu/3.091/www/pt/pert1.html which does not give a definition

Ratomic = {{1, 2.08}, {3, 1.55}, {4, 1.12}, {5, 0.98}, {6, 0.91}, {7, 
          0.92}, {8, 0.65}, {9, 0.57}, {10, 0.51}, {11, 1.9}, {12, 1.6}, {13, 
          1.43}, {14, 1.32}, {15, 1.28}, {16, 1.27}, {17, 0.97}, {18, 
          0.88}, {19, 2.35}, {20, 1.97}, {21, 1.62}, {22, 1.45}, {23, 
          1.34}, {24, 1.3}, {25, 1.35}, {26, 1.26}, {27, 1.25}, {28, 
          1.24}, {29, 1.28}, {30, 1.38}, {31, 1.41}, {32, 1.37}, {33, 
          1.39}, {34, 1.4}, {35, 1.12}, {36, 1.03}, {37, 2.48}, {38, 
          2.15}, {39, 1.78}, {40, 1.6}, {41, 1.46}, {42, 1.39}, {43, 
          1.36}, {44, 1.34}, {45, 1.34}, {46, 1.37}, {47, 1.44}, {48, 
          1.71}, {49, 1.66}, {50, 1.62}, {51, 1.59}, {52, 1.42}, {53, 
          1.32}, {54, 1.24}, {55, 2.67}, {56, 2.22}, {57, 1.38}, {58, 
          1.81}, {59, 1.82}, {60, 1.82}, {62, 1.81}, {63, 1.99}, {64, 
          1.8}, {65, 1.8}, {66, 1.8}, {67, 1.79}, {68, 1.78}, {69, 
          1.77}, {70, 1.94}, {71, 1.75}, {72, 1.67}, {73, 1.49}, {74, 
          1.41}, {75, 1.37}, {76, 1.35}, {77, 1.36}, {78, 1.39}, {79, 
          1.46}, {80, 1.6}, {81, 1.71}, {82, 1.75}, {83, 1.7}, {84, 
          1.67}, {85, 1.45}, {86, 1.34}, {87, 2.7}, {88, 2.33}, {89, 
          1.88}, {90, 1.8}, {91, 1.61}, {92, 1.38}, {93, 1.3}, {94, 
          1.51}, {95, 1.84}}

data = 
      ListPlot[Ratomic, PlotStyle -> PointSize[0.01], 
        PlotJoined -> True, PlotRange -> {0, 3}]

Show[size, data]

Show[size2, data]