hw #6

Q1 G8.8

Only the term from the nucleus-electron coulomb interaction is modified by the substitution He → H- or Li+. Therefore, we can start working straight from equation 8.33, making the substitution $Z-2 \to Z-Z_o$, where $Z_o$ is the charge of the nucleus in question (+1 for hydrogen, +3 for lithium). Minimizing this function with respect to $Z$ we find: $$ Z_{\text{min}} = Z_o - 5/16 $$ So that the effective charge that minimizes the trial wavefunction is some constant shift from the full nuclear charge. This is consistent with our interpretation of the modification $Z_o \to Z_{\text{min}}$ resulting one electron screening the other. The resulting upper bounds on the ground state energies for these systems are obtained, as usual, by evaluating the expectation values of $\hat{H}$ at $Z_{\text{min}}$: $$ E_{\text{min}} = E_1 \left( 2 Z_o^2 - \frac{5}{4} \left( Z_o - 5/32 \right) \right) $$ Where $E_1 \approx -13.6\text{eV}$ is the ground state energy of the (neutral) hydrogen atom. For the atomic hydrogen anion then we get an upper bound of $$ E_{\text{H}^-} = \frac{121}{127} E_1 $$ As Griffiths points out, this is greater than (mind the sign!) the energy of the hydrogen atom, so our trial wavefunction is not able to establish whether the anionic species even exists.

Q2 G8.9

If we examine the (scaled) direct integral $D(R)$ $$ \int d\tau \frac{\left( +e \right)\left[ |\psi(\bm{x})|^2 \left(-e\right) \right]}{4\pi \epsilon_o |\bm{x}-\bm{x'}|} $$ where $\bm{x'}$ is the location of the opposite nucleus, we see that it is mathematically equivalent to the coulombic potential energy between a point charge of charge $+e$ and a charge distribution given by charge density $\rho(\bm{x})\equiv (-e) \left| \psi(\bm{x}) \right|^2$ (see, e.g., Fitzpatrick's notes). Knowing this, we are then free to use whatever tricks we remember from electrostatics to help us compute the integral.

We begin by conceptually breaking up the full charge distribution $-e|\psi(\bm{x})|^2$ into a Russian doll of spherical shells of thickness $dr$ and radii $0, dr, 2dr, \ldots, \infty$. If $dr$ is small enough, then the charge of the shell with radius $r$ is approximately $4\pi r^2 dr\left( (-e) |\psi(\bm{x})|^2 \right)$.

In intro physics we learned using Gauss' law that a hollow sphere of charge $q$ and radius $r$ produces an electric field outside the sphere equal to the field generated by a point charge of charge $q$ located at the center of the sphere, i.e. the field $\bm{E}(r';r,q)$ points radially outward and at a distance $r'$ from the center of the sphere has a magnitude equal to $$ |E(r';r,q)| = \begin{cases} \frac{q}{4 \pi \epsilon_o r'^2} & r' \gt r \\ 0 & r' \lt r \end{cases} $$ The electric potential $\Phi(r';r,q)$ should, by the same token, look like a point charge outside the shell and should be constant inside, i.e. $$ \Phi(r';r,q) = \frac{q}{4 \pi \epsilon_o} \begin{cases} \frac{1}{r'} & r' \gt r \\ \frac{1}{r} & r' \lt r \end{cases} $$ To determine the full electric potential generated by the electron's charge distribution, then we can, by the principle of superposition, integrate $$ d\phi(r;R) \equiv \Phi\left(R;r,4\pi r^2 dr \left[ \left(-e\right) | \psi(r) |^2 \right]\right) $$ from $r = 0 \to \infty$. If we multiply the result by the nuclear charge $+e$ we get the direct integral, i.e. $$ D(R) = +e \int d\phi(r;R) = +e \int_0^\infty \Phi\left(R;r,4\pi r^2 dr \left[ \left(-e\right) | \psi(r) |^2 \right]\right) $$ Defining $z \equiv r/a$ and $Z \equiv R/a$, the integral becomes $$ D(R) = - \frac{e^2}{\pi\epsilon_o a} \int_0^{\infty} dz \, z^2 e^{-2z} \cdot \begin{cases} \frac{1}{Z} & Z \gt z \\ \frac{1}{z} & Z \lt z \end{cases} $$ We can simplify a bit further defining $y \equiv -2z$ and $Y \equiv -2Z$. We then get $$ D(R) = - \frac{e^2}{4\pi\epsilon_o a} \int_0^{-\infty} dy \, y^2 e^{+y} \cdot \begin{cases} \frac{1}{Y} & Y \gt y \\ \frac{1}{y} & Y \lt y \end{cases} $$ or equivalently $$ D(R) = - \frac{e^2}{4\pi\epsilon_o a} \left( \frac{1}{Y} \int_0^Y dy \, y^2 e^{y} + \int_Y^{-\infty} dy \, y e^{y} \right) $$ This form is amenable to integration by parts, since $e_y dy = d(e^y)$. e.g. for the second integral we note \begin{align} & y e^y dy \\ = & y d(e^y) \\ = & d( y e^y ) - (dy) e_y \\ = & d( y e^y ) - d(e_y) \\ = & d\left( e^y ( y - 1 ) \right) \end{align} and then for the first integral \begin{align} & y^2 e^y dy \\ = & y^2 d(e^y) \\ = & d( y^2 e^y ) - d(y^2) e^y \\ = & d( y^2 e^y ) - 2 y e^y dy \\ = & d( y^2 e^y ) - 2 y d(e^y) \\ = & d( y^2 e^y ) - 2 \left( d( y e^y ) - (dy) e^y \right) \\ = & d( y^2 e^y ) - 2 \left( d( y e^y ) - d(e^y) \right) \\ = & d\left( e^y \left( y^2 - 2 \left( y - 1 \right) \right) \right) \end{align} so, at last $$ D(R) = - \frac{e^2}{4\pi\epsilon_o a} \left( \frac{1}{Y} \left. \left( e^y \left( y^2 - 2 \left( y - 1 \right) \right) \right) \right|^Y_0 + \left. \left( e^y \left( y - 1 \right) \right) \right|^{-\infty}_Y \right) $$ Evaluating, collecting terms, and substituting, we arrive at the result $$ D(R) = \frac{a}{R}-e^{-2R/a}\left( 1 + \frac{a}{R} \right) $$ The exchange integral was pretty unpleasant. I wrote it out on paper and uploaded a scan here.

Q3 G8.11

I wrote a script in python which finds the minimum using newton's method. My results are: So we see that many (≈ 8) vibrational levels fit in the well.

Here's the code I used to find the minimum and compute the final quantities 
      
import math
proton_electron_mass_ratio = 1836.152
initial_guess = 1.0
dx = 1.E-7
error = 1.E-7
def F(x):
    return (
        -1.0 + 2.0 / x * (
            (
                (
                    1.0 - 2.0 / 3.0 * x**2
                ) * math.exp(-1.0 * x) + (
                    1.0 + x
                ) * math.exp(-2.0 * x)
            ) / (
                1.0 + (
                    1.0 + x + 1.0 / 3.0 * x**2
                ) * math.exp(-1.0 * x)
            )
        )
    )

def dFdx(x):
    return (F(x+.5*dx) - F(x-.5*dx))/dx
def d2Fdx2(x):
    return (dFdx(x+.5*dx) - dFdx(x-.5*dx))/dx
def format_number(guess):
    return (
        r'%.'+str(int(-math.log10(error))) + 'e'
    ) % guess
guess = initial_guess
while True:
    print 'current guess:', format_number(guess)
    slope = dFdx(guess)
    dslope = d2Fdx2(guess)
    new_guess = guess - slope / dslope
    if abs(new_guess - guess) < error:
        break
    guess = new_guess
print 'final result:', format_number(new_guess), '+-', '%.0e' % dx, '(in bohr radii)'
print 'well depth: (in hartrees)', format_number(
    (F(new_guess)-F(1000.0))/2.0
)
print 'hbar omega (in hartrees):', format_number(
    math.sqrt(d2Fdx2(new_guess)/proton_electron_mass_ratio)
)
Here is a graphical illustration of newton's method applied to this problem:
Notice how the algorithm's relative rate of convergence increases as our guess approaches the derivative's zero crossing. On the last iteration we also start to run into the limitations of the software's floating point arithmatic. Note that the expression for the vibrational quantum $\hbar \omega$, in units of hartrees, takes the simplified to the form $$ \hbar \omega = \sqrt{\frac{m}{M}F''(r_{\text{eq}})} $$ Where $\frac{m}{M}$ is the electron-proton mass ratio (≈ 1836, the year of the Alamo), and $F''(r_{\text{eq}})$ is the second derivative of $F$ evaluated at the equilibrium separation, in units of $\frac{E_1}{a_o^2}$. Since $F''(r_{\text{eq}})$ is of order 1, this tells us that the vibrational levels of a molecule have separations roughly $\sqrt{\frac{m}{M}}$ times lower than the separation of the electronic levels. Can you work out the scale of rotational levels?

footnotes

11 hartree = 2 rydbergs = 2 $E_1 \approx 13.6$eV