Midterm (221b)

1)

a)

Following the "technical remark" in the lecture notes on the Thomas-Fermi equation, we rewrite χ[x] = (1 + (4/3) (x^(3/2))) y[x] and solve for y.  The initial condition y'[0]=-1.588071 is found by trial and error, demanding that the function χ decay smoothly to zero at infinity, which is the condition for a neutral atom.  In practice, the function is senstive at least to the 5th digit of y'[0], so finding this by trial and error is a real bear.  I let you get away with using the value y'[0]=-1.588071 given in the lecture notes but took off a couple of points if you didn't discuss where it comes from.

tf = 
      NDSolve[{(1 + (4/3) (x^(3/2))) y''[x] + 
              4  (x^(1/2)) y'[x] + 
              If[x == 0, 
                0, (1/(x^(1/2))) 
                  y[x] (1 - (((1 + (4/3) (x^(3/2))))^(3/2)) 
                        y[x]^(1/2))] == 0,  
          y[0] == 1,  y'[0] == -1.588071}, 
        y, {x, 0, 30}]

χ[x_] := 
      Evaluate[(1 + (4/3) (x^(3/2))) y[x] /. tf]

Plot[χ[x],  {x, 0, 20}]

b)

The point here is the heuristic argument given in the lecture notes "Many Body Problems II", p. 6:  The Thomas-Fermi function χ describes the screening of the nuclear charge.  As expected, the electrons which stay in the outer part of the atom (electrons with large l) see smaller χ and hence a more-screened nuclear charge, so their binding-energy is smaller (less negative) than electrons which probe the interior of the atom (electrons with small l).  Quantitative support would follow from solving the single-electron  radial wave equation in the effective potential energy \(TraditionalForm\`\(V\_eff\)(
          r) = \(-\(Ze\^2\/r\)\) \(χ(
              r)\) + \(\(\^2\) \(l(l + 1)\)\)\/\(2  mr\^2\)\) to find the energies for different values of l.  If necessarly, the radial Schrodinger equation can be integrated numerically, using trial and error to find energy eigenvalues by requiring that the bound-state  wavefunction vanish at infinity.

A similar question which is easier to analyze quantitatively is what different l values can appear for different values of Z:  For a single electron to be bound in the effective potential energy \(TraditionalForm\`\(V\_eff\)(
          r) = \(-\(Ze\^2\/r\)\) \(χ(
              r)\) + \(\(\^2\) \(l(l + 1)\)\)\/\(2  mr\^2\)\), its energy E must satisfy (Subscript[V, eff])(r) < E < 0 for some range in r.  This means that for some x,  (Subscript[V, eff])(x)<0 , or 2 
           (Z^(2/3)) b x χ(x) > (Subscript[a, B]) l(l + 1), where  b  =  (1/2) (((((3  π)/4)))^(2/3)) 
          (Subscript[a, B]) and r =  (Z^(-1/
                3)) b x  as usual.  For a given l, Z must be sufficiently large that this inequality is satisfied for some x.  From the graph...

Plot[x χ[x],  {x, 0, 20}]

...you can see that x χ(x) has a maximum value of about .49.  Plugging in this value and the expression for b, the inequality becomes roughly   .87  (Z^(2/3))  > l(l + 1) .  For l = 0, Z can be anything, but for l=1, Z must be at least (((2/ .87))^(3/2)) -> 4, etc.  Following this pattern gives the following table, which tells us something about how transition elements appear in the periodic table:

l		(Subscript[Z, min])(observed)          (Subscript[Z, min])(Thomas - Fermi)	
0		1			1
1		5			4
2		21			19
3		58			52

I was talking with Prof. Commins earlier today, and he says that with an improved version of the Thomas-Fermi  model that takes account of electron-electron  correlations, you can find the better values (Subscript[Z, min])= .2, 5.4, 21, 58.

c)

The total binding energy is a sum of kinetic and potential energies.  The kinetic energy of a degenerate fermi gas has known expression in terms of particle number density ρ(r), namely \(TraditionalForm\`E(
          r) = \ \(\(3\^\(5/3\)\) \(π\^\(4/3\)\) \(\^2\) \(ρ(r)\)\^\(5/3\)\)\/\(\(10\)\(\ \)\(m\)\(\ \)\)\).  The potential energy can be calculated in two pieces, an electron-electron  interaction piece, \(TraditionalForm\`\(1\/2\) \(∫\(\(\(\(e\^2\) ρ\ \((r)\) \(ρ(r')\)\)\/\(\(|\)\(r\&⇀ - 
                    r\&⇀'\)\(|\)\)\) \(d\^3\) \(r\&⇀\) \(d\^3\) \(r\&⇀'\)\(\ \)\)\)\), and an electron-nucleus  piece, -∫e ρ(r) (Ze/r).   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.

The kinetic energy density of a totally degenerate fermi gas is \(TraditionalForm\`E(
          r) = \ \(\(3\^\(5/3\)\) \(π\^\(4/3\)\) \(\^2\) \(ρ(r)\)\^\(5/3\)\)\/\(\(10\)\(\ \)\(m\)\(\ \)\)\), and we know the spatial density as a function of r in terms of the Thomas-Fermi function, \(TraditionalForm\`ρ(
          r) = \(\((2  me)\)\^\(3/2\)\/\(3 \( π\^2\) \^3\)\) \((\(Ze\/r\) \(χ(r)\))\)\^\(3/2\)\).  Integrating E(r) over 4  π (r^2) d r=4  π  ((b^3)/Z) (x^2) d x gives the total kinetic energy in the atom.  Putting everything together, \(E\_kinetic = \(e\^2\/a\_B\) \(∫\(2\/\(5  π\)\) \(\((\(32  b\)\/a\_B)\)\^\(1/2\)\) Z\^\(7/3\)\)\)((χ^(5/2))/(x^(1/2)))dx .  Extracting the prefactor ((e^2)/(Subscript[a, B])) (Z^(7/3)), the kinetic energy is

kinetic =  
      NIntegrate[
        Evaluate[(2/(5  π)) (((32 (1/2) (((((3  π)/4)))^(2/3))))^(1/2)) (x^(-1/
                  2)) (((χ[x]))^(5/2))],  {x, 0, 30}]

For future use, we can rewrite the integral \(TraditionalForm\`∫\ χ\^\(5/2\)\/x\^\(1/2\) = \(\(-5\) \(∫χ\^\(3/2\)\ χ'\ x\^\(1/2\)\) = \(-5\) \(∫χ'\ χ'' x\)\)\), 
the first equality from integrating by parts (the boundary term vanishing by the b.c. on χ at infinity), the second equality from the Thomas-Fermi  equation χ'' = ((χ^(3/2))/(x^(1/2))).  Integrating the r.h.s by parts and rearranging the resulting equation shows that \(TraditionalForm\`∫\ χ\^\(5/2\)\/x\^\(1/2\) = \ \(5\/2\) \(\(∫\)\(\ \)\(χ'\^2\)\(\ \)\)\), so that 
(Subscript[E, kinetic])= \(TraditionalForm\`\(e\^2\/a\_B\) \(\(\(Z\^\(7/
                    3\)\)(\(32  b\)\/a\_B)\)\^\(1/2\)\) 
        1\/π\ \(∫\ \(χ'\^2\) dx\)\) .

The first part of the potential energy is (Subscript[E, elec - elec])= \(TraditionalForm\`\(1\/2\) \(∫\(\(\(e\^2\) ρ\ \((r)\) \(ρ(r')\)\)\/\(\(|\)\(r\&⇀ - 
                    r\&⇀'\)\(|\)\)\) \(d\^3\) \(r\&⇀\) \(d\^3\) r\&⇀'\)\).  Plugging in the Thomas-Fermi  functions where appropriate, writing (Overscript[r, ⇀])  =  
        (Overscript[x, ⇀]) (Z^(-1/3)) b, and consolidating constants,
 (Subscript[E, elec - elec])= \(\(e\^2\/a\_B\) \(\(Z\^\(7/3\)\)(\(32  b\)\/a\_B)\)\^\(1/2\)\)((((b/(Subscript[a, B]))))^(3/
                2))((Sqrt[2])/(18  (π^4))) \(∫\ \(\(\(\(χ(x)\)\^\(3/
                          2\)\) \(χ(x')\)\^\(3/2\)\)\/\((xx')\)\^\(3/2\)\) \(1\/\@\(x\^2 + x'\^2 - 2  xx' cosθ\_12\)\) \(d\^3\) \(x\&⇀\) \(d\^3\) x\&⇀'\).  
This can be computed numerically, but the final result makes it worth an exact computation:  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:

\(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, plugging in for ((((b/(Subscript[a, B]))))^(3/2)), and comparing with the 2nd form we derived for (Subscript[E, kinetic]), you can check that we have found exactly (Subscript[E, elec - elec]) =  (1/3) (Subscript[E, kinetic]).

The second part of the potential energy is (Subscript[E, elec - nuc]) = -∫e ρ(
                r) (Ze/r).  Plugging in and using the Thomas-Fermi  equation, this becomes
(Subscript[E, elec - nuc])= -(2/(3  π)) ((e^2)/(Subscript[a, B])) (((Z^(7/
                    3))(((32  b)/(Subscript[a, B]))))^(1/
              2)) ∫χ'' (x) dx = (2/(3  π)) ((e^2)/(Subscript[a, B])) (((Z^(7/
                    3))(((32  b)/(Subscript[a, B]))))^(1/2)) χ' (0).
We can of course plug in numerically for our value of χ'(0), again extracting a prefactor ((e^2)/(Subscript[a, B])) (Z^(7/3)):

elecnuc =  (2/(3  π)) (((32 (1/2) (((((3  π)/4)))^(2/3))))^(1/2)) D[χ[x], x] /. x -> 0

Numerically, (Subscript[E, elec - nuc])= -(7/3) (Subscript[E, kinetic]):

elecnuc/kinetic

But it would be nice to determine (Subscript[E, elec - nuc]) analytically to go with our computation of (Subscript[E, elec - elec]).  Further,  we could use an analytical computation as a self-consistency  check on our value of χ'(0).  One way follows:  Compute
(Subscript[E, wrong - potential]) = -(1/2) ∫e ρ(
              r) ϕ(r)  (d^3) r
 = \(TraditionalForm\`\(-2\) π\ e \(\((2  me)\)\^\(3/2\)\/\(3 \( π\^2\) \^3\)\) \(∫\ \(\((\(Ze\/r\) \(χ(r)\))\)\^\(3/
                  2\)\) \(Ze\/r\) \(χ(r)\)\ \(r\^2\) dr\)\) 
 = \(TraditionalForm\`\(-\(e\^2\/a\_B\)\) \(∫\ \(1\/\(3  π\)\) \(\((\(32  b\)\/a\_B)\)\^\(1/2\)\) \(Z\^\(7/
                  3\)\) \(χ\^\(5/2\)\/x\^\(1/2\)\) dx\)\) 
 = -(5/6) (Subscript[E, kinetic]).
We know from the reasoning at top that (Subscript[E, wrong - potential])= (Subscript[E, elec - elec]) + (1/2) (Subscript[E, elec - nuc]).  Our previous computation of (Subscript[E, elec - elec]) =  (1/3) (Subscript[E, kinetic]), along with the result (Subscript[E, wrong - potential]) = -(5/6) 
          (Subscript[E, kinetic]), then implies (Subscript[E, elec - nuc]) = -(7/3) (Subscript[E, kinetic]) exactly.

The grand total energy is (Subscript[E, total])= (Subscript[E, kinetic]) + (Subscript[E, elec - elec]) + 
        (Subscript[E, elec - nuc])= -(Subscript[E, kinetic])->  - .769 (Z^(7/3)) in atomic units.  The (Z^(7/3)) dependence is a definite prediction of the Thomas-Fermi  model (in contrast, say, to the E  ∝  A behavior of nuclei).  The Thomas-Fermi  model lends itself well to  Z-scaling  arguments like this and the one discussed in part (b).

In electronvolts the total binding energy is (Subscript[E, total])≈ -20.9 (Z^(7/3)) eV.  Several of you compared your results with the experimental values.  I wasn't smart enough to copy down your data, but I have one number from Elvis Dieguez via Prof. Murayama:  E(Copper) ≈  45 038.42 eV.  The Thomas-Fermi model gives

.769 (((29))^(7/3)) 27.21

which is unstunning.  Again referring to Prof. Commins, the numerical factor .769 should only be trusted within 10% (or worse, it would seem in our case), but the (Z^(7/3)) behavior matches experiment quite well.

The total potential energy is (Subscript[E, elec - elec]) + (Subscript[E, elec - nuc])= -2 (Subscript[E, kinetic]), which looks exactly like the virial theorem result for a particle in a 1/r potential!

Recall that the virial theorem for a particle moving in a potential V with kinetic energy K says  \(TraditionalForm\`2  K\&_\)= \(TraditionalForm\`\(r\&⇀·∇V\)\&_\), where the bar denotes time average.  For the 1/r potential, this becomes \(TraditionalForm\`2  K\&_\)=-\(TraditionalForm\`V\&_\).  In quantum mechanics time averages become expectation values, 2 ⟨ 
            K ⟩ =  ⟨ 
          (Overscript[r, ⇀])·∇V ⟩.  What's wierd is that in our problem, we have many electrons with a collective kinetic energy and a complicated potential -((Ze^2)/r) χ(r), so it's not obvious how or if these virial theorems apply.

There is a many-body virial  theorem from statistical mechanics, and it turns out  to be nothing but the ensemble average of the old single-particle  results from mechanics.  Here's how it goes (ref. Pathria Statistical Mechanics §3.7):   Let (Subscript[q, i]) be any of the 3N spatial coordinates for an N particle system, let dω be the element of phase space density, and consider the ensemble average (note new meaning for ⟨ ⟩) of (Subscript[q, i]) ∂(H/∂)(Subscript[q, j]) in the canonical ensemble:

\(TraditionalForm\`⟨\ 
          q\_i\ ∂H\/∂q\_j⟩\  = \ \(∫q\_i\ \(∂H\/∂q\_j\) \(e\^\(-βH\)\) dω\)\/\(∫\(e\^\(-βH\)\) dω\)\)

	= -(1/β)\(TraditionalForm\`\(∫q\_i\ \(∂\/∂q\_j\) \(e\^\(-βH\)\) dω\)\/\(∫\(e\^\(-βH\)\) dω\)\)
	
	= (1/β)(Subscript[δ, ij]),

the last equality following by integration by parts.  Now take i=j and sum over all i.  Using 1/β = kT, 

\(TraditionalForm\`\(\(⟨∑\ 
          q\_i\ ∂H\/∂q\_i⟩\)\(\ \)\)\)= 3NkT.

For our problem, the (Subscript[q, i]) dependence in the Hamiltonian comes only from the interparticle 
\(TraditionalForm\`1\/\(\(|\)\(r\&⇀\_i\  - \ r\&⇀\_j\)\(|\)\)\)  Coulomb potentials, so ⟨∑ 
          (Subscript[q, i]) ∂(H/∂)(Subscript[q, i])⟩ is really just -⟨ 
          V ⟩, the ensemble average of  interaction energies summed over all pairs of charges.  The right hand side, 3NkT, is twice the average kinetic energy.  In other words, while the total potential is the complicated Thomas-Fermi  function, the virial theorem only cares about the underlying interactions which are the usual 1/r potentials of mechanics.  Averaged over the ensemble, -⟨ V ⟩ = 
        2 ⟨ T ⟩.  Still it is rather mysterious and wonderful to me that the energies we computed above with electrostatics obey this statistical virial theorem.

On further reflection, the ordinary mechanical virial theorem probably holds in a similar way, time averaged and summed over all particles.  It's the fundamental 1/r interactions that govern the virial relation.  A few of you invoked the virial theorem in doing this problem, and if what I've said here was clear to you before, I owe you a few points on your midterms.

d)

The integral is basically calculated on p. 9 of the lecture notes  "Many Body Problems II," though for us χ((Subscript[x, 0])) ≠  0 (with R = (Z^(-1/3)) b (Subscript[x, 0])).  We have Z - 1 = Z((Subscript[x, 0]) χ' ((Subscript[x, 0])) - χ((Subscript[x, 0])) + 
          1), or Z = -\(TraditionalForm\`1\/\(\(x\_0\) χ' \((x\_0)\) - χ(x\_0)\)\).  This defines Z implicitly as a function of (Subscript[x, 0]), and we already know R =  .47  (Z^(-1/3)) (Subscript[x, 0]) (b ≈   .47 angstrom).  This calls for a parametric plot:

Z[t_] := -(1/(t χ'[t] - χ[t]));  
    R[t_] :=  .47  Z[t]^(-1/3) t

ParametricPlot[{Z[t][[1]], R[t][[1]]}, {t, 0, 
        20}]