## Tidal effect for objects in freefall near the Earth’s surface

Reference: Moore, Thomas A., A General Relativity Workbook, University Science Books (2013) – Chapter 18; Problem P18.1.

We’ve seen that in Newtonian physics, the tidal effect produces a relative acceleration between two objects in free fall above a sphere of mass ${M}$ given by

 $\displaystyle \frac{d^{2}n^{x}}{dt^{2}}$ $\displaystyle =$ $\displaystyle -\frac{GM}{r^{3}}n^{x}\ \ \ \ \ (1)$ $\displaystyle \frac{d^{2}n^{y}}{dt^{2}}$ $\displaystyle =$ $\displaystyle -\frac{GM}{r^{3}}n^{y}\ \ \ \ \ (2)$ $\displaystyle \frac{d^{2}n^{z}}{dt^{2}}$ $\displaystyle =$ $\displaystyle 2\frac{GM}{r^{3}}n^{z} \ \ \ \ \ (3)$

where ${\mathbf{n}}$ is the relative separation vector of the two objects, and the radial direction is taken to be along the ${z}$ axis.

To illustrate this, we can consider a situation where we have two masses separated initially by a vertical distance of 1 m and then released from rest near the Earth’s surface. How long will it take for the distance between these objects to increase by 1 nanometre due to the tidal effect?

We’ll start with the approximation that ${r}$ is a constant equal to the radius of the Earth, or ${r=6.378\times10^{6}\mbox{ m}}$. We can check the answer for consistency with this assumption after we’re done. We’ll also need the mass of the Earth: ${M=5.98\times10^{24}\mbox{ kg}}$ and of course ${G=6.67\times10^{-11}}$ in MKS units. We’re thus faced with the differential equation:

$\displaystyle \ddot{n^{z}}=3.075\times10^{-6}n^{z} \ \ \ \ \ (4)$

This has the general solution

$\displaystyle n^{z}\left(t\right)=Ae^{0.00175t}+Be^{-0.00175t} \ \ \ \ \ (5)$

where the numerical factor in the exponent is ${\sqrt{3.075\times10^{-6}}}$. From initial conditions

 $\displaystyle n^{z}\left(0\right)$ $\displaystyle =$ $\displaystyle 1=A+B\ \ \ \ \ (6)$ $\displaystyle \dot{n^{z}}\left(0\right)$ $\displaystyle =$ $\displaystyle 0=0.00175\left(A-B\right) \ \ \ \ \ (7)$

Therefore

$\displaystyle A=B=\frac{1}{2} \ \ \ \ \ (8)$

We want the time ${t}$ such that ${n^{z}\left(t\right)=1+10^{-9}\mbox{ m}}$ or in other words:

 $\displaystyle \frac{1}{2}\left(e^{0.00175t}+e^{-0.00175t}\right)-1$ $\displaystyle =$ $\displaystyle 10^{-9}\ \ \ \ \ (9)$ $\displaystyle \cosh0.00175t-1$ $\displaystyle =$ $\displaystyle 10^{-9} \ \ \ \ \ (10)$

Since the difference is so small, we can expand the argument of the cosh in a Taylor series and keep the leading term:

 $\displaystyle \cosh0.00175t-1$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(0.00175t\right)^{2}+\mathcal{O}\left(t^{4}\right)\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 10^{-9}\ \ \ \ \ (12)$ $\displaystyle t$ $\displaystyle =$ $\displaystyle 0.0255\mbox{ s} \ \ \ \ \ (13)$

In that time, the objects would fall a distance of

$\displaystyle d=\frac{1}{2}gt^{2}=0.003\mbox{ m}$

or about 3 mm, so the assumption of ${r}$ being constant at the Earth’s radius seems reasonable.

For the separation to reach 1 mm, the same calculation yields a time of around 25.5 sec, which would amount to a distance of (assuming the acceleration due to gravity ${g}$ is constant and ignoring air resistance) of around 3.2 km, so here the assumption of a constant ${r}$ is a bit shakier, considering we’re dealing with such a small tidal effect.

## Newtonian tidal effect

Reference: Moore, Thomas A., A General Relativity Workbook, University Science Books (2013) – Chapter 18; Box 18.1.

As a prelude to analyzing geodesic deviation in general relativity, we’ll look at the tidal effect in Newtonian physics. First, though, what is geodesic deviation? A geodesic curve is a curve followed by an object in the absence of any forces (we’re not considering gravity as a force here since its effects define the curvature of spacetime through which the object moves). Unless spacetime is flat, two objects initially following parallel geodesics will accelerate towards or away from each other.

This effect occurs in Newtonian physics as well, although it’s not analyzed by considering spacetime to be curved. One illustration of this is to consider a box in freefall above the earth. We’ll put an observer at the centre of the box, and place balls above, below and to either side of the observer. Because the gravitational force (we can think of gravity as a force if we’re using Newtonian physics) due to the Earth is radial and falls off according to the inverse square formula, the ball at the top of the box experiences a force (per unit mass) less than that of the observer, who in turn feels a smaller force than the ball at the bottom of the box. Thus the ball at the top accelerates more slowly than the observer, who in turn accelerates more slowly than the ball at the bottom. To the observer, these two balls appear to be accelerating away from him in opposite directions.

Due to the radial direction of the gravitational force, the two balls on either side of the observer will feel a small component of the force pulling them towards the centre of the box (although most of the force is directed downwards, of course). As a result, the observer will see these two balls accelerating towards him. The five objects inside the box start out on parallel geodesics, but each of them follows its own geodesic, which deviates from the initial parallel path.

Qualitatively, we can see that it is this effect which causes the tides in the Earth’s oceans. The Earth is in freefall about the Earth-Moon centre of mass, so points on the side of the Earth opposite the Moon have a smaller acceleration than the Earth’s centre, which in turn has a smaller acceleration than points on the side closest to the Moon. Points on the surfaces of the Earth lying along the normal to the line connecting the Earth and Moon experience accelerations towards the centre of the Earth. This gives rise to high tides on the near and far points on the Earth, and low tides in between.

To make the argument quantitative, we can introduce the gravitational potential

$\displaystyle \Phi=-\frac{GM}{r} \ \ \ \ \ (1)$

The acceleration, or force per unit mass, is the negative gradient of ${\Phi}$ which we can write using index notation and the flat space metric ${\eta^{ij}}$ (which is just ${\eta^{ij}=\delta^{ij}}$ here):

$\displaystyle \frac{d^{2}x^{i}}{dt^{2}}=-\eta^{ij}\left.\frac{\partial\Phi}{\partial x^{j}}\right|_{\vec{x}} \ \ \ \ \ (2)$

All quantities are evaluated at the object’s position ${\vec{x}\left(t\right)}$. Note that everywhere in this post, the summation over indices is over only the spatial coordinates.

Now suppose we have a second object which starts off an infinitesimal separation ${\vec{n}\left(t\right)}$ away from the first object, so its position is ${x^{i}\left(t\right)+n^{i}\left(t\right)}$. The first object could be our observer in the box, and the second object could be one of the balls. To see how geodesic deviation (or the tidal effect) works, we’d like to find the relative acceleration, which is ${d^{2}n^{i}/dt^{2}}$.

The acceleration of the second object is

$\displaystyle \frac{d^{2}\left(x^{i}+n^{i}\right)}{dt^{2}}=-\eta^{ij}\left.\frac{\partial\Phi}{\partial x^{j}}\right|_{\vec{x}+\vec{n}} \ \ \ \ \ (3)$

where this time the gradient is evaluated at ${\vec{x}+\vec{n}}$. If ${\vec{n}}$ is small, we can expand this gradient in a Taylor series about ${\vec{x}}$:

$\displaystyle \left.\frac{\partial\Phi}{\partial x^{j}}\right|_{\vec{x}+\vec{n}}=\left.\frac{\partial\Phi}{\partial x^{j}}\right|_{\vec{x}}+n^{i}\left.\frac{\partial^{2}\Phi}{\partial x^{i}\partial x^{j}}\right|_{\vec{x}}+... \ \ \ \ \ (4)$

Putting this into 3 we get, ignoring higher order terms:

 $\displaystyle \frac{d^{2}\left(x^{i}+n^{i}\right)}{dt^{2}}$ $\displaystyle =$ $\displaystyle -\eta^{ij}\left(\left.\frac{\partial\Phi}{\partial x^{j}}\right|_{\vec{x}}+n^{k}\left.\frac{\partial^{2}\Phi}{\partial x^{k}\partial x^{j}}\right|_{\vec{x}}\right)\ \ \ \ \ (5)$ $\displaystyle \frac{d^{2}x^{i}}{dt^{2}}+\frac{d^{2}n^{i}}{dt^{2}}$ $\displaystyle =$ $\displaystyle -\eta^{ij}\left.\frac{\partial\Phi}{\partial x^{j}}\right|_{\vec{x}}-\eta^{ij}n^{k}\left.\frac{\partial^{2}\Phi}{\partial x^{k}\partial x^{j}}\right|_{\vec{x}}\ \ \ \ \ (6)$ $\displaystyle \frac{d^{2}n^{i}}{dt^{2}}$ $\displaystyle =$ $\displaystyle -\eta^{ij}n^{k}\left.\frac{\partial^{2}\Phi}{\partial x^{k}\partial x^{j}}\right|_{\vec{x}} \ \ \ \ \ (7)$

where in the last line, we used 2. This is the Newtonian tidal deviation equation, expressed using index notation.

We can use this equation to work out the relative acceleration in the three rectangular directions for the falling box above, taking the ${z}$ direction to be radially outward. Using index notation, we have

$\displaystyle r=\sqrt{\eta_{mn}x^{m}x^{n}} \ \ \ \ \ (8)$

The gradient terms can be written as

 $\displaystyle \frac{\partial\Phi}{\partial x^{j}}$ $\displaystyle =$ $\displaystyle \frac{d\Phi}{dr}\frac{\partial r}{\partial x^{j}}\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{GM}{r^{2}}\frac{\eta_{jn}x^{n}+\eta_{mj}x^{m}}{2\sqrt{\eta_{mn}x^{m}x^{n}}}\ \ \ \ \ (10)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{GM}{r^{3}}\eta_{jn}x^{n} \ \ \ \ \ (11)$

where we got the last line using the symmetry of the metric ${\eta_{ij}=\eta_{ji}}$.

The second derivatives are then:

 $\displaystyle \partial_{k}\partial_{j}\Phi$ $\displaystyle =$ $\displaystyle -\frac{3GM}{r^{4}}\frac{\partial r}{\partial x^{k}}\eta_{jn}x^{n}+\frac{GM}{r^{3}}\eta_{jk}\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{3GM}{r^{4}}\frac{\eta_{km}x^{m}}{r}\eta_{jn}x^{n}+\frac{GM}{r^{3}}\eta_{jk}\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{3GM}{r^{5}}\eta_{km}\eta_{jn}x^{n}x^{m}+\frac{GM}{r^{3}}\eta_{jk} \ \ \ \ \ (14)$

Putting this into 7, and using ${\eta^{ij}\eta_{jn}=\delta_{n}^{i}}$:

 $\displaystyle \frac{d^{2}n^{i}}{dt^{2}}$ $\displaystyle =$ $\displaystyle -\eta^{ij}n^{k}\left[-\frac{3GM}{r^{5}}\eta_{km}\eta_{jn}x^{n}x^{m}+\frac{GM}{r^{3}}\eta_{jk}\right]\ \ \ \ \ (15)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3GM}{r^{5}}\eta_{km}n^{k}x^{m}x^{i}-\frac{GM}{r^{3}}n^{i} \ \ \ \ \ (16)$

If we take the reference object to be at the position ${x=y=0}$, ${z=r}$, then the first term is non-zero only when all the coordinates are ${z}$ coordinates so

 $\displaystyle \frac{d^{2}n^{x}}{dt^{2}}$ $\displaystyle =$ $\displaystyle -\frac{GM}{r^{3}}n^{x}\ \ \ \ \ (17)$ $\displaystyle \frac{d^{2}n^{y}}{dt^{2}}$ $\displaystyle =$ $\displaystyle -\frac{GM}{r^{3}}n^{y}\ \ \ \ \ (18)$ $\displaystyle \frac{d^{2}n^{z}}{dt^{2}}$ $\displaystyle =$ $\displaystyle 2\frac{GM}{r^{3}}n^{z} \ \ \ \ \ (19)$

Thus the acceleration in the transverse directions (${x}$ and ${y}$) is opposite to the relative displacement so that the objects accelerate towards each other, and in the radial (${z}$) direction the acceleration is away from each other.

## Quantum dots

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 7.20.

An unusual problem involving the variational principle is as follows. Suppose a particle is constrained to move in two dimensions in a cross-shaped region. We can specify the cross shape by recognizing that it is symmetric in each octant of the ${xy}$ plane. In the first octant (bounded by the ${x}$ axis and the line ${y=x}$) the allowed region is bounded by the ${x}$ axis, the line ${y=x}$ and the horizontal line ${y=+a}$. This region is then reflected about the line ${y=x}$ to get the allowed region in the second octant, and then the total region is the first quadrant is replicated in the other three quadrants to get the cross.

What kinds of energy levels are allowed in such a system? We can consider first the limiting case where the particle is far out along the arm of the cross extending towards ${+x}$. Out here, we have essentially an infinite square well of width ${2a}$ in the ${y}$ direction, and a free particle in the ${x}$ direction. That is, we can write the Schrödinger equation as

$\displaystyle -\frac{\hbar^{2}}{2m}\left(\frac{d^{2}}{dx^{2}}+\frac{d^{2}}{dy^{2}}\right)\psi=E\psi \ \ \ \ \ (1)$

and use separation of variables so that ${\psi=X\left(x\right)Y\left(y\right)}$. The ${y}$ equation is just that of the infinite square well, so its energy levels are given by

$\displaystyle E_{y}=\frac{\left(n\pi\hbar\right)^{2}}{2m\left(2a\right)^{2}} \ \ \ \ \ (2)$

The energy due to the free particle contribution in the ${x}$ direction must be positive, so any particle that can escape to infinity must have an energy greater than ${E_{y}}$; in other words, a particle cannot escape to infinity (that is, it’s in a bound state) if its energy is lower than the ground state (${n=1}$) of the square well:

$\displaystyle E<\frac{\left(\pi\hbar\right)^{2}}{8ma^{2}} \ \ \ \ \ (3)$

Clearly, trying to calculate the exact energy levels for a potential with such an odd shape would be very difficult, but we can use the variational principle to get an upper bound on this energy. The trial function is

$\displaystyle \psi=\begin{cases} A\left(1-\frac{\left|xy\right|}{a^{2}}\right)e^{-\alpha} & \left|x\right|\le a\mbox{ and }\left|y\right|\le a\\ A\left(1-\frac{\left|x\right|}{a}\right)e^{-\alpha\left|y\right|/a} & \left|x\right|\le a\mbox{ and }\left|y\right|>a\\ A\left(1-\frac{\left|y\right|}{a}\right)e^{-\alpha\left|x\right|/a} & \left|x\right|>a\mbox{ and }\left|y\right|\le a\\ 0 & \mbox{otherwise} \end{cases} \ \ \ \ \ (4)$

The parameter ${\alpha}$ is the variational parameter. The first line is in the square of side ${2a}$ centred on the origin, the second line is in the top and bottom arms of the cross and the third line is in the left and right arms. The wave function is of course zero outside the cross, since the potential is infinite there and the particle cannot exist there.

First, we need to find ${A}$ from normalization. We can use the symmetry of the problem and of ${\psi}$ to simplify the integrals, since the problem is symmetric in all 8 octants. Thus we can integrate over the first octant and multiply the result by 8. Therefore

 $\displaystyle 8\left[\int_{0}^{a}\int_{y}^{a}+\int_{0}^{a}\int_{a}^{\infty}\right]\psi^{2}dx\; dy$ $\displaystyle =$ $\displaystyle 1\ \ \ \ \ (5)$ $\displaystyle \int_{0}^{a}\int_{y}^{a}\left(1-\frac{xy}{a^{2}}\right)^{2}e^{-2\alpha}dx\; dy+\int_{0}^{a}\int_{a}^{\infty}\left(1-\frac{y}{a}\right)^{2}e^{-2\alpha x/a}dx\; dy$ $\displaystyle =$ $\displaystyle \frac{1}{8A^{2}}\ \ \ \ \ (6)$ $\displaystyle A$ $\displaystyle =$ $\displaystyle \frac{3e^{\alpha}}{2a}\sqrt{\frac{2\alpha}{11\alpha+6}} \ \ \ \ \ (7)$

Now to work out ${\left\langle H\right\rangle }$ we need the derivatives of ${\psi}$. We need to be careful at the boundaries between the regions, since although ${\psi}$ is continuous there, its first derivative isn’t, so the second derivative will contain delta functions. We get, again considering only the first octant:

$\displaystyle \frac{1}{A}\frac{\partial\psi}{\partial x}=\begin{cases} -\frac{y}{\alpha^{2}}e^{-\alpha} & xa \end{cases} \ \ \ \ \ (8)$

$\displaystyle \frac{1}{A}\frac{\partial^{2}\psi}{\partial x^{2}}=\begin{cases} 0 & xa \end{cases} \ \ \ \ \ (9)$

The first derivative is discontinuous along the line ${x=a}$, so there we get

$\displaystyle \frac{1}{A}\frac{\partial\psi}{\partial x}=\begin{cases} -\frac{y}{\alpha^{2}}e^{-\alpha} & x=a_{-}\\ -\frac{\alpha}{a}\left(1-\frac{y}{a}\right)e^{-\alpha} & x=a_{+} \end{cases} \ \ \ \ \ (10)$

so we get a step change at ${x=a}$ in the amount of

$\displaystyle \left.\frac{\partial\psi}{\partial x}\right|_{+}-\left.\frac{\partial\psi}{\partial x}\right|_{-}=A\left[-\frac{\alpha}{a}\left(1-\frac{y}{a}\right)+\frac{y}{\alpha^{2}}\right]e^{-\alpha} \ \ \ \ \ (11)$

The derivative of a step function is a delta function, so at ${x=a}$ we have

$\displaystyle \left.\frac{\partial^{2}\psi}{\partial x^{2}}\right|_{x=a}=A\left[-\frac{\alpha}{a}\left(1-\frac{y}{a}\right)+\frac{y}{\alpha^{2}}\right]e^{-\alpha}\delta\left(x-a\right)$

and the complete derivative is then

$\displaystyle \frac{1}{A}\frac{\partial^{2}\psi}{\partial x^{2}}=\begin{cases} \left[-\frac{\alpha}{a}\left(1-\frac{y}{a}\right)+\frac{y}{\alpha^{2}}\right]e^{-\alpha}\delta\left(x-a\right) & x\le a\\ \frac{\alpha^{2}}{a^{2}}\left(1-\frac{y}{a}\right)e^{-\alpha x/a} & x>a \end{cases} \ \ \ \ \ (12)$

In the ${y}$ direction, we get, for ${x

$\displaystyle \frac{1}{A}\frac{\partial\psi}{\partial y}=\begin{cases} -\frac{x}{a^{2}}e^{-\alpha} & y>0\\ \frac{x}{a^{2}}e^{-\alpha} & y<0 \end{cases} \ \ \ \ \ (13)$

so for ${y\ne0}$, we have

$\displaystyle \frac{\partial^{2}\psi}{\partial y^{2}}=0 \ \ \ \ \ (14)$

Because of the discontinuity at ${y=0}$ we get another delta function, so the total second derivative is, for ${x

$\displaystyle \frac{\partial^{2}\psi}{\partial y^{2}}=-\frac{2x}{a^{2}}e^{-\alpha}\delta\left(y\right) \ \ \ \ \ (15)$

For ${x>a}$

$\displaystyle \frac{1}{A}\frac{\partial\psi}{\partial y}=\begin{cases} -\frac{1}{a^{2}}e^{-\alpha x/a} & y>0\\ \frac{1}{a^{2}}e^{-\alpha x/a} & y<0 \end{cases} \ \ \ \ \ (16)$

Again, the second derivative is zero on both sides of the ${x}$ axis, but there is a delta function on the axis:

$\displaystyle \frac{\partial^{2}\psi}{\partial y^{2}}=-\frac{2A}{a^{2}}e^{-\alpha x/a}\delta\left(y\right) \ \ \ \ \ (17)$

We can now put all this together to calculate ${\left\langle H\right\rangle }$. First, the term excluding the delta functions. Since both second derivatives are zero for ${x, there is only the contribution from ${x>a}$:

$\displaystyle H_{1}=-\frac{A^{2}\hbar^{2}}{2m}\int_{0}^{a}\int_{a}^{\infty}\frac{\alpha^{2}}{a^{2}}\left(1-\frac{y}{a}\right)^{2}e^{-2\alpha x/a}dx\; dy \ \ \ \ \ (18)$

Next, the term resulting from the delta function in the ${x}$ derivative. We set ${x=a}$ and then integrate over ${y}$:

$\displaystyle H_{2}=-\frac{A^{2}\hbar^{2}}{2m}e^{-2\alpha}\int_{0}^{a}\left(1-\frac{y}{a}\right)\left[-\frac{\alpha}{a}\left(1-\frac{y}{a}\right)+\frac{y}{\alpha^{2}}\right]dy \ \ \ \ \ (19)$

Finally, the terms resulting from the delta functions in the ${y}$ derivative. We set ${y=0}$ and integrate over ${x}$:

$\displaystyle H_{3}=\frac{1}{2}\left[\frac{A^{2}\hbar^{2}}{ma^{2}}e^{-2\alpha}\int_{0}^{a}x\; dx-\frac{A^{2}\hbar^{2}}{ma}\int_{a}^{\infty}e^{-2\alpha x/a}dx\right] \ \ \ \ \ (20)$

The factor of ${\frac{1}{2}}$ in ${H_{3}}$ comes from the fact that the ${x}$ axis is shared between two octants, so only half the integral over the delta function contributes to the first octant.

We can evaluate these integrals using Maple (or by hand; none of them is difficult, but there’s a lot of calculation) to find, after simplifying:

 $\displaystyle \left\langle H\right\rangle$ $\displaystyle =$ $\displaystyle 8\left(H_{1}+H_{2}+H_{3}\right)\ \ \ \ \ (21)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3\hbar^{2}}{ma^{2}}\left(\frac{\alpha^{2}+2\alpha+3}{6+11\alpha}\right) \ \ \ \ \ (22)$

We can find the value of ${\alpha}$ that minimizes this by taking the derivative as usual, and we get

$\displaystyle \alpha_{min}=\frac{1}{11}\left(-6\pm\sqrt{267}\right) \ \ \ \ \ (23)$

We require ${\alpha>0}$ in order for the integrals above to converge, so we must take the positive root. Plugging this back into ${\left\langle H\right\rangle }$ we find

$\displaystyle \left\langle H\right\rangle =\frac{1.058\hbar^{2}}{ma^{2}}<\frac{\pi^{2}\hbar^{2}}{8ma^{2}}=\frac{1.234\hbar^{2}}{ma^{2}} \ \ \ \ \ (24)$

The upper bound given by the variational principle is thus below the minimum energy at which a particle can escape, so there is at least one bound state in this system.

## Fusion with a muon-deuteron system

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 7.19.

When we examined the hydrogen molecule ion ${\mbox{H}_{2}^{+}}$ using the variational principle, we found that the lowest energy occurred when the protons were separated by ${R=2.4a}$, where ${a}$ is Bohr radius. It is this separation which causes the problem in attempting to achieve controlled nuclear fusion, since we need to get the two nuclei close enough for the short range nuclear force to overcome the electrical repulsion. A possible alternative is to use a deuteron molecule ion with the electron replaced by a muon, whose mass is 207 times the mass of the electron.

The only difference from the original calculation is in the masses of the particles. In this case, the Bohr radius is

$\displaystyle a_{\mu}=\frac{4\pi\epsilon_{0}\hbar^{2}}{\mu_{d\mu}e^{2}}$

where ${\mu_{d\mu}}$ is the reduced mass of the muon and deuteron. Taking the deuteron mass to be twice that of the proton, which is in turn 1833 times the mass of the electron, so

$\displaystyle \mu_{d\mu}=\frac{\left(207m_{e}\right)\left(2\times1833m_{e}\right)}{207m_{e}+2\times1833m_{e}}=196m_{e}$

Therefore, the Bohr radius is reduced by a factor of 196, meaning that the two deuterons are separated by ${R=2.4a_{\mu}=6.73\times10^{-13}\mbox{ m}}$. Having the two nuclei much closer together should make fusion easier to achieve.

## Variational principle and the hydrogen ion: two parameters

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 7.18.

An interesting variant on the variational principle used to calculate an upper bound on the ground state of the hydrogen ion is one proposed in 1944 by Chandrasekhar. It uses two adjustable parameters instead of the one we’ve used so far. The exact hamiltonian for the ${\mbox{H}^{-}}$ ion is

$\displaystyle H=-\frac{\hbar^{2}}{2m}\left(\nabla_{1}^{2}+\nabla_{2}^{2}\right)-\frac{e^{2}}{4\pi\epsilon_{0}}\left[\frac{1}{r_{1}}+\frac{1}{r_{2}}-\frac{1}{\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|}\right] \ \ \ \ \ (1)$

Here we are using two independent spatial coordinates ${\mathbf{r}_{1}}$ and ${\mathbf{r}_{2}}$, one for each electron. The terms in the square brackets are the interaction terms between the electrons and the nucleus and between the two electrons.

We use as a test function a combination of hydrogen ground states in the form

$\displaystyle \psi=A\left[\psi_{1}\left(r_{1}\right)\psi_{2}\left(r_{2}\right)+\psi_{1}\left(r_{2}\right)\psi_{2}\left(r_{1}\right)\right] \ \ \ \ \ (2)$

where

$\displaystyle \psi_{i}\left(r\right)=\sqrt{\frac{Z_{i}^{3}}{\pi a^{3}}}e^{-Z_{i}r/a} \ \ \ \ \ (3)$

is a normalized ground state hydrogen wave function and ${Z_{1},\; Z_{2}}$ are the parameters. The idea is that each of the two electrons experiences a different amount of shielding due to differing distances from the nucleus. Equation 2 is written as a symmetric function of ${r_{1}}$ and ${r_{2}}$ implying that the spin state is antisymmetric in order to give the two electrons an overall antisymmetric function, as required by the Pauli exclusion principle.

To find the optimal values of the ${Z_{i}}$ we can follow the usual procedure in the variational principle, although the calculations are quite tedious. One approach is simply to plug the equations as they stand into software such as Maple and let it grind through the calculations. However, to get the answer in Griffiths’s book requires a bit more work.

First, we observe that ${\psi_{i}\left(r_{j}\right)}$ is an eigenfunction of

$\displaystyle H_{ij}\equiv-\frac{\hbar^{2}}{2m}\nabla_{j}^{2}-\frac{e^{2}}{4\pi\epsilon_{0}}\frac{Z_{i}}{r_{j}} \ \ \ \ \ (4)$

with eigenvalue (energy) ${Z_{i}^{2}E_{1}}$ where ${E_{1}=-13.6\mbox{ eV}}$ is the ground state energy of a hydrogen atom. Thus we can write the original hamiltonian 1 as

 $\displaystyle H$ $\displaystyle =$ $\displaystyle H_{11}+H_{22}+\frac{e^{2}}{4\pi\epsilon_{0}}\left[\frac{Z_{1}-1}{r_{1}}+\frac{Z_{2}-1}{r_{2}}\right]+\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|}\ \ \ \ \ (5)$ $\displaystyle$ $\displaystyle =$ $\displaystyle H_{12}+H_{21}+\frac{e^{2}}{4\pi\epsilon_{0}}\left[\frac{Z_{2}-1}{r_{1}}+\frac{Z_{1}-1}{r_{2}}\right]+\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|} \ \ \ \ \ (6)$

In calculating ${H\psi}$, we can use the fact that

 $\displaystyle \left(H_{11}+H_{22}\right)\psi_{1}\left(r_{1}\right)\psi_{2}\left(r_{2}\right)$ $\displaystyle =$ $\displaystyle E_{1}\left(Z_{1}^{2}+Z_{2}^{2}\right)\psi_{1}\left(r_{1}\right)\psi_{2}\left(r_{2}\right)\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(H_{12}+H_{21}\right)\psi_{1}\left(r_{2}\right)\psi_{2}\left(r_{1}\right) \ \ \ \ \ (8)$

To calculate ${\left\langle H\right\rangle }$ we therefore need to calculate the means of the remaining terms in 5 and 6.

First, however, we need to calculate the normalization constant ${A}$ in 2. We get

$\displaystyle A^{2}\int\psi^{2}d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}=1 \ \ \ \ \ (9)$

Doing the integrals with Maple, we get

$\displaystyle 2A^{2}\,{\frac{{Z_{{1}}}^{6}+15\,{Z_{{1}}}^{4}{Z_{{2}}}^{2}+84\,{Z_{{2}}}^{3}{Z_{{1}}}^{3}+6\,{Z_{{2}}}^{5}Z_{{1}}+15\,{Z_{{2}}}^{4}{Z_{{1}}}^{2}+{Z_{{2}}}^{6}+6\,{Z_{{1}}}^{5}Z_{{2}}}{{Z_{{1}}}^{6}+6\,{Z_{{1}}}^{5}Z_{{2}}+15\,{Z_{{1}}}^{4}{Z_{{2}}}^{2}+20\,{Z_{{2}}}^{3}{Z_{{1}}}^{3}+15\,{Z_{{2}}}^{4}{Z_{{1}}}^{2}+6\,{Z_{{2}}}^{5}Z_{{1}}+{Z_{{2}}}^{6}}}=1 \ \ \ \ \ (10)$

This rather frightening expression can be simplified by factoring the numerator and denominator. First, the numerator:

 $\displaystyle {Z_{{1}}}^{6}+15\,{Z_{{1}}}^{4}{Z_{{2}}}^{2}+84\,{Z_{{2}}}^{3}{Z_{{1}}}^{3}+6\,{Z_{{2}}}^{5}Z_{{1}}+15\,{Z_{{2}}}^{4}{Z_{{1}}}^{2}+{Z_{{2}}}^{6}+6\,{Z_{{1}}}^{5}Z_{{2}}$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (11)$ $\displaystyle \left({Z_{{2}}}^{2}+6\, Z_{{1}}Z_{{2}}+{Z_{{1}}}^{2}\right)\left({Z_{{2}}}^{4}+14\,{Z_{{2}}}^{2}{Z_{{1}}}^{2}+{Z_{{1}}}^{4}\right)$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (12)$ $\displaystyle \left[\left(Z_{1}+Z_{2}\right)^{2}+4Z_{1}Z_{2}\right]\left[\left(Z_{1}+Z_{2}\right)^{4}-4Z_{1}^{3}Z_{2}+8Z_{1}^{2}Z_{2}^{2}-4Z_{1}Z_{2}^{3}\right]$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (13)$ $\displaystyle \left[\left(Z_{1}+Z_{2}\right)^{2}+4Z_{1}Z_{2}\right]\left[\left(Z_{1}+Z_{2}\right)^{4}+16Z_{1}^{2}Z_{2}^{2}-4Z_{1}^{3}Z_{2}-8Z_{1}^{2}Z_{2}^{2}-4Z_{1}Z_{2}^{3}\right]$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (14)$ $\displaystyle \left[\left(Z_{1}+Z_{2}\right)^{2}+4Z_{1}Z_{2}\right]\left[\left(Z_{1}+Z_{2}\right)^{4}+16Z_{1}^{2}Z_{2}^{2}-4Z_{1}Z_{2}\left(Z_{1}^{2}+2Z_{1}Z_{2}+Z_{2}^{2}\right)\right]$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (15)$ $\displaystyle \left[\left(Z_{1}+Z_{2}\right)^{2}+4Z_{1}Z_{2}\right]\left[\left(Z_{1}+Z_{2}\right)^{4}+16Z_{1}^{2}Z_{2}^{2}-4Z_{1}Z_{2}\left(Z_{1}+Z_{2}\right)^{2}\right]$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (16)$ $\displaystyle \left(x^{2}+y^{2}\right)\left[x^{4}+y^{4}-x^{2}y^{2}\right]$ $\displaystyle =$ $\displaystyle x^{6}+y^{6} \ \ \ \ \ (17)$

where we’ve defined

 $\displaystyle x$ $\displaystyle \equiv$ $\displaystyle Z_{1}+Z_{2}\ \ \ \ \ (18)$ $\displaystyle y$ $\displaystyle \equiv$ $\displaystyle 2\sqrt{Z_{1}Z_{2}} \ \ \ \ \ (19)$

The denominator comes out to

$\displaystyle {Z_{{1}}}^{6}+6\,{Z_{{1}}}^{5}Z_{{2}}+15\,{Z_{{1}}}^{4}{Z_{{2}}}^{2}+20\,{Z_{{2}}}^{3}{Z_{{1}}}^{3}+15\,{Z_{{2}}}^{4}{Z_{{1}}}^{2}+6\,{Z_{{2}}}^{5}Z_{{1}}+{Z_{{2}}}^{6}=x^{6} \ \ \ \ \ (20)$

so

$\displaystyle A^{2}=\frac{x^{6}}{2\left(x^{6}+y^{6}\right)} \ \ \ \ \ (21)$

Now we need to grind through the calculations for the terms in 5 and 6. First, we’ll do the electron-electron interaction term:

 $\displaystyle \left\langle V_{ee}\right\rangle$ $\displaystyle =$ $\displaystyle \frac{e^{2}A^{2}}{4\pi\epsilon_{0}}\int\left[\psi_{1}\left(r_{1}\right)\psi_{2}\left(r_{2}\right)+\psi_{1}\left(r_{2}\right)\psi_{2}\left(r_{1}\right)\right]^{2}\frac{d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}}{\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|}\ \ \ \ \ (22)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{e^{2}A^{2}}{4\pi\epsilon_{0}}\int\left[\psi_{1}\left(r_{1}\right)\psi_{2}\left(r_{2}\right)+\psi_{1}\left(r_{2}\right)\psi_{2}\left(r_{1}\right)\right]^{2}\frac{d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}}{\sqrt{r_{1}^{2}+r_{2}^{2}-2r_{1}r_{2}\cos\theta_{2}}}\ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{e^{2}A^{2}}{2\times4\pi\epsilon_{0}a}{\frac{{y}^{2}\left(5\,{Z_{{2}}}^{3}Z_{{1}}+5\,{Z_{{1}}}^{3}Z_{{2}}+28\,{Z_{{2}}}^{2}{Z_{{1}}}^{2}+{Z_{{1}}}^{4}+{Z_{{2}}}^{4}\right)}{\left({Z_{{1}}}^{5}+5\,{Z_{{1}}}^{4}Z_{{2}}+10\,{Z_{{1}}}^{3}{Z_{{2}}}^{2}+10\,{Z_{{1}}}^{2}{Z_{{2}}}^{3}+5\,{Z_{{2}}}^{4}Z_{{1}}+{Z_{{2}}}^{5}\right)}}\ \ \ \ \ (24)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -E_{1}A^{2}\frac{{y}^{2}\left(5\,{Z_{{2}}}^{3}Z_{{1}}+5\,{Z_{{1}}}^{3}Z_{{2}}+28\,{Z_{{2}}}^{2}{Z_{{1}}}^{2}+{Z_{{1}}}^{4}+{Z_{{2}}}^{4}\right)}{x^{5}} \ \ \ \ \ (25)$

To reduce the numerator:

 $\displaystyle 5\,{Z_{{2}}}^{3}Z_{{1}}+5\,{Z_{{1}}}^{3}Z_{{2}}+28\,{Z_{{2}}}^{2}{Z_{{1}}}^{2}+{Z_{{1}}}^{4}+{Z_{{2}}}^{4}$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (26)$ $\displaystyle x^{4}+{Z_{{2}}}^{3}Z_{{1}}+{Z_{{1}}}^{3}Z_{{2}}+22\,{Z_{{2}}}^{2}{Z_{{1}}}^{2}$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (27)$ $\displaystyle x^{4}+Z_{1}Z_{2}\left(Z_{1}^{2}+Z_{2}^{2}+2Z_{1}Z_{2}+20Z_{1}Z_{2}\right)$ $\displaystyle =$ $\displaystyle x^{4}+\frac{y^{2}}{4}\left(x^{2}+5y^{2}\right) \ \ \ \ \ (28)$

Therefore

 $\displaystyle \left\langle V_{ee}\right\rangle$ $\displaystyle =$ $\displaystyle -E_{1}A^{2}y^{2}\frac{4x^{4}+y^{2}\left(x^{2}+5y^{2}\right)}{4x^{5}}\ \ \ \ \ (29)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{1}{8}{\frac{E_{{1}}x{y}^{2}\left(4\,{x}^{4}+5\,{y}^{4}+{x}^{2}{y}^{2}\right)}{{x}^{6}+{y}^{6}}} \ \ \ \ \ (30)$

Now for the electron-nucleus terms in 5 and 6.

 $\displaystyle I_{1}$ $\displaystyle =$ $\displaystyle \frac{e^{2}}{4\pi\epsilon_{0}}\int\psi\left[\frac{Z_{1}-1}{r_{1}}+\frac{Z_{2}-1}{r_{2}}\right]\psi_{1}\left(r_{1}\right)\psi_{2}\left(r_{2}\right)d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}\ \ \ \ \ (31)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -2E_{1}A^{2}\frac{N}{x^{5}} \ \ \ \ \ (32)$

The numerator is

 $\displaystyle N$ $\displaystyle =$ $\displaystyle {Z_{{2}}}^{7}+5\,{Z_{{1}}}^{6}Z_{{2}}-84\,{Z_{{2}}}^{3}{Z_{{1}}}^{3}+11\,{Z_{{1}}}^{2}{Z_{{2}}}^{5}-6\,{Z_{{1}}}^{5}Z_{{2}}+47\,{Z_{{1}}}^{4}{Z_{{2}}}^{3}-{Z_{{1}}}^{6}+11\,{Z_{{1}}}^{5}{Z_{{2}}}^{2}\nonumber$ $\displaystyle$ $\displaystyle$ $\displaystyle -15\,{Z_{{1}}}^{4}{Z_{{2}}}^{2}-{Z_{{2}}}^{6}-6\,{Z_{{2}}}^{5}Z_{{1}}+47\,{Z_{{1}}}^{3}{Z_{{2}}}^{4}+5\, Z_{{1}}{Z_{{2}}}^{6}+{Z_{{1}}}^{7}-15\,{Z_{{2}}}^{4}{Z_{{1}}}^{2} \ \ \ \ \ (33)$

We can break this down by picking out the terms where the exponents sum to 6 (denoted by ${N_{6}}$) and then the terms where they sum to 7 (${N_{7}}$). The terms summing to 6 are:

 $\displaystyle N_{6}$ $\displaystyle =$ $\displaystyle -{Z_{{2}}}^{6}-6\,{Z_{{1}}}^{5}Z_{{2}}-{Z_{{1}}}^{6}-15\,{Z_{{1}}}^{4}{Z_{{2}}}^{2}-6\,{Z_{{2}}}^{5}Z_{{1}}-84\,{Z_{{2}}}^{3}{Z_{{1}}}^{3}-15\,{Z_{{2}}}^{4}{Z_{{1}}}^{2}\ \ \ \ \ (34)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\left({Z_{{2}}}^{2}+6\, Z_{{1}}Z_{{2}}+{Z_{{1}}}^{2}\right)\left({Z_{{2}}}^{4}+14\,{Z_{{2}}}^{2}{Z_{{1}}}^{2}+{Z_{{1}}}^{4}\right)\ \ \ \ \ (35)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -x^{6}-y^{6} \ \ \ \ \ (36)$

using 17.

The terms summing to 7 give:

 $\displaystyle N_{7}$ $\displaystyle =$ $\displaystyle {Z_{{1}}}^{7}+47\,{Z_{{1}}}^{3}{Z_{{2}}}^{4}+5\, Z_{{1}}{Z_{{2}}}^{6}+{Z_{{2}}}^{7}+11\,{Z_{{1}}}^{5}{Z_{{2}}}^{2}+47\,{Z_{{1}}}^{4}{Z_{{2}}}^{3}+11\,{Z_{{1}}}^{2}{Z_{{2}}}^{5}+5\,{Z_{{1}}}^{6}Z_{{2}}\ \ \ \ \ (37)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(Z_{{2}}+Z_{{1}}\right)\left({Z_{{2}}}^{6}+4\,{Z_{{2}}}^{5}Z_{{1}}+7\,{Z_{{2}}}^{4}{Z_{{1}}}^{2}+40\,{Z_{{2}}}^{3}{Z_{{1}}}^{3}+7\,{Z_{{1}}}^{4}{Z_{{2}}}^{2}+4\,{Z_{{1}}}^{5}Z_{{2}}+{Z_{{1}}}^{6}\right) \ \ \ \ \ (38)$

We can write the sixth degree polynomial as

 $\displaystyle {Z_{{2}}}^{6}+4\,{Z_{{2}}}^{5}Z_{{1}}+7\,{Z_{{2}}}^{4}{Z_{{1}}}^{2}+40\,{Z_{{2}}}^{3}{Z_{{1}}}^{3}+7\,{Z_{{1}}}^{4}{Z_{{2}}}^{2}+4\,{Z_{{1}}}^{5}Z_{{2}}+{Z_{{1}}}^{6}$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (39)$ $\displaystyle \left(Z_{1}+Z_{2}\right)^{6}-2Z_{1}^{5}Z_{2}-2Z_{1}Z_{2}^{5}-8Z_{1}^{4}Z_{2}^{2}-8Z_{1}^{2}Z_{2}^{4}+20Z_{1}^{3}Z_{2}^{3}$ $\displaystyle =$ $\displaystyle ...\ \ \ \ \ (40)$ $\displaystyle x^{6}-\frac{y^{2}}{2}\left(Z_{1}^{4}+Z_{2}^{4}+4Z_{1}^{3}Z_{2}+4Z_{1}Z_{2}^{3}-10Z_{1}^{2}Z_{2}^{2}\right)$ $\displaystyle =$ $\displaystyle x^{6}-\frac{y^{2}}{2}\left(x^{4}-y^{4}\right)\ \ \ \ \ (41)$ $\displaystyle N_{7}$ $\displaystyle =$ $\displaystyle x\left[x^{6}-\frac{y^{2}}{2}\left(x^{4}-y^{4}\right)\right] \ \ \ \ \ (42)$

Putting it together:

 $\displaystyle I_{1}$ $\displaystyle =$ $\displaystyle -2E_{1}A^{2}\frac{x\left[x^{6}-\frac{y^{2}}{2}\left(x^{4}-y^{4}\right)\right]-x^{6}-y^{6}}{x^{5}}\ \ \ \ \ (43)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -E_{1}x\frac{x\left[x^{6}-\frac{y^{2}}{2}\left(x^{4}-y^{4}\right)\right]-x^{6}-y^{6}}{x^{6}+y^{6}} \ \ \ \ \ (44)$

The other integral is

$\displaystyle I_{2}=\frac{e^{2}}{4\pi\epsilon_{0}}\int\psi\left[\frac{Z_{2}-1}{r_{1}}+\frac{Z_{1}-1}{r_{2}}\right]\psi_{1}\left(r_{2}\right)\psi_{2}\left(r_{1}\right)d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2} \ \ \ \ \ (45)$

which turns out to be equal to ${I_{1}}$. So finally

 $\displaystyle \left\langle H\right\rangle$ $\displaystyle =$ $\displaystyle E_{1}\left(Z_{1}^{2}+Z_{2}^{2}\right)+I_{1}+I_{2}+\left\langle V_{ee}\right\rangle \ \ \ \ \ (46)$ $\displaystyle$ $\displaystyle =$ $\displaystyle E_{1}\left(Z_{1}^{2}+Z_{2}^{2}\right)-2E_{1}x\frac{x\left[x^{6}-\frac{y^{2}}{2}\left(x^{4}-y^{4}\right)\right]-x^{6}-y^{6}}{x^{6}+y^{6}}-\frac{1}{8}{\frac{E_{{1}}x{y}^{2}\left(4\,{x}^{4}+5\,{y}^{4}+{x}^{2}{y}^{2}\right)}{{x}^{6}+{y}^{6}}} \ \ \ \ \ (47)$

Using ${Z_{1}^{2}+Z_{2}^{2}=x^{2}-y^{2}/2}$, we can put everything over a common denominator and cancel terms to get

$\displaystyle \left\langle H\right\rangle =\frac{E_{1}}{x^{6}+y^{6}}\left[-x^{8}+2x^{7}+\frac{1}{2}x^{6}y^{2}-\frac{1}{2}x^{5}y^{2}-\frac{1}{8}x^{3}y^{4}+\frac{11}{8}xy^{6}-\frac{1}{2}y^{8}\right] \ \ \ \ \ (48)$

To find the minimum of this, we can find the maximum of ${\left\langle H\right\rangle /E_{1}}$ (since ${E_{1}<0}$) and we can use Maple’s Maximize function to do this numerically. We find

 $\displaystyle x_{min}$ $\displaystyle =$ $\displaystyle 1.32245\ \ \ \ \ (49)$ $\displaystyle y_{min}$ $\displaystyle =$ $\displaystyle 1.08505\ \ \ \ \ (50)$ $\displaystyle \left\langle H\right\rangle _{min}$ $\displaystyle =$ $\displaystyle 1.0266E_{1} \ \ \ \ \ (51)$

These values of ${x_{min}}$ and ${y_{min}}$ correspond to

 $\displaystyle Z_{1-min}$ $\displaystyle =$ $\displaystyle 1.03923\ \ \ \ \ (52)$ $\displaystyle Z_{2-min}$ $\displaystyle =$ $\displaystyle 0.28322 \ \ \ \ \ (53)$

Thus the upper bound on the energy is slightly lower than ${E_{1}}$ indicating that the ${\mbox{H}^{-}}$ ion is stable.

## Rubber band helium

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 7.17.

We’ve looked at the helium atom using the variational principle. Although the helium atom using the correct Coulomb potential cannot be solved exactly, a variant known as ‘rubber band helium’ can be. Here we use simple harmonic oscillator potentials. The hamiltonian is then:

$\displaystyle H=-\frac{\hbar^{2}}{2m}\left(\nabla_{1}^{2}+\nabla_{2}^{2}\right)+\frac{1}{2}m\omega^{2}\left(r_{1}^{2}+r_{2}^{2}\right)-\frac{\lambda}{4}m\omega^{2}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|^{2} \ \ \ \ \ (1)$

By introducing a change of variables, we can decouple the hamiltonian. Let

 $\displaystyle \mathbf{u}$ $\displaystyle \equiv$ $\displaystyle \frac{1}{\sqrt{2}}\left(\mathbf{r}_{1}+\mathbf{r}_{2}\right)\ \ \ \ \ (2)$ $\displaystyle \mathbf{v}$ $\displaystyle \equiv$ $\displaystyle \frac{1}{\sqrt{2}}\left(\mathbf{r}_{1}-\mathbf{r}_{2}\right) \ \ \ \ \ (3)$

The gradient operators then transform as

 $\displaystyle \nabla_{u}$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{2}}\left(\nabla_{1}+\nabla_{2}\right)\ \ \ \ \ (4)$ $\displaystyle \nabla_{v}$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{2}}\left(\nabla_{1}-\nabla_{2}\right)\ \ \ \ \ (5)$ $\displaystyle \nabla_{u}^{2}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(\nabla_{1}^{2}+\nabla_{2}^{2}+2\nabla_{1}\cdot\nabla_{2}\right)\ \ \ \ \ (6)$ $\displaystyle \nabla_{v}^{2}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(\nabla_{1}^{2}+\nabla_{2}^{2}-2\nabla_{1}\cdot\nabla_{2}\right)\ \ \ \ \ (7)$ $\displaystyle \nabla_{u}^{2}+\nabla_{v}^{2}$ $\displaystyle =$ $\displaystyle \nabla_{1}^{2}+\nabla_{2}^{2} \ \ \ \ \ (8)$

For the potential terms, we have

 $\displaystyle u^{2}+v^{2}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left[r_{1}^{2}+r_{2}^{2}+2\mathbf{r}_{1}\cdot\mathbf{r}_{2}+r_{1}^{2}+r_{2}^{2}-2\mathbf{r}_{1}\cdot\mathbf{r}_{2}\right]\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle r_{1}^{2}+r_{2}^{2}\ \ \ \ \ (10)$ $\displaystyle \left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|^{2}$ $\displaystyle =$ $\displaystyle r_{1}^{2}+r_{2}^{2}-2\mathbf{r}_{1}\cdot\mathbf{r}_{2}\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 2v^{2} \ \ \ \ \ (12)$

Thus the hamiltonian separates:

$\displaystyle H=\left[-\frac{\hbar^{2}}{2m}\nabla_{u}^{2}+\frac{1}{2}m\omega^{2}u^{2}\right]+\left[-\frac{\hbar^{2}}{2m}\nabla_{v}^{2}+\frac{1}{2}m\omega^{2}\left(1-\lambda\right)v^{2}\right] \ \ \ \ \ (13)$

which is the sum of two 3-d harmonic oscillators. The exact ground state energy of this system are then just the sum of the two separate oscillator energies:

$\displaystyle E_{0}=\frac{3}{2}\hbar\omega+\frac{3}{2}\hbar\omega\sqrt{1-\lambda} \ \ \ \ \ (14)$

To test the variational principle for this potential, we can start with the (known) ground state wave function for the 3-d harmonic oscillator as the test function.

$\displaystyle \psi=\left(\frac{m\omega}{\pi\hbar}\right)^{3/2}e^{-m\omega\left(r_{1}^{2}+r_{2}^{2}\right)/2\hbar} \ \ \ \ \ (15)$

This function is an eigenfunction of the first two terms in 1 with energy ${3\hbar\omega}$ so we have

$\displaystyle \left\langle H\right\rangle =3\hbar\omega+\left\langle V_{\lambda}\right\rangle \ \ \ \ \ (16)$

where

 $\displaystyle \left\langle V_{\lambda}\right\rangle$ $\displaystyle =$ $\displaystyle -\frac{\lambda}{4}m\omega^{2}\left(\frac{m\omega}{\pi\hbar}\right)^{3}\int e^{-m\omega\left(r_{1}^{2}+r_{2}^{2}\right)/\hbar}\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|^{2}d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{\lambda}{4}m\omega^{2}\left(\frac{m\omega}{\pi\hbar}\right)^{3}\int e^{-m\omega\left(r_{1}^{2}+r_{2}^{2}\right)/\hbar}\left(r_{1}^{2}+r_{2}^{2}-2r_{1}r_{2}\cos\theta_{2}\right)d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2} \ \ \ \ \ (18)$

The term with ${\cos\theta_{2}}$ integrates to zero when we do the ${\theta_{2}}$ integral, so we’re left with two Gaussian integrals and we get

$\displaystyle \left\langle V_{\lambda}\right\rangle =-\frac{3}{4}\lambda\hbar\omega \ \ \ \ \ (19)$

Plugging this back into 16 we get

$\displaystyle \left\langle H\right\rangle =3\hbar\omega\left(1-\frac{\lambda}{4}\right) \ \ \ \ \ (20)$

This is actually the Taylor expansion with respect to ${\lambda}$ of 14 up to first order.

## Variational principle and the electron in a magnetic field

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 7.16.

As an example of the two-state system we looked at earlier, we can consider an electron in a magnetic field. The hamiltonian for a spin-1/2 particle such as the electron in a magnetic field is

$\displaystyle \mathsf{H}=-\gamma\mathbf{B}\cdot\mathsf{S} \ \ \ \ \ (1)$

where ${\mathsf{S}}$ is the spin matrix and ${\gamma=-e/m}$ is the gyromagnetic ratio of the electron.

We start with a uniform magnetic field in the ${z}$ direction ${\mathbf{B}=B_{z}\hat{\mathbf{z}}}$, for which the hamiltonian is

$\displaystyle \mathsf{H}=\frac{e}{m}B_{z}S_{z}=\frac{eB_{z}\hbar}{2m}\left[\begin{array}{cc} 1 & 0\\ 0 & -1 \end{array}\right] \ \ \ \ \ (2)$

The energies are

$\displaystyle E_{b,a}=\pm\frac{eB_{z}\hbar}{2m} \ \ \ \ \ (3)$

(taking ${E_{b}>E_{a}}$ as before).

We now introduce a perturbation by turning on a weak field in the ${x}$ direction ${\mathbf{B}'=B_{x}\hat{\mathbf{x}}.}$ The perturbation in the hamiltonian is

$\displaystyle H'=\frac{eB_{x}\hbar}{2m}\left[\begin{array}{cc} 0 & 1\\ 1 & 0 \end{array}\right] \ \ \ \ \ (4)$

This has the same form as the perturbation in the previous problem:

$\displaystyle H'=\left[\begin{array}{cc} 0 & h\\ h & 0 \end{array}\right] \ \ \ \ \ (5)$

with

$\displaystyle h=\frac{eB_{x}\hbar}{2m} \ \ \ \ \ (6)$

Using second order perturbation theory, the the perturbations on the two energies are

 $\displaystyle E_{a2}$ $\displaystyle =$ $\displaystyle \frac{h^{2}}{E_{a}-E_{b}}=-\frac{\hbar eB_{x}^{2}}{4mB_{z}}\ \ \ \ \ (7)$ $\displaystyle E_{b2}$ $\displaystyle =$ $\displaystyle \frac{h^{2}}{E_{b}-E_{a}}=\frac{\hbar eB_{x}^{2}}{4mB_{z}} \ \ \ \ \ (8)$

The variational calculation gave the exact answer, which is

 $\displaystyle \left\langle H\right\rangle _{min}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(E_{a}+E_{b}-\sqrt{\left(E_{a}-E_{b}\right)^{2}+4h^{2}}\right)\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar e}{2m}\sqrt{B_{x}^{2}+B_{z}^{2}} \ \ \ \ \ (10)$

That is, it is the same as 3 except we now have the magnitude of the full magnetic field.

## Variational principle with a two-state hamiltonian

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 7.15.

Suppose we have a system with just two possible energies and corresponding eigenstates, which we’ll call ${\psi_{a}}$ with energy ${E_{a}}$ and ${\psi_{b}}$ with energy ${E_{b}}$, with ${\left\langle \left.a\right|b\right\rangle =0}$, ${\left\langle \left.a\right|a\right\rangle =\left\langle \left.b\right|b\right\rangle =1}$ ${E_{a}. Now we turn on a perturbation ${H'}$ which has the matrix elements

$\displaystyle H'=\left[\begin{array}{cc} 0 & h\\ h & 0 \end{array}\right] \ \ \ \ \ (1)$

The total hamiltonian is now ${H=H_{0}+H'}$, where ${H_{0}}$ is the unperturbed hamiltonian. The matrix elements of ${H}$ are then

$\displaystyle H=\left[\begin{array}{cc} E_{a} & h\\ h & E_{b} \end{array}\right] \ \ \ \ \ (2)$

so the exact perturbed energies are the eigenvalues of this matrix, which are

$\displaystyle E^{\prime}=\frac{1}{2}\left(E_{a}+E_{b}\pm\sqrt{\left(E_{a}-E_{b}\right)^{2}+4h^{2}}\right) \ \ \ \ \ (3)$

Now we can apply perturbation theory to this problem. Since the diagonal matrix elements of ${H^{\prime}}$ are both zero, the first order perturbation is also zero. The second order perturbation is

$\displaystyle E_{n2}=\sum_{j\ne n}\frac{\left|\left\langle j0\right|H^{\prime}\left|n0\right\rangle \right|^{2}}{E_{n0}-E_{j0}} \ \ \ \ \ (4)$

This gives for the perturbations on the two energies:

 $\displaystyle E_{a2}$ $\displaystyle =$ $\displaystyle \frac{h^{2}}{E_{a}-E_{b}}\ \ \ \ \ (5)$ $\displaystyle E_{b2}$ $\displaystyle =$ $\displaystyle \frac{h^{2}}{E_{b}-E_{a}} \ \ \ \ \ (6)$

If we expand 3 in a Taylor series, these are second order terms in the series.

Finally, we can use the variational principle with the trial function

$\displaystyle \psi=\left(\cos\phi\right)\psi_{a}+\left(\sin\phi\right)\psi_{b}=\left[\begin{array}{c} \cos\phi\\ \sin\phi \end{array}\right] \ \ \ \ \ (7)$

We can calculate ${\left\langle H\right\rangle }$ as follows:

 $\displaystyle \left\langle H\right\rangle$ $\displaystyle =$ $\displaystyle \psi^{T}H\psi\ \ \ \ \ (8)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left[\begin{array}{cc} \cos\phi & \sin\phi\end{array}\right]\left[\begin{array}{cc} E_{a} & h\\ h & E_{b} \end{array}\right]\left[\begin{array}{c} \cos\phi\\ \sin\phi \end{array}\right]\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle E_{a}\cos^{2}\phi+E_{b}\sin^{2}\phi+2h\sin\phi\cos\phi \ \ \ \ \ (10)$

The variable parameter here is ${\phi}$ so we take the derivative with respect to it and set to zero to get the minimum energy:

 $\displaystyle \frac{d\left\langle H\right\rangle }{d\phi}$ $\displaystyle =$ $\displaystyle \left(E_{b}-E_{a}\right)\left(2\sin\phi\cos\phi\right)+2h\left(\cos^{2}\phi-\sin^{2}\phi\right)\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(E_{b}-E_{a}\right)\sin2\phi+2h\cos2\phi=0\ \ \ \ \ (12)$ $\displaystyle \tan2\phi_{min}$ $\displaystyle =$ $\displaystyle -\frac{2h}{E_{b}-E_{a}}\ \ \ \ \ (13)$ $\displaystyle \sin2\phi_{min}$ $\displaystyle =$ $\displaystyle -\frac{2h}{E_{b}-E_{a}}\cos2\phi_{min}\ \ \ \ \ (14)$ $\displaystyle \cos2\phi_{min}$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{1+\tan^{2}2\phi_{min}}}\ \ \ \ \ (15)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left(1+\frac{4h^{2}}{\left(E_{b}-E_{a}\right)^{2}}\right)^{-1/2}\ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{E_{b}-E_{a}}{\sqrt{\left(E_{b}-E_{a}\right)^{2}+4h^{2}}} \ \ \ \ \ (17)$

We can express ${\left\langle H\right\rangle _{min}}$ using trig identities:

 $\displaystyle \left\langle H\right\rangle _{min}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(1+\cos2\phi_{min}\right)E_{a}+\frac{1}{2}\left(1-\cos2\phi_{min}\right)E_{b}+h\sin2\phi_{min}\ \ \ \ \ (18)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left(E_{a}+E_{b}-\sqrt{\left(E_{a}-E_{b}\right)^{2}+4h^{2}}\right) \ \ \ \ \ (19)$

which is exactly the lower of the two exact energies 3. The variational principle gives the exact answer because the trial function is the exact eigenfunction, with ${\phi_{min}}$ giving the components of the two unperturbed eigenfunctions that make up the perturbed eigenfunction.

## Variational principle and the Yukawa potential

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 7.14.

A result from quantum field theory is the Yukawa potential which, when applied to the Coulomb potential, gives

$\displaystyle V=-\frac{e^{2}}{4\pi\epsilon_{0}}\frac{e^{-\mu r}}{r} \ \ \ \ \ (1)$

where ${\mu=m_{\gamma}c/\hbar}$ and ${m_{\gamma}}$ is the mass of the particle mediating the force, which in the case of electromagnetism is the photon. Thus if the photon did have a non-zero mass, this exponential factor would appear in the potential.

We can use the variational principle to estimate the effect of this potential on the ground state of the hydrogen atom. We can try as a trial function

$\displaystyle \psi=Ae^{-br/a} \ \ \ \ \ (2)$

where ${a}$ is the Bohr radius and ${b}$ is the parameter we can vary to find the minimum upper bound on the energy. Most of the calculations from here onwards were done using Maple, as they get a bit messy in places.

First, we normalize it:

 $\displaystyle A^{2}\int e^{-2br/a}r^{2}\sin\theta d\phi d\theta dr$ $\displaystyle =$ $\displaystyle 1\ \ \ \ \ (3)$ $\displaystyle A$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{\pi}}\left(\frac{b}{a}\right)^{3/2} \ \ \ \ \ (4)$

We now find ${H\psi=\left(T+V\right)\psi}$:

 $\displaystyle T\psi$ $\displaystyle =$ $\displaystyle -\frac{\hbar^{2}}{2m}\nabla^{2}\psi\ \ \ \ \ (5)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}\left(2a-br\right)\sqrt{a}b^{5/2}e^{-br/a}}{2\pi^{2}\epsilon_{0}a^{2}r}\ \ \ \ \ (6)$ $\displaystyle V\psi$ $\displaystyle =$ $\displaystyle -\frac{e^{2}}{4\pi^{3/2}\epsilon_{0}}\frac{e^{-\mu r}}{r}\left(\frac{b}{a}\right)^{3/2}e^{-br/a} \ \ \ \ \ (7)$

We can now find ${\left\langle H\right\rangle }$:

 $\displaystyle \left\langle T\right\rangle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}\sqrt{a}b^{5/2}}{2\pi^{2}\epsilon_{0}a^{2}}\int\left(2a-br\right)e^{-br/a}r\sin\theta d\phi d\theta dr\ \ \ \ \ (8)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}b^{2}}{2a^{2}m}\ \ \ \ \ (9)$ $\displaystyle \left\langle V\right\rangle$ $\displaystyle =$ $\displaystyle -\frac{e^{2}}{4\pi^{3/2}\epsilon_{0}}\left(\frac{b}{a}\right)^{3/2}\int e^{-\mu r}e^{-br/a}r\sin\theta d\phi d\theta dr\ \ \ \ \ (10)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\frac{e^{2}b^{3}}{\pi\epsilon_{0}a\left(4b^{2}+4b\mu a+\mu^{2}a^{2}\right)}\ \ \ \ \ (11)$ $\displaystyle \left\langle H\right\rangle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}b^{2}}{2a^{2}m}-\frac{e^{2}b^{3}}{\pi\epsilon_{0}a\left(4b^{2}+4b\mu a+\mu^{2}a^{2}\right)} \ \ \ \ \ (12)$

At this stage, we’d take the derivative w.r.t. ${b}$ and set to zero to find ${\left\langle H\right\rangle _{min}}$. As you can see, the derivative is quite complicated, so we’ll do this by software. After substituting ${x\equiv\mu a}$ we get

$\displaystyle \frac{d\left\langle H\right\rangle }{db}=\frac{b\hbar^{2}}{ma^{2}}-\frac{3b^{2}e^{2}}{\pi\epsilon_{0}a\left(4bx+x^{2}+4b^{2}\right)}+\frac{b^{3}e^{2}\left(8b+4x\right)}{\pi\epsilon_{0}a\left(4bx+x^{2}+4b^{2}\right)^{2}} \ \ \ \ \ (13)$

For ${\mu=0}$, ${b=1}$ since that just gives us the exact hydrogen wave function for the normal Coulomb potential. Since ${x=\mu a\ll1}$, we’d expect ${b}$ to be very close to 1, so we can write ${b=1+y}$ where ${y\ll1}$. Therefore we can try expanding the expression in a Taylor series, first in ${x}$ about ${x=0}$ and then in ${y}$ about ${y=0}$, and saving only terms that are second order in products of ${x}$ and ${y}$. First, we expand around ${x=0}$:

$\displaystyle \frac{d\left\langle H\right\rangle }{db}=\frac{b\hbar^{2}}{ma^{2}}-\frac{e^{2}}{4\pi\epsilon_{0}a}+\frac{3e^{2}}{16\pi\epsilon_{0}ab^{2}}x^{2}+\mathcal{O}\left(x^{3}\right) \ \ \ \ \ (14)$

Now we substitute ${b=1+y}$ and expand about ${y=0}$, saving only terms where the combined powers of ${x}$ and ${y}$ are ${\le2}$:

$\displaystyle \frac{d\left\langle H\right\rangle }{db}\approx\frac{\hbar^{2}}{ma^{2}}-\frac{e^{2}}{4\pi\epsilon_{0}a}+\frac{\hbar^{2}}{ma^{2}}y+\frac{3e^{2}}{16\pi\epsilon_{0}a}x^{2} \ \ \ \ \ (15)$

Now we can set this to zero and solve for ${y}$ (since we’re trying to find ${b_{min}}$):

$\displaystyle b_{min}=1+y_{min}=1+\frac{4me^{2}a-16\pi\hbar^{2}\epsilon_{0}-3e^{2}max^{2}}{16\pi\epsilon_{0}\hbar^{2}} \ \ \ \ \ (16)$

We can substitute this back into 12 and then take the Taylor expansion of the result around ${x=0}$. The result is:

$\displaystyle \left\langle H\right\rangle _{min}=-\frac{me^{4}}{32\pi^{2}\hbar^{2}\epsilon_{0}^{2}}+\frac{e^{2}}{4\pi\epsilon_{0}a}x-\frac{3\hbar^{2}}{4ma^{2}}x^{2}+\mathcal{O}\left(x^{3}\right) \ \ \ \ \ (17)$

Substituting in for the Bohr radius:

$\displaystyle a=\frac{4\pi\epsilon_{0}\hbar^{2}}{me^{2}} \ \ \ \ \ (18)$

we get

 $\displaystyle \left\langle H\right\rangle _{min}$ $\displaystyle =$ $\displaystyle -\frac{me^{4}}{32\pi^{2}\hbar^{2}\epsilon_{0}^{2}}\left(1-2x+\frac{3}{2}x^{2}+\mathcal{O}\left(x^{3}\right)\right)\ \ \ \ \ (19)$ $\displaystyle$ $\displaystyle =$ $\displaystyle E_{1}\left(1-2x+\frac{3}{2}x^{2}+\mathcal{O}\left(x^{3}\right)\right) \ \ \ \ \ (20)$

where ${E_{1}}$ is the exact ground state energy for hydrogen.

## Variational principle and the hydrogen atom

References: Griffiths, David J. (2005), Introduction to Quantum Mechanics, 2nd Edition; Pearson Education – Problem 7.13.

Another example of the variational principle this time to estimate the ground state of the hydrogen atom. This is the first example we’ve done in three dimensions, although since it’s spherically symmetric, the integral isn’t much more complicated.

The trial function here is

$\displaystyle \psi=Ae^{-br^{2}} \ \ \ \ \ (1)$

First, we normalize it:

 $\displaystyle A^{2}\int e^{-2br^{2}}r^{2}\sin\theta d\phi d\theta dr$ $\displaystyle =$ $\displaystyle 1\ \ \ \ \ (2)$ $\displaystyle A$ $\displaystyle =$ $\displaystyle \left(\frac{2b}{\pi}\right)^{3/4} \ \ \ \ \ (3)$

We now find ${H\psi}$:

 $\displaystyle H\psi$ $\displaystyle =$ $\displaystyle -\frac{\hbar^{2}}{2m}\nabla^{2}\psi-\frac{e^{2}}{4\pi\epsilon_{0}r}\psi\ \ \ \ \ (4)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\left(\frac{2b}{\pi}\right)^{3/4}e^{-br^{2}}\left[\frac{2b^{2}r^{2}-3b}{m}+\frac{e^{2}}{4\pi\epsilon_{0}r}\right] \ \ \ \ \ (5)$

We can now find ${\left\langle H\right\rangle }$:

 $\displaystyle \left\langle H\right\rangle$ $\displaystyle =$ $\displaystyle -\left(\frac{2b}{\pi}\right)^{3/2}\int e^{-2br^{2}}\left[\frac{2b^{2}r^{2}-3b}{m}+\frac{e^{2}}{4\pi\epsilon_{0}r}\right]r^{2}\sin\theta d\phi d\theta dr\ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3\hbar^{2}}{2m}b-\frac{e^{2}}{\sqrt{2}\pi^{3/2}\epsilon_{0}}\sqrt{b} \ \ \ \ \ (7)$

Taking the derivative w.r.t. ${b}$ and setting to zero gives

$\displaystyle b_{min}=\frac{m^{2}e^{4}}{18\hbar^{4}\pi^{3}\epsilon_{0}^{2}} \ \ \ \ \ (8)$

Substituting back into 7 gives

 $\displaystyle \left\langle H\right\rangle _{min}$ $\displaystyle =$ $\displaystyle -\frac{me^{4}}{12\hbar^{2}\pi^{3}\epsilon_{0}^{2}}\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -11.6\mbox{ eV} \ \ \ \ \ (10)$

This is significantly above the correct value of ${-13.6\mbox{ eV}}$.