Physics 221B, Solution to Midterm Problem #1:
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 it extremely sensitive to χ '(0), so we need to take all of the digits suggested in the notes. In fact, I had to add a few digits to χ '(0) to make it well-behaved all the way upto x = 100.
After getting used to dealing with Mathematica's InterpolatingFunction[ ] and plotting we get the result familiar from the notes.

f[x_] := 
      If[x == 0, 
        0, y[x]/(x^(1/2)) (1 - (((1 + 4/3*(x^(3/2))))^(3/2))*
                y[x]^(1/2))]

sol = 
      NDSolve[{(1 + 4/3*(x^(3/2))) y''[x] + 
              4*(x^(1/2)) y'[x] + f[x] == 
            0, y'[0] == -1.5880708425,  y[0] == 1}, 
        y, {x, 0, 100}]
χ[x_] := (1 + 4/3*(x^(3/2))) y[x]

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

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.
A more tentative answer (that I will not do) can be to (numerically) solve the Schrodinger equation for an electron in the potential V(r)=-(Ze/r) χ(r)+(((()^2)l(l + 1))/((r^2)))for different values of l.

c) The Total Binding Energy

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

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)) 
          (Subscript[a, B]) we get
(Subscript[E, k]) = (3/5) (((Z^(7/3)) (e^2))/b)∫((χ^(5/2))/(Sqrt[x]))dx.
Doing the integral with Mathematica gives

3/5 NIntegrate[
        Evaluate[χ[x]^(5/2)/(Sqrt[x]) /. sol], {x, 0, 100}]

%/((1/2) (((((3  π)/4)))^(2/3)))

We get (Subscript[E, k]) = 0.680602 (( (Z^(7/3)) (e^2))/b) = 
                0.768745(Z^(7/3))((e^2)/(Subscript[a, B])).

This result is certainly enough, but it will be useful later to play with the integral above. Note that \(TraditionalForm\`∫\ χ\^\(5/2\)\/x\^\(1/2\) = \(\([2  χ\^\(5/2\)\ x\^\(1/2\)]\)\_0\%∞ - 
            5 \(∫χ\^\(3/2\)\ χ'\ x\^\(1/2\)\) = \(\(0 - 
              5 \(∫χ\^\(3/2\)\ χ'\ x\^\(1/2\)\)\)\(=\)\)\)\), 
where the last equality from the Thomas-Fermi  equation χ'' = ((χ^(3/2))/(x^(1/2))).  Integrating the r.h.s by parts once more gives 
\(TraditionalForm\`\(-5\) \(∫χ'\ χ'' 
              x\) = \(-\([\(5\/2\) \(\((χ')\)\^2\) x]\)\_0\%∞\) + \(5\/2\) \(∫\ χ'\^2\)\)
and the boundary term drops once more. The resulting equation shows that \(TraditionalForm\`∫\ χ\^\(5/2\)\/x\^\(1/2\) = \ \(5\/2\) \(\(∫\)\(\ \)\(χ'\^2\)\(\ \)\)\), so that the kinetic energy is
(Subscript[E, k]) = (3/5) (((Z^(7/3)) (e^2))/b)\(TraditionalForm\`∫\ χ'\^2\).

Potential Energy:

The potential energy should be separated into the electron-nuclear interaction and the electron-electron repulsion. Treating the electron cloud as a charge density ρ=\(TraditionalForm\`\(Z\^2\/\(4  π\ b\^3\)\) \((χ\/x)\)\^\(3/2\)\) , (for positive χ) we can write 
(Subscript[E, n - e])=-∫d(Overscript[r, ⇀]) eρ(r)(Ze/r)
(Subscript[E, e - e])=(1/2) ∫  ((e ρ(r) e ρ(
                      r'))/(| (Overscript[r, ⇀])  -  
                      (Overscript[r, ⇀])' |))(d^3)\(r\&⇀\)(d^3)(Overscript[r, ⇀])'
where the factor of 1/2 in (Subscript[E, e - e]) is there to prevent double counting.
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 above.

Lets start with the electron-nuclear term. Changing variables in the integral above from r to x and writing ρ in terms of χ
(Subscript[E, n - e])=-(((Z^(7/3)) (e^2))/b)∫dx (((χ(x))^(3/2))/(Sqrt[x])).

NIntegrate[
      Evaluate[χ[x]^(3/2)/(Sqrt[x]) /. sol], {x, 0, 100}]

%/((1/2) (((((3  π)/4)))^(2/3)))

We get \(E\_\(n - e\) = \(\(-1.58805\) \(\(Z\^\(7/3\)\) e\^2\)\/b = \(-1.79371\)\)\)(Z^(7/3))((e^2)/(Subscript[a, B])).

Moving on to the electron-electron interaction, plugging Thomas-Fermi functions where appropriate and changing integration variables we get
(Subscript[E, e - e]) = ((e^2)/(32  (π^2) b)) 
                (Z^(7/3)) \(∫\ \(\(\(\(χ(x)\)\^\(3/
                          2\)\) \(χ(x')\)\^\(3/2\)\)\/\((xx')\)\^\(3/2\)\) \(1\/\@\(x\^2 + x'\^2 - 2  xx' cosθ\_12\)\) \(d\^3\) \(x\&⇀\) \(d\^3\) x\&⇀'\).
Here is one possible way to solve this integral. Choose the x' axis so that \(TraditionalForm\`\(d\^3\) x\&⇀' -> \ 
        2  π\ dcos\ \(θ\_12\) \(x'\^2\) dx'\), which makes the dcos(Subscript[θ, 12]) integral easy: 
  \(TraditionalForm\`∫\ \ dcosθ\_12\/\@\(x\^2 + x'\^2 - 2  xx' cosθ\_12\)\) = \(TraditionalForm\`\((\(-\@\((x - x')\)\^2\) + \@\((x + x')\)\^2)\)\/xx'\)= (2/x if x>x' or 2/x' if x<x').
  Because of this result the dx' integral must be split into two parts, ((∫_0)^x) and ((∫_x)^∞), but the latter piece can be rewritten identically to the former:  Integrating ((∫_0)^∞) 
          dx ((∫_x)^∞) dx' is the same as integrating ((∫_0)^∞) 
            dx' ((∫_0)^x)'dx , since they are equivalent ways to add up the shaded area in the triangle below:

\(Show[
      Graphics[Polygon[{{0, 0}, {10\^6, 10\^6}, {0, 10\^6}}], 
        Axes -> True, 
        Ticks -> {{{10\^6, "∞"}}, {{10\^6, "∞"}}}, AxesLabel -> {x, x'}]]\)

After the above adjustment, relabel variables in the second term to give exactly the first term again.  In all, this leaves 32  (π^2) ((∫_0)^∞) 
            dx ((∫_0)^x) 
                  dx' ((χ(x))^(3/2)) ((χ(x'))^(3/
                      2)) (Sqrt[(x/x)']).  Using the Thomas-Fermi equation and integrating by parts several times,
32  (π^2) ((∫_0)^∞) 
            dx ((∫_0)^x) 
                  dx' ((χ(x))^(3/2)) ((χ(x'))^(3/
                      2)) (Sqrt[(x/x)'])
= 32  (π^2) ((∫_0)^∞) 
            dx ((∫_0)^x) dx' χ'' (x) χ(
                  x') x'
= \(TraditionalForm\`32 \( π\^2\) \(∫\_0\^∞ dx\ [
            x' χ' \((x')\)\( | \_0\)\^x\(-\(∫\_0\^x 
                      dx'\ χ' \((x')\)\)\)]\)\)
= \(TraditionalForm\`32 \( π\^2\) \(∫\_0\^∞ 
              dx\ \ χ'' \((x)\) χ' \((x)\)\ x\)\  - \ χ'' \((x)\)\ [\(χ(x')\)\(\( | \_0\)\^x\)]\)
= \(TraditionalForm\`32 \( π\^2\) \(∫\_0\^∞ 
              dx\ \ χ'' \((x)\) χ' \((x)\)\ x\)\  - \ χ'' \((x)\) \(χ(x)\)\  + \ χ'' \((x)\)\)
= \(TraditionalForm\`\(\(32\) \(\( π\^2\)[\(x\/2\) χ' \((x)\)\^2\( | \_0\)\^∞\(-\(1\/2\)\) \(∫\_0\^∞ 
                    dx\ \ χ' \((x)\)\^2\) - \(χ(
                  x)\) χ' \((x)\)\( | \_0\)\^∞\(-χ'\) \((x)\)\( | \_0\)\^∞]\)\(\ \)\)\)
= 16  (π^2) ((∫_0)^∞) 
            dx  χ' (((x))^2).
Restoring the constants out front, we get (Subscript[E, e - e]) = (1/2) ((e^2)/b) (Z^(7/
                3)) ((∫_0)^∞) 
              dx  χ' (((x))^2). Recalling our result for the kinetic energy from above, we get exactly (Subscript[E, elec - elec]) =  (1/3) 
            (Subscript[E, kinetic]) = 0.226867  ((e^2)/b) (Z^(7/3)) = 
            0.25624828233699587`  ((e^2)/(Subscript[a, B])) (Z^(7/3)).

Total Binding Energy:

Putting everything together gives
(Subscript[E, p]) = (-1.79371 + 
                0.256248) ((e^2)/(Subscript[a, B])) 
            (Z^(7/3)) = -1.566843 ((e^2)/(Subscript[a, B])) 
              (Z^(7/3)) = -42.63 eV *(Z^(7/3))
(Subscript[E, k]) = 0.766745  ((e^2)/(Subscript[a, B])) (Z^(7/3)) = 
          20.91 eV *(Z^(7/3))
(Subscript[E, bind]) = 21.71  eV *(Z^(7/3))

Examining our result we notice that the virial theorem is very nearly satisfied, \(TraditionalForm\`\(\(|\)\(E\_p\)\(|\)\(\(/\)\(\(|\)\(E\_k\)\(|\)\)\)\) = 2.03\). We are lead to suspect that the virial theorem hold exactly and the small deviation from 2 is the result of our use of numerical methods. Doing a little extra work we can show that the virial theorem holds exactly.

d)

We want to solve the equation Z - 1=(∫_0^R)d\(r\&_\) ρ(r) and get a plot of R(Z).
(∫_0^R)d\(r\&_\) ρ(r)=4π(∫_0^X)((((bZ^(-1/3))))^3)dx (x^2) ((Z^2)/(4  π (b^3)))(((((χ(x))/x)))^(3/2))=Z(∫_0^X)dx ((χ(x))^(3/
                        2))(Sqrt[x])
Setting the r.h.s equal to Z - 1 we get

Z=(([1-(∫_0^X)dx (((((χ(x))/x)))^(3/2))])^(-1))
with((  R = (Z^(-1/3)) 
                      b X=[1-(∫_0^X)dx (((((χ(x))/x)))^(3/2))])^(1/3))b X

The equations above are in the form we want. Once we specify X we can easily find Z and R. All we need to do is let X vary over a range of values. We can easily do this with Mathematica's parametric plot function.

f[X_] := 
      NIntegrate[
        Evaluate[χ[x]^(3/2) (Sqrt[x]) /. sol[[1]]], {x, 0, 
          X}]

ParametricPlot[{(((1 - f[X]))^(-1)), (((1 - f[X]))^(1/3)) 
            b X /. b -> 0.5}, {X, 0, 25}, 
      AxesLabel -> {Z, R[(Overscript[A, o])]}]

This plot demonstrates the weak dependence of R on Z. As Z changes form 10 to 80 R changes only by ~20%. Another possible way to get this plot is using Eq. (39) in the lecture notes.