Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Curious problem using the variational method to find the ground state energy of the Harmonic oscillator.

Posted by peeterjoot on October 3, 2011

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]


Recap. Variational method to find the ground state energy.

Problem 3 of section 24.4 in the text [1] is an interesting one. It asks to use the variational method to find the ground state energy of a one dimensional harmonic oscillator Hamiltonian.

Somewhat unexpectedly, once I take derivatives equate to zero, I find that the variational parameter beta becomes imaginary?

I tried this twice on paper and pencil, both times getting the same thing. This seems like a noteworthy problem, and one worth reflecting on a bit.

Recap. The variational method.

Given any, not necessarily normalized wavefunction, with a series representation specified using the energy eigenvectors for the space

\begin{aligned}{\lvert {\psi} \rangle} = \sum_m c_{m} {\lvert {\psi_m} \rangle},\end{aligned} \hspace{\stretch{1}}(3.1)


\begin{aligned}H {\lvert {\psi_m} \rangle} = E_m {\lvert {\psi_m} \rangle},\end{aligned} \hspace{\stretch{1}}(3.2)


\begin{aligned}\left\langle{{\psi_m}} \vert {{\psi_n}}\right\rangle = \delta_{mn}.\end{aligned} \hspace{\stretch{1}}(3.3)

We can perform an energy expectation calculation with respect to this more general state

\begin{aligned}{\langle {\psi} \rvert} H {\lvert {\psi} \rangle} &= \sum_m c_{m}^{*} {\langle {\psi_m} \rvert} H\sum_n c_{n} {\lvert {\psi_n} \rangle} \\ &=\sum_m c_{m}^{*} {\langle {\psi_m} \rvert}\sum_n c_{n} E_n {\lvert {\psi_n} \rangle} \\ &=\sum_{m,n} c_{m}^{*} c_n E_n \left\langle{{\psi_m}} \vert {{\psi_n}}\right\rangle \\ &\sum_{m} {\left\lvert{c_{m}}\right\rvert}^2 E_m \\ &\ge\sum_{m} {\left\lvert{c_{m}}\right\rvert}^2 E_0&=E_0 \left\langle{{\psi}} \vert {{\psi}}\right\rangle\end{aligned}

This allows us to form an estimate of the ground state energy for the system, by using any state vector formed from a superposition of energy eigenstates, by simply calculating

\begin{aligned}E_0 \le \frac{{\langle {\psi} \rvert} H {\lvert {\psi} \rangle}}{ \left\langle{{\psi}} \vert {{\psi}}\right\rangle }.\end{aligned} \hspace{\stretch{1}}(3.4)

One of the examples in the text is to use this to find an approximation of the ground state energy for the Helium atom Hamiltonian

\begin{aligned}H = -\frac{\hbar^2}{2m} \left( \boldsymbol{\nabla}_1^2+\boldsymbol{\nabla}_1^2\right) - 2 e^2 \left( \frac{1}{{r_1}} + \frac{1}{{r_2}} \right) + \frac{e^2}{{\left\lvert{\mathbf{r}_1 - \mathbf{r}_2}\right\rvert}}.\end{aligned} \hspace{\stretch{1}}(3.5)

This calculation is performed using a trial function that was a solution of the interaction free Hamiltonian

\begin{aligned}\phi = \frac{Z^3}{pi a_0^3} e^{-Z (r_1 + r_2)/a_0 }.\end{aligned} \hspace{\stretch{1}}(3.6)

This is despite the fact that this is not a solution to the interaction Hamiltonian. The end result ends up being pretty close to the measured value (although there is a pesky error in the book that appears to require a compensating error somewhere else).

Part of the variational technique used in that problem, is to allow Z to vary, and then once the normalized expectation is computed, set the derivative of that equal to zero to calculate the trial wavefunction as a parameter of Z that has the lowest energy eigenstate for a function of that form. We find considering the Harmonic oscillator that this final variation does not necessarily produce meaningful results.

The Harmonic oscillator variational problem.

The problem asks for the use of the trial wavefunction

\begin{aligned}\phi = e^{-\beta {\left\lvert{x}\right\rvert}},\end{aligned} \hspace{\stretch{1}}(4.7)

to perform the variational calculation above for the Harmonic oscillator Hamiltonian, which has the one dimensional position space representation

\begin{aligned}H = -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + \frac{1}{{2}} m \omega^2 x^2.\end{aligned} \hspace{\stretch{1}}(4.8)

We can find the normalization easily

\begin{aligned}\left\langle{{\phi}} \vert {{\phi}}\right\rangle &= \int_{-\infty}^\infty e^{- 2 \beta {\left\lvert{x}\right\rvert}} dx \\ &= 2 \frac{1}{{2 \beta}} \int_{0}^\infty e^{- 2 \beta x} 2 \beta dx \\ &= 2 \frac{1}{{2 \beta}} \int_{0}^\infty e^{- u} du \\ &= \frac{1}{{\beta}}\end{aligned}

After a bit of calculation, using integration by parts, we find for the energy expectation

\begin{aligned}{\langle {\phi} \rvert} H {\lvert {\phi} \rangle} = \int_{-\infty}^\infty dxe^{- \beta {\left\lvert{x}\right\rvert}} \left( -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + \frac{1}{{2}} m \omega^2 x^2 \right)e^{- \beta {\left\lvert{x}\right\rvert}} = -\frac{\beta \hbar^2}{2m} + \frac{m \omega^2}{4 \beta^3}.\end{aligned} \hspace{\stretch{1}}(4.9)

Note that evaluating this integral requires the origin to be avoided where the derivative of {\left\lvert{x}\right\rvert} becomes undefined. So one has to evaluate this in the limit as \int_{-\infty}^\infty = \int_{\infty}^{-\epsilon} + \int_\epsilon^\infty (which is an easy way to do it anyways since the absolute value can be eliminated by doubling the integral). Because of the curious end result, I also verified my calculation using Mathematica.

So, our ground state energy estimation, parametrized by \beta becomes

\begin{aligned}E[\beta] = -\frac{\beta^2 \hbar^2}{2m} + \frac{m \omega^2}{4 \beta^2}.\end{aligned} \hspace{\stretch{1}}(4.10)


Observe that if we set the derivative of this equal to zero to find the “best” beta associated with this trial function

\begin{aligned}0 = \frac{\partial {E}}{\partial {\beta}} = -\frac{\beta \hbar^2}{2m} - \frac{m \omega^2}{2 \beta^3}\end{aligned} \hspace{\stretch{1}}(4.11)

we find that the parameter beta that best minimizes this ground state energy function is complex with value

\begin{aligned}\beta^2 = \pm \frac{i m \omega}{\sqrt{2} \hbar}.\end{aligned} \hspace{\stretch{1}}(4.12)

So, it appears that we can’t minimize 4.10 to find a best ground state energy estimate associated with the trial function 4.7. We do however, know the exact ground state energy \hbar \omega/2 for the Harmonic oscillator. Is is possible to show that for all \beta^2 we have

\begin{aligned}\frac{\hbar \omega}{2} \le -\frac{ \beta^2 \hbar^2}{2m} + \frac{m \omega^2}{4 \beta^2}\end{aligned} \hspace{\stretch{1}}(4.13)

? This inequality would be expected if we can assume that the trial wavefunction has a Fourier series representation utilizing the actual energy eigenfunctions for the system.

Some reflection.

I think that it is notable that I don’t believe the trial wave function for this problem lies in the span of the Hilbert space that describes the solutions to the Harmonic oscillator. Another thing of possible interest is the trouble near the origin for this wave function, when operated on by P^2/2m. Observe that Mathematica also required some special hand holding to deal with the origin.

I had initially thought that part of the value of this variational method was that we can use it despite not even knowing what the exact solution is (and in the case of the Helium atom, I believe it was stated in class that an exact closed form solution is not even known). This makes me wonder what restrictions must be imposed on the trial solutions to get a meaningful answer from the variational calculation?

Having started with a wavefunction that is probably not representable in the solution space is likely the bigger problem here. We probably need to adjust the treatment to account for that. Suppose we have

\begin{aligned}{\lvert {\phi} \rangle} = \sum_n c_n {\lvert {\psi_n} \rangle} + c_\perp {\lvert {\psi_\perp} \rangle}.\end{aligned} \hspace{\stretch{1}}(5.14)

where {\lvert {\psi_\perp} \rangle} is unknown, and presumed not orthogonal to any of the energy eigenkets. We can still calculate the norm of the trial function

\begin{aligned}\left\langle{{\phi}} \vert {{\phi}}\right\rangle&=\sum_{n,m} \left\langle{{ c_n \psi_n + c_\perp \psi_\perp}} \vert {{ c_m \psi_m + c_\perp \psi_\perp}}\right\rangle \\ &=\sum_n {\left\lvert{c_n}\right\rvert}^2 + c_n^{*} c_\perp \left\langle{{\psi_n}} \vert {{\psi_\perp}}\right\rangle+ c_n c_\perp^{*} \left\langle{{\psi_\perp}} \vert {{\psi_n}}\right\rangle+ {\left\lvert{c_\perp}\right\rvert}^2\left\langle{{\psi_\perp}} \vert {{\psi_\perp}}\right\rangle \\ &=\left\langle{{\psi_\perp}} \vert {{\psi_\perp}}\right\rangle +\sum_n {\left\lvert{c_n}\right\rvert}^2 + 2 \text{Real} \left(c_n^{*} c_\perp \left\langle{{\psi_n}} \vert {{\psi_\perp}}\right\rangle \right).\end{aligned}

Similarly we can calculate the energy expectation for this unnormalized state and find

\begin{aligned}{\langle {\phi} \rvert} H {\lvert {\phi} \rangle}&=\sum_{n,m} {\langle { c_n \psi_n + c_\perp \psi_\perp} \rvert} H {\lvert { c_m \psi_m + c_\perp \psi_\perp} \rangle} \\ &=\sum_n {\left\lvert{c_n}\right\rvert}^2 E_n+ c_n^{*} c_\perp E_n\left\langle{{\psi_n}} \vert {{\psi_\perp}}\right\rangle+ c_n c_\perp^{*} E_n \left\langle{{\psi_\perp}} \vert {{\psi_n}}\right\rangle+ {\left\lvert{c_\perp}\right\rvert}^2{\langle {\psi_\perp} \rvert} H {\lvert {\psi_\perp} \rangle} \end{aligned}

Our normalized energy expectation is therefore the considerably messier

\begin{aligned}\begin{aligned}\frac{{\langle {\phi} \rvert} H {\lvert {\phi} \rangle}}{\left\langle{{\phi}} \vert {{\phi}}\right\rangle}&=\frac{\sum_n {\left\lvert{c_n}\right\rvert}^2 E_n+ c_n^{*} c_\perp E_n\left\langle{{\psi_n}} \vert {{\psi_\perp}}\right\rangle+ c_n c_\perp^{*} E_n \left\langle{{\psi_\perp}} \vert {{\psi_n}}\right\rangle+ {\left\lvert{c_\perp}\right\rvert}^2{\langle {\psi_\perp} \rvert} H {\lvert {\psi_\perp} \rangle} }{\left\langle{{\psi_\perp}} \vert {{\psi_\perp}}\right\rangle +\sum_m {\left\lvert{c_m}\right\rvert}^2 + 2 \text{Real} \left(c_m^{*} c_\perp \left\langle{{\psi_m}} \vert {{\psi_\perp}}\right\rangle \right)} \\ &\ge \frac{\sum_n {\left\lvert{c_n}\right\rvert}^2 E_0+ c_n^{*} c_\perp E_n\left\langle{{\psi_n}} \vert {{\psi_\perp}}\right\rangle+ c_n c_\perp^{*} E_n \left\langle{{\psi_\perp}} \vert {{\psi_n}}\right\rangle+ {\left\lvert{c_\perp}\right\rvert}^2{\langle {\psi_\perp} \rvert} H {\lvert {\psi_\perp} \rangle} }{\left\langle{{\psi_\perp}} \vert {{\psi_\perp}}\right\rangle +\sum_m {\left\lvert{c_m}\right\rvert}^2 + 2 \text{Real} \left(c_m^{*} c_\perp \left\langle{{\psi_m}} \vert {{\psi_\perp}}\right\rangle \right)}\end{aligned}\end{aligned} \hspace{\stretch{1}}(5.15)

With a requirement to include the perpendicular cross terms the norm doesn’t just cancel out, leaving us with a clean estimation of the ground state energy. In order to utilize this variational method, we implicitly have an assumption that the \left\langle{{\psi_\perp}} \vert {{\psi_\perp}}\right\rangle and \left\langle{{\psi_m}} \vert {{\psi_\perp}}\right\rangle terms in the denominator are sufficiently small that they can be neglected.

For this harmonic oscillator problem, it would be interesting to calculate that Fourier remainder explicitly.


[1] BR Desai. Quantum mechanics with basic field theory. Cambridge University Press, 2009.

8 Responses to “Curious problem using the variational method to find the ground state energy of the Harmonic oscillator.”

  1. Peeter

    I ran through your calculations and I just wanted to point out a useful for trick integrating quantities such as absolute values within Mathematica.

    As you know, an expression of the form Integrate[Abs[f[x]], {x, -Infinity, +Infinity}] usually needs to be split up into Integrate[f[-x], {x, -Infinity, 0} + Integrate[f[x], {x, 0, Infinity}].

    To get around this, you can simply use Sqrt[f[x]^2] for the absolute value. This has all the properties of the absolute value and Mathematica can handle this much more cleanly than Abs[f[x]]. All the analytic derivatives of this slightly modified square root are well defined.

  2. […] Curious problem using the variational method to find the ground state energy of the Harmonic os… […]

  3. Ken said

    Hey Peeter

    I think your solution is in-correct. The double derivative of the trial wave function has a delta function. So, when you work out the expectation value of the kinetic part of the hamiltonian, you get a positive answer, which is turn makes your beta positive.
    So, the derivative of the wave function gives you a step function and the derivative of the step function is a delta function…

  4. Tom said

    I know this is way after the fact but following your link to section was very helpful after I was trying to solve Sakurai 5.20. Heavisides…Who would have thought! Thank-you!

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: