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 (cont.)

Posted by peeterjoot on October 7, 2011

Picking up where this problem was last abandoned.

[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.)]

Calculating the Fourier terms

Our wave function, with \beta=1 is plotted in figure (\ref{fig:expMinusBetsAbsX})

\caption{Exponential trial function with absolute exponential die off.}

The zeroth order fitting using the gaussian exponential is found to be

\begin{aligned}\psi_0(x) = \sqrt{2 \beta} \text{erfc}\left(\frac{\beta }{\sqrt{2} \alpha }\right)e^{- \alpha^2 x^2/2 +\beta^2/(2 \alpha^2)} \end{aligned} \hspace{\stretch{1}}(6.16)

With \alpha = \beta = 1, this is plotted in figure (\ref{fig:expMinusBetsAbsXfirstOrderFitting}) and can be seen to match fairly well

\caption{First ten orders, fitting harmonic oscillator wavefunctions to this trial function.}

The higher order terms get small fast, but we can see in figure (\ref{fig:expMinusBetsAbsXtenthOrderFitting}), where a tenth order fitting is depicted that it would take a number of them to get anything close to the sharp peak that we have in our exponential trial function.

\caption{Tenth order harmonic oscillator wavefunction fitting.}

Note that the all the brakets of even orders in n with the trial function are zero, which is why the tenth order approximation is only a sum of six terms.

Details for this \href{\

The question of interest is why if we can approximate the trial function so nicely (except at the origin) even with just a first order approximation (polynomial times gaussian functions where the polynomials are Hankel functions), and we can get an exact value for the lowest energy state using the first order approximation of our trial function, why do we get garbage if we include enough terms that the peak is sharp? It must therefore be important to consider the origin, but how do we give some meaning to the derivative of the absolute value function? The key, supplied when asking in office hours for the course, is to express the absolute value function in terms of Heavyside step functions, for which the derivative can be identified as the delta function.

Correcting, treating the origin this way.

Here’s how we can express the absolute value function using the Heavyside step

\begin{aligned}{\left\lvert{x}\right\rvert} = x \theta(x) - x \theta(-x),\end{aligned} \hspace{\stretch{1}}(7.17)

where the step function is zero for x  0 as plotted in figure (\ref{fig:stepFunction}).


Expressed this way, with the identification \theta'(x) = \delta(x), we have for the derivative of the absolute value function

\begin{aligned}{\left\lvert{x}\right\rvert}' &= x' \theta(x) - x' \theta(-x) + x \theta'(x) - x \theta'(-x) \\ &= \theta(x) - \theta(-x) + x \delta(x) + x \delta(-x) \\ &= \theta(x) - \theta(-x) + x \delta(x) + x \delta(x) \\ \end{aligned}

Observe that we have our expected unit derivative for x > 0, and -1 derivative for x < 0. At the origin our \theta contributions vanish, and we are left with

\begin{aligned}{\left.{{{\left\lvert{x}\right\rvert}' }}\right\vert}_{{x=0}}= 2 {\left.{{x \delta(x)}}\right\vert}_{{x=0}} \\ \end{aligned}

We’ve got zero times infinity here, so how do we give meaning to this? As with any delta functional, we’ve got to apply it to a well behaved (square integrable) test function f(x) and integrate. Doing so we have

\begin{aligned}\int_{-\infty}^\infty dx {\left\lvert{x}\right\rvert}' f(x) &= 2 \int_{-\infty}^\infty dx x \delta(x) f(x) \\ &= 2 (0) f(0)\end{aligned}

This equals zero for any well behaved test function f(x). Since the delta function only picks up the contribution at the origin, we can therefore identify {\left\lvert{x}\right\rvert}' as zero at the origin.

Using the same technique, we can express our trial function in terms of steps

\begin{aligned}\psi = e^{-\beta {\left\lvert{x}\right\rvert}} = \theta(x) e^{-\beta x} + \theta(-x) e^{\beta x}.\end{aligned} \hspace{\stretch{1}}(7.18)

This we can now take derivatives of, even at the origin, and find

\begin{aligned}\psi' &= \theta'(x) e^{-\beta x} + \theta'(-x) e^{\beta x} -\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \\ &= \delta(x) e^{-\beta x} - \delta(-x) e^{\beta x} -\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \\ &= \not{{\delta(x) e^{-\beta x} - \delta(x) e^{\beta x}}} -\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \\ &= \beta \left( -\theta(x) e^{-\beta x} + \theta(-x) e^{\beta x} \right)\end{aligned}

Taking second derivatives we find

\begin{aligned}\psi''&= \beta \left( -\theta'(x) e^{-\beta x} + \theta'(-x) e^{\beta x} +\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \right) \\ &= \beta \left( -\delta(x) e^{-\beta x} - \delta(-x) e^{\beta x} +\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \right) \\ &= \beta^2 \psi - 2 \beta \delta(x)\end{aligned}

Now application of the Hamiltonian operator on our trial function gives us

\begin{aligned}H \psi = -\frac{\hbar^2}{2m} \left( \beta^2 \psi - 2 \delta(x) \right) + \frac{1}{{2}} m \omega^2 x^2 \psi,\end{aligned} \hspace{\stretch{1}}(7.19)


\begin{aligned}{\langle {\psi} \rvert} H {\lvert {\psi} \rangle} &= \int_{-\infty}^\infty \left( -\frac{\hbar^2 \beta^2}{2m} + \frac{1}{{2}} m \omega^2 x^2 \right) e^{-2 \beta {\left\lvert{x}\right\rvert} }+ \frac{\hbar^2 \beta}{m}\int_{-\infty}^\infty \delta(x)e^{- \beta {\left\lvert{x}\right\rvert}} \\ &=-\frac{\beta \hbar^2}{2m} + \frac{m \omega^2}{4 \beta^3} + \frac{\hbar^2 \beta}{m} \\ &=\frac{\beta \hbar^2}{2m} + \frac{m \omega^2}{4 \beta^3}.\end{aligned}

Normalized we have

\begin{aligned}E[\beta] = \frac{{\langle {\psi} \rvert} H {\lvert {\psi} \rangle}}{\left\langle{{\psi}} \vert {{\psi}}\right\rangle} = \frac{\beta^2 \hbar^2}{2m} + \frac{m \omega^2}{4 \beta^2}.\end{aligned} \hspace{\stretch{1}}(7.20)

This is looking much more promising. We’ll have the sign alternation that we require to find a positive, non-complex, value for \beta when E[\beta] is minimized. That is

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

so the extremum is found at

\begin{aligned}\beta^4 = \frac{m^2 \omega^2}{2 \hbar^2}.\end{aligned} \hspace{\stretch{1}}(7.22)

Plugging this back in we find that our trial function associated with the minimum energy (unnormalized still) is

\begin{aligned}\psi = e^{-\sqrt{\frac{m \omega x^2}{\sqrt{2} \hbar}}},\end{aligned} \hspace{\stretch{1}}(7.23)

and that energy, after substitution, is

\begin{aligned}E[\beta_{\text{min}}] = \frac{\hbar \omega}{2} \sqrt{2}\end{aligned} \hspace{\stretch{1}}(7.24)

We have something that’s 1.4 \times the true ground state energy, but is at least a ball park value. However, to get this result, we have to be very careful to treat our point of singularity. A derivative that we’d call undefined in first year calculus, is not only defined, but required, for this treatment to work!


2 Responses to “Curious problem using the variational method to find the ground state energy of the Harmonic oscillator (cont.)”

  1. Hi Peeter,

    I was working through your algebra from your previous post and I noticed you picked up an extra minus sign somewhere, which led to an imaginary number after you took a square root. You arrive at the same energy as before, but this time around you didn’t the factor of i. I double checked my work on paper and through Mathematica, and I can confirm that the minimum energy is indeed what you find. Solving the problem through the standard methods as before — By splitting the integral into separate portions (or equivalently integrating over one half of the whole plane and doubling that result) — I was able to obtain the correct expression for beta squared which was \pm \frac{m \omega}{\sqrt{2} \hbar}.

    Specifically, within the pdf you have available on this post I see the error is between equations 9 and 10. You’ll see that when you differentiate 9 with respect to beta, one of the minus signs from equation 9 disappeared. If you fix that, you get the same exact results onwards as you did here.

    You can try the following Mathematica expression and you’ll see what I mean:

    Solve[D[(\[Beta]^2 \[HBar]^2)/(2 m) + (m \[Omega]^2)/(
    4 \[Beta]^2), \[Beta]] == 0, \[Beta]]

    • peeterjoot said

      The extra minus sign in (4.10) (or (9) in the pdf) is a typo. If you look back to just before ‘Note that evaluating this’ in the previous post (or just before ‘A naive evaluation of this integral’ in the pdf which is now a bit different), you’ll see there’s only one negative sign, and it appears I’ve corrected things on the next line. This still leaves the imaginary trouble, only resolved by the delta function handling of the absolute value’s derivative at the origin.

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: