# Every attractive 1-dimensional potential has a bound state

Reference: Lecture by Barton Zwiebach in MIT course 8.05.1x, week 1, Deep Dive 2 (not in the PDF notes).

Shankar, R. (1994), Principles of Quantum Mechanics, Plenum Press. Section 5.2, Exercise 5.2.2b.

An interesting application of the variational principle in quantum mechances is the following theorem:

Theorem Every 1-dimensional attractive potential has at least one bound state.

To prove this, we need first to define what we mean by an attractive potential ${V\left(x\right)}$. ${V\left(x\right)}$ must satisfy the following conditions:

• ${V\left(x\right)\rightarrow0}$ as ${x\rightarrow\pm\infty}$.
• ${V\left(x\right)<0}$ everywhere.
• ${V\left(x\right)}$ is piecewise continuous. This means that it may have a finite number of jump discontinuities.

One possible form for ${V\left(x\right)}$ is as shown:

This is a particularly simple potential that satisfies the above conditions. We could introduce a few step functions, multiple local maxima and minima, and so on, provided we don’t violate any of the 3 conditions above.

Since ${V\left(x\right)<0}$ everywhere, we can write it as

$\displaystyle V\left(x\right)=-\left|V\left(x\right)\right| \ \ \ \ \ (1)$

What we would like to prove is that for any hamiltonian of the form

$\displaystyle H=-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dx^{2}}-\left|V\left(x\right)\right| \ \ \ \ \ (2)$

the ground state ${E_{0}}$ is a bound state, that is

$\displaystyle E_{0}<0 \ \ \ \ \ (3)$

We can apply the variational principle, which states

If ${\psi}$ is any normalized function and ${H}$ is a hamiltonian, then the ground state energy ${E_{0}}$ of this hamiltonian has an upper bound given by

$\displaystyle E_{0}\le\left\langle \psi\left|H\right|\psi\right\rangle \equiv\left\langle H\right\rangle \ \ \ \ \ (4)$

The use of the variational principle to prove the above theorem involves a bit of a convoluted argument, but the mathematics involved is fairly simple. Our goal is to find some wave function ${\psi_{\alpha}}$ (where ${\alpha}$ is some parameter that we can vary) so that

$\displaystyle E_{0}\le\left\langle \psi_{\alpha}\left|H\right|\psi_{\alpha}\right\rangle =\left\langle H\right\rangle _{\psi_{\alpha}}<0 \ \ \ \ \ (5)$

From 2 we have

 $\displaystyle \left\langle \hat{H}\right\rangle _{\psi_{\alpha}}$ $\displaystyle =$ $\displaystyle \int dx\;\psi_{\alpha}\left(x\right)\hat{H}\psi_{\alpha}\left(x\right)\ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left\langle T\right\rangle _{\psi_{\alpha}}-\left\langle \left|V\left(x\right)\right|\right\rangle _{\psi_{\alpha}} \ \ \ \ \ (7)$

where

 $\displaystyle \left\langle T\right\rangle _{\psi_{\alpha}}$ $\displaystyle =$ $\displaystyle -\int dx\;\psi_{\alpha}\left(x\right)\frac{\hbar^{2}}{2m}\frac{d^{2}}{dx^{2}}\psi_{\alpha}\left(x\right)\ \ \ \ \ (8)$ $\displaystyle \left\langle \left|V\left(x\right)\right|\right\rangle _{\psi_{\alpha}}$ $\displaystyle =$ $\displaystyle \int dx\;\psi_{\alpha}\left(x\right)\left|V\left(x\right)\right|\psi_{\alpha}\left(x\right) \ \ \ \ \ (9)$

We can integrate 8 by parts once to get

 $\displaystyle \left\langle T\right\rangle _{\psi_{\alpha}}$ $\displaystyle =$ $\displaystyle -\int dx\;\psi_{\alpha}\left(x\right)\frac{\hbar^{2}}{2m}\frac{d^{2}}{dx^{2}}\psi_{\alpha}\left(x\right)\ \ \ \ \ (10)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -\left.\frac{\hbar^{2}}{2m}\psi_{\alpha}\left(x\right)\frac{d}{dx}\psi_{\alpha}\left(x\right)\right|_{-\infty}^{\infty}+\frac{\hbar^{2}}{2m}\int dx\;\left(\frac{d}{dx}\psi_{\alpha}\left(x\right)\right)^{2}\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}}{2m}\int_{-\infty}^{\infty}dx\;\left(\frac{d}{dx}\psi_{\alpha}\left(x\right)\right)^{2} \ \ \ \ \ (12)$

where we invoke the usual requirement that ${\psi_{\alpha}}$ and its first derivative vanish at infinity.

We therefore see that since the integrand in the last line is always positive (we’re assuming that ${\psi_{\alpha}}$ is not zero everywhere), that ${\left\langle T\right\rangle _{\psi_{\alpha}}>0}$. Likewise, from 9, ${\left\langle \left|V\left(x\right)\right|\right\rangle _{\psi_{\alpha}}>0}$. Thus in order that ${\left\langle H\right\rangle _{\psi_{\alpha}}<0}$, we must have

$\displaystyle \left\langle T\right\rangle _{\psi_{\alpha}}<\left\langle \left|V\left(x\right)\right|\right\rangle _{\psi_{\alpha}} \ \ \ \ \ (13)$

To get any further, we need to choose a test function ${\psi_{\alpha}\left(x\right)}$. We’ll pick (because it works!)

$\displaystyle \psi_{\alpha}=\left(\frac{\alpha}{\pi}\right)^{1/4}e^{-\frac{1}{2}\alpha x^{2}} \ \ \ \ \ (14)$

The factor of ${\left(\frac{\alpha}{\pi}\right)^{1/4}}$ is required so that ${\psi_{\alpha}}$ is normalized. The integral in 12 can be done using standard methods; I’ll just use Maple, and we find

$\displaystyle \left\langle T\right\rangle _{\psi_{\alpha}}=\frac{\hbar^{2}\alpha}{4m} \ \ \ \ \ (15)$

The integral 9 of course can’t be done exactly if we don’t know what ${V}$ is, so we have just

$\displaystyle \left\langle \left|V\left(x\right)\right|\right\rangle _{\psi_{\alpha}}=\int dx\;\psi_{\alpha}^{2}\left(x\right)\left|V\left(x\right)\right| \ \ \ \ \ (16)$

(No need for modulus signs around ${\psi_{\alpha}}$ since the function 14 is real.) To progress further, we need to start invoking some inequalities to get where we want to go. The argument consists of several steps, so watch carefully as we go along.

From 13 through 15 we have to show that we can satisfy the condition

$\displaystyle \frac{\left\langle \left|V\left(x\right)\right|\right\rangle _{\psi_{\alpha}}}{\left\langle T\right\rangle _{\psi_{\alpha}}}=\frac{4m}{\hbar^{2}\sqrt{\pi}}\frac{1}{\sqrt{\alpha}}\int_{-\infty}^{\infty}e^{-\alpha x^{2}}\left|V\left(x\right)\right|dx>1 \ \ \ \ \ (17)$

Since ${V}$ is arbitrary subject to the 3 conditions above, the only thing we can legitimately fiddle with is the value of ${\alpha}$. We can see that if we choose ${\alpha}$ small enough, we should be able to satisfy this inequality, since for small ${\alpha}$, the ${1/\sqrt{\alpha}}$ term gets large, while the ${e^{-\alpha x^{2}}}$ term in the integrand is bounded between 0 and 1. We need to find some upper limit for ${\alpha}$.

In what follows, you’ll need to refer to the following diagram:

First, we choose some point ${x_{0}}$ at which ${V\left(x_{0}\right)}$ is continuous (that is, we ensure that ${x_{0}}$ isn’t at one of the points where ${V\left(x\right)}$ has a discontinuity, or jump). The value of ${V\left(x_{0}\right)}$ is defined as ${-2v_{0}}$ where ${v_{0}>0}$. Because ${V\rightarrow0}$ at ${x\rightarrow\pm\infty}$, there must be points ${x_{1}}$ and ${x_{2}}$ on either side of ${x_{0}}$ where ${V}$ has the value ${-v_{0}}$ (actually, I’m not sure this is strictly true, because, as ${V}$ is allowed a few jumps, it might jump over the point where it’s equal to ${-v_{0}}$. However, as the number of jumps is required to be finite, there must be some points ${x_{1}}$ and ${x_{2}}$ on either side of ${x_{0}}$ where ${V}$ attains a value that is between ${-2v_{0}}$ and 0, and I think the argument below still works if we choose those points instead.)

Now for the first inequality. We know that, because the integrand is positive

$\displaystyle \int_{-\infty}^{\infty}e^{-\alpha x^{2}}\left|V\left(x\right)\right|dx>\int_{x_{1}}^{x_{2}}e^{-\alpha x^{2}}\left|V\left(x\right)\right|dx \ \ \ \ \ (18)$

Second inequality: in the interval ${x_{1}}$ to ${x_{2}}$, ${\left|V\left(x\right)\right|>v_{0}}$ (see the diagram!), so we have

$\displaystyle \int_{x_{1}}^{x_{2}}e^{-\alpha x^{2}}\left|V\left(x\right)\right|dx>v_{0}\int_{x_{1}}^{x_{2}}e^{-\alpha x^{2}}dx \ \ \ \ \ (19)$

The last integral has no closed form solution, but we know that in the interval ${x_{1}}$ to ${x_{2}}$

$\displaystyle e^{-\alpha x^{2}}>e^{-\alpha\max\left(x_{1}^{2},x_{2}^{2}\right)} \ \ \ \ \ (20)$

Therefore

 $\displaystyle v_{0}\int_{x_{1}}^{x_{2}}e^{-\alpha x^{2}}dx$ $\displaystyle >$ $\displaystyle v_{0}\int_{x_{1}}^{x_{2}}e^{-\alpha\max\left(x_{1}^{2},x_{2}^{2}\right)}dx\ \ \ \ \ (21)$ $\displaystyle$ $\displaystyle =$ $\displaystyle v_{0}\left(x_{2}-x_{1}\right)e^{-\alpha\max\left(x_{1}^{2},x_{2}^{2}\right)} \ \ \ \ \ (22)$

Now suppose we choose ${\alpha}$ to be

$\displaystyle \alpha<\frac{1}{\max\left(x_{1}^{2},x_{2}^{2}\right)} \ \ \ \ \ (23)$

Then

$\displaystyle e^{-\alpha\max\left(x_{1}^{2},x_{2}^{2}\right)}>e^{-1} \ \ \ \ \ (24)$

We can now summarize as follows:

$\displaystyle \int_{-\infty}^{\infty}e^{-\alpha x^{2}}\left|V\left(x\right)\right|dx>v_{0}\left(x_{2}-x_{1}\right)e^{-1} \ \ \ \ \ (25)$

provided we choose ${\alpha}$ according to 23. Plugging this back into 17 we have

$\displaystyle \frac{\left\langle \left|V\left(x\right)\right|\right\rangle _{\psi_{\alpha}}}{\left\langle T\right\rangle _{\psi_{\alpha}}}>\frac{4m}{\hbar^{2}\sqrt{\pi}}\frac{v_{0}\left(x_{2}-x_{1}\right)}{e}\frac{1}{\sqrt{\alpha}} \ \ \ \ \ (26)$

This expression will now be greater than 1 provided that

 $\displaystyle \sqrt{\alpha}$ $\displaystyle <$ $\displaystyle \frac{4m}{\hbar^{2}\sqrt{\pi}}\frac{v_{0}\left(x_{2}-x_{1}\right)}{e}\ \ \ \ \ (27)$ $\displaystyle \alpha$ $\displaystyle <$ $\displaystyle \left[\frac{4m}{\hbar^{2}\sqrt{\pi}}\frac{v_{0}\left(x_{2}-x_{1}\right)}{e}\right]^{2} \ \ \ \ \ (28)$

Comparing 23 and 28, we see that we can satisfy both conditions if we take

$\displaystyle \alpha<\min\left\{ \frac{1}{\max\left(x_{1}^{2},x_{2}^{2}\right)},\left[\frac{4m}{\hbar^{2}\sqrt{\pi}}\frac{v_{0}\left(x_{2}-x_{1}\right)}{e}\right]^{2}\right\} \ \ \ \ \ (29)$

This condition depends on ${x_{1}}$ and ${x_{2}}$ but that doesn’t matter, since both quantities in the RHS of 29 are positive, so there is always some positive value of ${\alpha}$ that satisfies the condition. In other words, going right back to 17 and then to 7, we can always find a value of ${\alpha}$ so that ${\left\langle H\right\rangle <0}$ which means that the ground state of ${H}$ must be negative, which makes it a bound state.

# 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}} \ \ \ \ \ (1)$

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} \ \ \ \ \ (2)$

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 \ldots\ \ \ \ \ (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 \ldots\ \ \ \ \ (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 \ldots\ \ \ \ \ (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 \ldots\ \ \ \ \ (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 \ldots\ \ \ \ \ (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 \ldots\ \ \ \ \ (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 \ldots\ \ \ \ \ (26)$ $\displaystyle x^{4}+{Z_{{2}}}^{3}Z_{{1}}+{Z_{{1}}}^{3}Z_{{2}}+22\,{Z_{{2}}}^{2}{Z_{{1}}}^{2}$ $\displaystyle =$ $\displaystyle \ldots\ \ \ \ \ (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 \ldots\ \ \ \ \ (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 \ldots\ \ \ \ \ (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 – Chapter 7, Post 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 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{\hbar^{2}\left(2b^{2}r^{2}-3b\right)}{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{\hbar^{2}\left(2b^{2}r^{2}-3b\right)}{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}}$.

# Variational principle and harmonic oscillator: a more general trial function

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

In an earlier problem we used the variational principle to estimate the ground state of the harmonic oscillator. The trial function there was

$\displaystyle \psi=\frac{A}{x^{2}+b^{2}} \ \ \ \ \ (1)$

We can generalize this by introducing another parameter ${n}$:

$\displaystyle \psi=\frac{A}{\left(x^{2}+b^{2}\right)^{n}} \ \ \ \ \ (2)$

As usual, we first normalize ${\psi}$:

$\displaystyle A^{2}\int_{-\infty}^{\infty}\frac{dx}{\left(x^{2}+b^{2}\right)^{2n}}=1 \ \ \ \ \ (3)$

As far as I know, there is no simple version of this integral, so we can use tables or Maple to work it out:

$\displaystyle \int_{-\infty}^{\infty}\frac{dx}{\left(x^{2}+b^{2}\right)^{2n}}=\frac{1}{b^{4n}}\frac{\sqrt{\pi}b\Gamma\left(2n-\frac{1}{2}\right)}{\Gamma\left(2n\right)} \ \ \ \ \ (4)$

where ${\Gamma\left(x\right)}$ is the gamma function. Therefore

$\displaystyle A=b^{2n}\left[\frac{\Gamma\left(2n\right)}{\sqrt{\pi}b\Gamma\left(2n-\frac{1}{2}\right)}\right]^{1/2} \ \ \ \ \ (5)$

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

 $\displaystyle \left\langle H\right\rangle$ $\displaystyle =$ $\displaystyle \left\langle \psi\left|H\right|\psi\right\rangle =\left\langle \psi\left|T+V\right|\psi\right\rangle \ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle A^{2}\int_{-\infty}^{\infty}\left[-\frac{\hbar^{2}}{2m}\frac{1}{\left(x^{2}+b^{2}\right)^{n}}\frac{d^{2}}{dx^{2}}\left(\frac{1}{\left(x^{2}+b^{2}\right)^{n}}\right)+\frac{m\omega^{2}x^{2}}{2\left(x^{2}+b^{2}\right)^{2n}}\right]dx\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}\left(16n^{3}-16n^{2}+3n\right)+b^{4}m^{2}\omega^{2}\left(4n+2\right)}{4mb^{2}\left(2n+1\right)\left(4n-3\right)} \ \ \ \ \ (8)$

where Maple was used to do the integrals and simplify the result.

We now take the derivative w.r.t. ${b}$ and set to zero to find ${\left\langle H\right\rangle _{min}}$:

 $\displaystyle b_{min}$ $\displaystyle =$ $\displaystyle \left[\frac{n\left(16n^{2}-16n+3\right)}{2\left(2n+1\right)}\right]^{1/4}\sqrt{\frac{\hbar}{m\omega}}\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left[\frac{n\left(4n-1\right)\left(4n-3\right)}{2\left(2n+1\right)}\right]^{1/4}\sqrt{\frac{\hbar}{m\omega}} \ \ \ \ \ (10)$

This gives an upper bound of

$\displaystyle \left\langle H\right\rangle _{min}=\sqrt{\frac{n\left(4n-1\right)}{2\left(2n+1\right)\left(4n-3\right)}}\hbar\omega \ \ \ \ \ (11)$

For ${n=1}$ this reduces to the solution we had earlier:

$\displaystyle \left\langle H\right\rangle _{n=1}=\frac{1}{\sqrt{2}}\hbar\omega \ \ \ \ \ (12)$

Also, as ${n\rightarrow\infty}$, this tends to the exact answer:

$\displaystyle \lim_{n\rightarrow\infty}\left\langle H\right\rangle =\frac{1}{2}\hbar\omega \ \ \ \ \ (13)$

We can use the corollary to estimate the first excited state’s energy. Since we know the exact ground state wave function ${\psi_{0}}$ of the harmonic oscillator is even (it’s a Gaussian), we can take as a trial function the odd function:

$\displaystyle \psi=\frac{Bx}{\left(x^{2}+b^{2}\right)^{n}} \ \ \ \ \ (14)$

Following the same procedure as above, we get for ${B}$:

 $\displaystyle B^{2}\int_{-\infty}^{\infty}\frac{x^{2}dx}{\left(x^{2}+b^{2}\right)^{2n}}$ $\displaystyle =$ $\displaystyle 1\ \ \ \ \ (15)$ $\displaystyle B$ $\displaystyle =$ $\displaystyle b^{2n}\left[\frac{2\Gamma\left(2n\right)}{\sqrt{\pi}b^{3}\Gamma\left(2n-\frac{3}{2}\right)}\right]^{1/2} \ \ \ \ \ (16)$

For the energy, we get

 $\displaystyle \left\langle H\right\rangle$ $\displaystyle =$ $\displaystyle \left\langle \psi\left|H\right|\psi\right\rangle =\left\langle \psi\left|T+V\right|\psi\right\rangle \ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle B^{2}\int_{-\infty}^{\infty}\left[-\frac{\hbar^{2}}{2m}\frac{x}{\left(x^{2}+b^{2}\right)^{n}}\frac{d^{2}}{dx^{2}}\left(\frac{x}{\left(x^{2}+b^{2}\right)^{n}}\right)+\frac{m\omega^{2}x^{4}}{2\left(x^{2}+b^{2}\right)^{2n}}\right]dx\ \ \ \ \ (18)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 3\frac{\hbar^{2}\left(16n^{3}-32n^{2}+15n\right)+b^{4}m^{2}\omega^{2}\left(4n+2\right)}{4mb^{2}\left(2n+1\right)\left(4n-5\right)} \ \ \ \ \ (19)$

Finding the ${b}$ that minimizes ${\left\langle H\right\rangle }$ gives

 $\displaystyle b_{min}$ $\displaystyle =$ $\displaystyle \left[\frac{n\left(16n^{2}-32n+15\right)}{2\left(2n+1\right)}\right]^{1/4}\sqrt{\frac{\hbar}{m\omega}}\ \ \ \ \ (20)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \left[\frac{n\left(4n-5\right)\left(4n-3\right)}{2\left(2n+1\right)}\right]^{1/4}\sqrt{\frac{\hbar}{m\omega}}\ \ \ \ \ (21)$ $\displaystyle \left\langle H\right\rangle _{min}$ $\displaystyle =$ $\displaystyle \sqrt{\frac{n\left(4n-3\right)}{2\left(2n+1\right)\left(4n-5\right)}}\hbar\omega \ \ \ \ \ (22)$

Again, for large ${n}$ we tend to the exact answer:

$\displaystyle \lim_{n\rightarrow\infty}\left\langle H\right\rangle _{min}=\frac{3}{2}\hbar\omega \ \ \ \ \ (23)$

To see why the limit of large ${n}$ gives the exact answer, we can use Maple’s limit function to find the limit of the trial functions for large ${n}$. We find (remembering to substitute for ${A}$ and ${B}$ by their expressions from above):

 $\displaystyle \lim_{n\rightarrow\infty}\frac{A}{\left(x^{2}+b^{2}\right)^{n}}$ $\displaystyle =$ $\displaystyle \left(\frac{m\omega}{\pi\hbar}\right)^{1/4}e^{-m\omega x^{2}/2\hbar}=\psi_{0}\ \ \ \ \ (24)$ $\displaystyle \lim_{n\rightarrow\infty}\frac{Bx}{\left(x^{2}+b^{2}\right)^{n}}$ $\displaystyle =$ $\displaystyle \left(\frac{m\omega}{\pi\hbar}\right)^{1/4}\sqrt{\frac{2m\omega}{\hbar}}xe^{-m\omega x^{2}/2\hbar}=\psi_{1} \ \ \ \ \ (25)$

That is, in the limit of large ${n}$, both trial functions tend to the exact wave functions.

# Variational principle and the harmonic oscillator – 2

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

Here’s an estimate of the first two energy levels of the harmonic oscillator using the variational principle. Our trial function is

$\displaystyle \psi=\begin{cases} A\cos\frac{\pi x}{a} & -\frac{a}{2}\le x\le\frac{a}{2}\\ 0 & \mbox{otherwise} \end{cases} \ \ \ \ \ (1)$

First, we normalize the wave function:

 $\displaystyle 1$ $\displaystyle =$ $\displaystyle \left|A\right|^{2}\int_{-a/2}^{a/2}\cos^{2}\frac{\pi x}{a}dx\ \ \ \ \ (2)$ $\displaystyle A$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2}{a}} \ \ \ \ \ (3)$

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

$\displaystyle \left\langle H\right\rangle =\left\langle \psi\left|H\right|\psi\right\rangle =\left\langle \psi\left|T+V\right|\psi\right\rangle \ \ \ \ \ (4)$

The kinetic energy term ${T}$ contains the second derivative of ${\psi}$ which contains delta functions at ${x=\pm\frac{a}{2}}$. However, since ${\psi}$ itself is zero at these two points, we don’t need to worry about these points, so we can just ignore the delta functions in what follows.

For the kinetic energy operator ${T}$ we have

 $\displaystyle T\left|\psi\right\rangle$ $\displaystyle =$ $\displaystyle -\frac{\hbar^{2}}{2m}\frac{d^{2}\psi}{dx^{2}}\ \ \ \ \ (5)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}\pi^{2}}{\sqrt{2}ma^{3/2}}\cos\frac{\pi x}{a} \ \ \ \ \ (6)$

The potential energy is

 $\displaystyle V\left|\psi\right\rangle$ $\displaystyle =$ $\displaystyle \frac{1}{2}m\omega^{2}x^{2}\psi\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{m\omega^{2}x^{2}}{\sqrt{2a}}\cos\frac{\pi x}{a} \ \ \ \ \ (8)$

Combining them we get

 $\displaystyle \left\langle H\right\rangle$ $\displaystyle =$ $\displaystyle \sqrt{\frac{2}{a}}\int_{-a/2}^{a/2}\left(\frac{\hbar^{2}\pi^{2}}{\sqrt{2}ma^{3/2}}\cos^{2}\frac{\pi x}{a}+\frac{m\omega^{2}x^{2}}{\sqrt{2a}}\cos^{2}\frac{\pi x}{a}\right)dx\ \ \ \ \ (9)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\left(\pi^{2}-6\right)m^{2}\omega^{2}a^{4}+12\pi^{4}\hbar^{2}}{24\pi^{2}ma^{2}} \ \ \ \ \ (10)$

To find the value of ${a}$ that minimizes ${\left\langle H\right\rangle }$ we take the derivative and set to zero as usual, and find that

$\displaystyle a_{min}=\frac{3^{1/4}\pi}{\left(\pi^{2}-6\right)^{1/4}}\sqrt{\frac{2\hbar}{m\omega}} \ \ \ \ \ (11)$

with the corresponding energy of

 $\displaystyle \left\langle H\right\rangle _{min}$ $\displaystyle =$ $\displaystyle \frac{\sqrt{3\left(\pi^{2}-6\right)}}{6}\hbar\omega\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle \approx$ $\displaystyle 0.569\hbar\omega\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle >$ $\displaystyle \frac{1}{2}\hbar\omega \ \ \ \ \ (14)$

The estimate here is quite close to the exact value.

We can now use the corollary to the variational principle to get an estimate of the first excited state. To do this, we need to find a trial function that is orthogonal to the true ground state wave function. We know that this is

$\displaystyle \psi_{0}=\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}e^{-m\omega x^{2}/2\hbar} \ \ \ \ \ (15)$

so, since it’s an even function, any odd function is orthogonal to it. We can therefore use

$\displaystyle \psi=B\sin\frac{\pi x}{a} \ \ \ \ \ (16)$

Following the same procedure as above, we first normalize ${\psi}$:

 $\displaystyle 1$ $\displaystyle =$ $\displaystyle \left|B\right|^{2}\int_{-a/2}^{a/2}\sin^{2}\frac{\pi x}{a}dx\ \ \ \ \ (17)$ $\displaystyle B$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{a}} \ \ \ \ \ (18)$

For the kinetic energy operator ${T}$ we have

 $\displaystyle T\left|\psi\right\rangle$ $\displaystyle =$ $\displaystyle -\frac{\hbar^{2}}{2m}\frac{d^{2}\psi}{dx^{2}}\ \ \ \ \ (19)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\hbar^{2}\pi^{2}}{2ma^{5/2}}\sin\frac{\pi x}{a} \ \ \ \ \ (20)$

The potential energy is

 $\displaystyle V\left|\psi\right\rangle$ $\displaystyle =$ $\displaystyle \frac{1}{2}m\omega^{2}x^{2}\psi\ \ \ \ \ (21)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{m\omega^{2}x^{2}}{2\sqrt{a}}\sin\frac{\pi x}{a} \ \ \ \ \ (22)$

Combining them we get

 $\displaystyle \left\langle H\right\rangle$ $\displaystyle =$ $\displaystyle \frac{1}{\sqrt{a}}\int_{-a/2}^{a/2}\left(\frac{\hbar^{2}\pi^{2}}{2ma^{5/2}}\sin^{2}\frac{\pi x}{a}+\frac{m\omega^{2}x^{2}}{2\sqrt{a}}\sin^{2}\frac{\pi x}{a}\right)dx\ \ \ \ \ (23)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{\left(2\pi^{2}-3\right)m^{2}\omega^{2}a^{4}+6\pi^{4}\hbar^{2}}{12ma^{2}} \ \ \ \ \ (24)$

To find the value of ${a}$ that minimizes ${\left\langle H\right\rangle }$ we take the derivative and set to zero as usual, and find that

 $\displaystyle a_{min}$ $\displaystyle =$ $\displaystyle \frac{6^{1/4}\pi}{\left(2\pi^{2}-3\right)^{1/4}}\sqrt{\frac{\hbar}{m\omega}}\ \ \ \ \ (25)$ $\displaystyle \left\langle H\right\rangle _{min}$ $\displaystyle =$ $\displaystyle \frac{\sqrt{6\left(2\pi^{2}-3\right)}}{6}\hbar\omega\ \ \ \ \ (26)$ $\displaystyle$ $\displaystyle \approx$ $\displaystyle 1.67\hbar\omega\ \ \ \ \ (27)$ $\displaystyle$ $\displaystyle >$ $\displaystyle \frac{3}{2}\hbar\omega \ \ \ \ \ (28)$

Again, we get a reasonable upper bound on the energy.