A.25 Details of the animations

This note explains how the wave packet animations of chapter 7.11 and 7.12 were obtained. If you want a better understanding of unsteady solutions of the Schrödinger equation and their boundary conditions, this is a good place to start. In fact, deriving such solutions is a popular item in quantum mechanics books for physicists.

First consider the wave packet of the particle in free space, as shown in chapter 7.11.1. An energy eigenfunction with energy $E$ in free space takes the general form

$$\psi_E = C_{\rm {f}} e^{{\rm i}p x/\hbar} + C_{\rm {b}} e^{-{\rm i}p x/\hbar} \qquad p=\sqrt{2mE}$$

where $p$ is the momentum of the particle and $C_{\rm {f}}$ and $C_{\rm {b}}$ are constants.

To study a single wave packet coming in from the far left, the coefficient $C_{\rm {b}}$ has to be set to zero. The reason was worked out in chapter 7.10: combinations of exponentials of the form $C_{\rm {b}}e^{-{\rm i}{p}x/\hbar}$ produce wave packets that propagate backwards in $x$, from right to left. Therefore, a nonzero value for $C_{\rm {b}}$ would add an unwanted second wave packet coming in from the far right.

Figure A.9: Example energy eigenfunction for the particle in free space.

With only the coefficient $C_{\rm {f}}$ of the forward moving part left, you may as well scale the eigenfunction so that $C_{\rm {f}}=1$, simplifying it to

$$\psi_E=e^{{\rm i}p x/\hbar}$$

A typical example is shown in figure A.9. Plus and minus the magnitude of the eigenfunction are shown in black, and the real part is shown in red. This wave function is an eigenfunction of linear momentum, with $p$ the linear momentum.

To produce a coherent wave packet, eigenfunctions with somewhat different energies $E$ have to be combined together. Since the momentum is given by $p=\sqrt{2mE}$, different energy means different momentum $p$; therefore the wave packet can be written as

$$\Psi(x,t) = \int_{{\rm all\ }p} c(p) e^{- {\rm i}E t/\hbar} \psi_E(x) {\,\rm d}p %$$ (A.205)

where $c(p)$ is some function that is only nonzero in a relatively narrow range of momenta $p$ around the nominal momentum. Except for that basic requirement, the choice of the function $c(p)$ is quite arbitrary. Choose some suitable function $c(p)$, then use a computer to numerically integrate the above integral at a large number of plot points and times. Dump the results into your favorite animation software and bingo, out comes the movie.

Figure A.10: Example energy eigenfunction for a particle entering a constant accelerating force field.

Next consider the animation of chapter 7.11.2, where the particle accelerates along a downward potential energy ramp starting from point A. A typical energy eigenfunction is shown in figure A.10. Since to the left of point A, the potential energy is still zero, in that region the energy eigenfunction is still of the form

$$\psi_E = C^{\rm {l}}_{\rm {f}} e^{{\rm i}p_{\rm {c}}^{\rm ... ...{ for }x< x_{\rm {A}} \qquad p_{\rm {c}}^{\rm {l}}=\sqrt{2mE}$$

where $p_{\rm {c}}^{\rm {l}}$ is the momentum that a classical particle of energy $E$ would have in the left region. (Quantum mechanics looks at the complete wave function, not just a single point of it, and would say that the momentum is uncertain.)

In this case, it can no longer be argued that the coefficient $C^{\rm {l}}_{\rm {b}}$ must be zero to avoid a packet entering from the far right. After all, the $C^{\rm {l}}_{\rm {b}}e^{-{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$ term does not extend to the far right anymore. To the right of point $A$, the potential changes linearly with position, and the exponentials are no longer valid.

In fact, it is known that the solution of the Hamiltonian eigenvalue problem in a region with a linearly varying potential is a combination of two weird functions Ai and Bi that are called the “Airy” functions. The bad news is that if you are interested in learning more about their properties, you will need an advanced mathematical handbook like [1] or at least look at addendum {A.27}. The good news is that free software to evaluate these functions and their first derivatives is readily available on the web. The general solution for a linearly varying potential is of the form

$$\psi_E = C_{\rm {B}} {\rm Bi}(\overline{x}) + C_{\rm {A}... ...2}}\frac{V-E}{V'} \quad V' \equiv \frac{{\rm d}V}{{\rm d}x}$$

Note that $(V-E)/V'$ is the $x$-position measured from the point where $V=E$. Also note that the cube root is negative, so that $\overline{x}$ is.

It may be deduced from the approximate analysis of addendum {A.26} that to prevent a second wave packet coming in from the far right, Ai and Bi must appear together in the combination ${\rm {Bi}}+{\rm i}{\rm {Ai}}$ as shown in figure A.10. The fact that no second packet comes in from the far right in the animation can be taken as an experimental confirmation of that result, so there seems little justification to go over the messy argument.

To complete the determination of the eigenfunction for a given value of $E$, the constants $C^{\rm {l}}_{\rm {f}}$, $C^{\rm {l}}_{\rm {b}}$ and $C^{\rm {r}}$ must still be determined. That goes as follows. For now, assume that $C^{\rm {r}}$ has the provisional value $c^{\rm {r}}=1$. Then provisional values $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$ for the other two constants may be found from the requirements that the left and right regions give the same values for $\psi_E$ and ${\rm d}\psi_E/{\rm d}{x}$ at the point A in figure A.10 where they meet:

$$\begin{eqnarray*} & \displaystyle c^{\rm {l}}_{\rm {f}}e^{{\rm i}p_{\rm {c}}... ...ine{x}_{\rm {A}})\right] \frac{{\rm d}\overline{x}}{{\rm d}x} \end{eqnarray*}$$

That is equivalent to two equations for the two constants $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$, since everything else can be evaluated, using the mentioned software. So $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$ can be found from solving these two equations.

As the final step, it is desirable to normalize the eigenfunction $\psi_E$ so that $C^{\rm {l}}_{\rm {f}}=1$. To do so, the entire provisional eigenfunction can be divided by $c^{\rm {l}}_{\rm {f}}$, giving $C^{\rm {l}}_{\rm {b}}=c^{\rm {l}}_{\rm {b}}/c^{\rm {l}}_{\rm {f}}$ and $C^{\rm {r}}=c^{\rm {r}}/c^{\rm {l}}_{\rm {f}}$. The energy eigenfunction has now been found. And since $C^{\rm {l}}_{\rm {f}}=1$, the $e^{{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$ term is exactly the same as the free space energy eigenfunction of the first example. That means that if the eigenfunctions $\psi_E$ are combined into a wave packet in the same way as in the free space case, (A.205) with $p$ replaced by $p_{\rm {c}}^{\rm {l}}$, the $e^{{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$ terms produce the exact same wave packet coming in from the far left as in the free space case.

For larger times, the $C^{\rm {l}}_{\rm {b}}e^{-{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$ terms produce a “reflected” wave packet that returns toward the far left. Note that $e^{-{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$ is the complex conjugate of $e^{{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$, and it can be seen from the unsteady Schrödinger equation that if the complex conjugate of a wave function is taken, it produces a reversal of time. Wave packets coming in from the far left at large negative times become wave packets leaving toward the far left at large positive times. However, the constant $C^{\rm {l}}_{\rm {b}}$ turns out to be very small in this case, so there is little reflection.

Figure A.11: Example energy eigenfunction for a particle entering a constant decelerating force field.

Next consider the animation of chapter 7.11.3, where the particle is turned back by an upward potential energy ramp. A typical energy eigenfunction for this case is shown in figure A.11. Unlike in the previous example, where the argument $\overline{x}$ of the Airy functions was negative at the far right, here it is positive. Table books that cover the Airy functions will tell you that the Airy function Bi blows up very strongly with increasing positive argument $\overline{x}$. Therefore, if the solution in the right hand region would involve any amount of Bi, it would locate the particle at infinite $x$ for all times. For a particle not at infinity, the solution in the right hand region can only involve the Airy function Ai. That function decays rapidly with positive argument $\overline{x}$, as seen in figure A.11.

The further determination of the energy eigenfunctions proceeds along the same lines as in the previous example: give $C^{\rm {r}}$ a provisional value $c^{\rm {r}}=1$, then compute $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$ from the requirements that the left and right regions produce the same values for $\psi$ and ${\rm d}\psi/{\rm d}{x}$ at the point A where they meet. Finally divide the eigenfunction by $c^{\rm {l}}_{\rm {f}}$. The big difference is that now $C^{\rm {l}}_{\rm {b}}$ is no longer small; $C^{\rm {l}}_{\rm {b}}$ turns out to be of unit magnitude just like $C^{\rm {l}}_{\rm {f}}$. It means that the incoming wave packet is reflected back completely.

Figure A.12: Example energy eigenfunction for the harmonic oscillator.

For the harmonic oscillator of chapter 7.11.4, the analysis is somewhat different. In particular, chapter 4.1.2 showed that the energy levels of the one-dimensional harmonic oscillator are discrete,

$$E_n = \frac{2n+1}2 \hbar \omega \mbox{ for $n = 0,1,2,\ldots$}$$

so that unlike the motions just discussed, the solution of the Schrödinger equation is a sum, rather than the integral (A.205),

$$\Psi(x,t) = \sum_{n=0}^\infty c_n e^{-{\rm i}E_n t/\hbar} h_n(x)$$

However, for large $n$ the difference between summation and integration is small.

Also, while the energy eigenfunctions $h_n(x)$ are not exponentials as for the free particle, for large $n$ they can be pairwise combined to approximate such exponentials. For example, eigenfunction $h_{50}$, shown in figure A.12, behaves near the center point much like a cosine if you scale it properly. Similarly, $h_{51}$ behaves much like a sine. A cosine plus ${\rm i}$ times a sine gives an exponential, according to the Euler formula (2.5). Create similar exponential combinations of eigenfunctions with even and odd values of $n$ for a range of $n$ values, and there are the approximate exponentials that allow you to create a wave packet that is at the center point at time $t=0$. In the animation, the range of $n$ values was centered around $n=50$, making the nominal energy hundred times the ground state energy. The exponentials degenerate over time, since their component eigenfunctions have slightly different energy, hence time evolution. That explains why after some time, the wave packet can return to the center point going the other way.

Figure A.13: Example energy eigenfunction for a particle encountering a brief accelerating force.

For the particle of chapter 7.12.1 that encounters a brief accelerating force, an example eigenfunction looks like figure A.13. In this case, the solution in the far right region is similar to the one in the far left region. However, there cannot be a term of the form $e^{-{\rm i}p_{\rm {c}}^{\rm {r}} x/\hbar}$ in the far right region, because when the eigenfunctions are combined, it would produce an unwanted wave packet coming in from the far right. In the middle region of linearly varying potential, the wave function is again a combination of the two Airy functions. The way to find the constants now has an additional step. First give the constant $C^{\rm {r}}$ of the far right exponential the provisional value $c^{\rm {r}}=1$ and from that, compute provisional values $c^{\rm {m}}_{\rm {A}}$ and $c^{\rm {m}}_{\rm {B}}$ by demanding that the Airy functions give the same values for $\psi$ and ${\rm d}\psi/{\rm d}{x}$ as the far-right exponential at point B, where they meet. Next compute provisional values $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$ by demanding that the far-left exponentials give the same values for $\psi$ and ${\rm d}\psi/{\rm d}{x}$ as the Airy functions at point A, where they meet. Finally, divide all the constants by $c^{\rm {l}}_{\rm {f}}$ to make $C^{\rm {l}}_{\rm {f}}=1$.

Figure A.14: Example energy eigenfunction for a particle tunneling through a barrier.

For the tunneling particle of chapter 7.12.2, an example eigenfunction is as shown in figure A.14. In this case, the solution in the middle part is not a combination of Airy functions, but of real exponentials. It is essentially the same solution as in the left and right parts, but in the middle region the potential energy is greater than the total energy, making $p_{\rm {c}}^{\rm {m}}=\sqrt{2m(E-V_{\rm {m}})}$ an imaginary number. Therefore the arguments of the exponentials become real when written in terms of the absolute value of the momentum $\vert p_{\rm {c}}^{\rm {m}}\vert=\sqrt{2m(V_{\rm {m}}-E)}$. The rest of the analysis is similar to that of the previous example.

Figure A.15: Example energy eigenfunction for tunneling through a delta function barrier.

For the particle tunneling through the delta function potential in chapter 7.12.2, an example energy eigenfunction is shown in figure A.15. The potential energy in this case is $V=\nu\delta(x-x_{\rm {A}})$, where $\delta(x-x_{\rm {A}})$ is a spike at point A that integrates to one and the strength $\nu$ is a chosen constant. In the example, $\nu$ was chosen to be $\sqrt{2\hbar^2E_{\rm {nom}}/m}$ with $E_{\rm {nom}}$ the nominal energy. For that strength, half the wave packet will pass through.

For a delta function potential, a modification must be made in the analysis as used so far. As figure A.15 illustrates, there are kinks in the energy eigenfunction at the location A of the delta function. The left and right expressions for the eigenfunction do not predict the same value for its derivative ${\rm d}\psi/{\rm d}{x}$ at point A. To find the difference, integrate the Hamiltonian eigenvalue problem from a point a very short distance $\varepsilon$ before point A to a point the same very short distance behind it:

$$- \frac{\hbar^2}{2m} \int_{x=x_{\rm {A}}-\varepsilon}^{x... ...{A}}-\varepsilon}^{x_{\rm {A}}+\varepsilon} E \psi {\,\rm d}x$$

The integral in the right hand side is zero because of the vanishingly small interval of integration. But the delta function spike in the left hand side integrates to one regardless of the small integration range, so

$$- \frac{\hbar^2}{2m} \frac{{\rm d}\psi}{{\rm d}x}\bigg\v... ...ilon}^{x_{\rm {A}}+\varepsilon} + \nu \psi(x_{\rm {A}}) = 0$$

For vanishingly small $\varepsilon$, ${\rm d}\psi/{\rm d}{x}$ at $x_{\rm {A}}+\varepsilon$ becomes what the right hand part of the eigenfunction gives for ${\rm d}\psi/{\rm d}{x}$ at $x_{\rm {A}}$, while ${\rm d}\psi/{\rm d}{x}$ at $x_{\rm {A}}-\varepsilon$ becomes what the left hand part gives for it. As seen from the above equation, the difference is not zero, but $2m\nu\psi(x_{\rm {A}})/\hbar^2$.

So the correct equations for the provisional constants are in this case

$$\begin{eqnarray*} & \displaystyle c^{\rm {l}}_{\rm {f}}e^{{\rm i}p_{\rm {c}}... ... c^{\rm {r}} e^{{\rm i}p_{\rm {c}}^{\rm {r}} x_{\rm {A}}/\hbar} \end{eqnarray*}$$

Compared to the analysis as used previously, the difference is the final term in the second equation that is added by the delta function.

The remainder of this note gives some technical details for if you are actually planning to do your own animations. It is a good idea to assume that the units of mass, length, and time are chosen such that $\hbar$ and the nominal energy are one, while the mass of the particle is one-half. That avoids having to guesstimate suitable values for all sorts of very small numbers. The Hamiltonian eigenvalue problem then simplifies to

$$- \frac{{\rm d}^2\psi}{{\rm d}x^2} + V \psi = E \psi$$

where the values of $E$ of interest cluster around 1. The nominal momentum will be one too. In those units, the length of the plotted range was one hundred in all but the harmonic oscillator case.

It should be noted that to select a good function $c(p)$ in (A.205) is somewhat of an art. The simplest idea would be to choose $c(p)$ equal to one in some limited range around the nominal momentum, and zero elsewhere, as in

$$c(p) = 1\quad \mbox{if } (1-r) p_{\rm nom} < p < (1+r)p_{\rm nom} \qquad c(p) = 0\quad \mbox{otherwise}$$

where $r$ is the relative deviation from the nominal momentum below which $c(p)$ is nonzero. However, it is know from Fourier analysis that the locations where $c(p)$ jumps from one to zero lead to lengthy wave packets when viewed in physical space. {D.45}. Functions $c(p)$ that do lead to nice compact wave packets are known to be of the form

$$c(p) = \exp\left(-\frac{(p-p_{\rm nom})^2}{r^2 p_{\rm nom}^2}\right)$$

And that is essentially the function $c(p)$ used in this study. The typical width of the momentum range was chosen to be $r=0.15$, or 15%, by trial and error. However, it is nice if $c(p)$ becomes not just very small, but exactly zero beyond some point, for one because it cuts down on the number of energy eigenfunctions that have to be evaluated numerically. Also, it is nice not to have to worry about the possibility of $p$ being negative in writing energy eigenfunctions. Therefore, the final function used was

$$c(p) = \exp\left( -\frac{(p-p_{\rm nom})^2}{r^2[p_{\rm... ...or } 0 < p < 2p_{\rm nom} \quad c(p) = 0\mbox{ otherwise}$$

The actual difference in numerical values is small, but it does make $c(p)$ exactly zero for negative momenta and those greater than twice the nominal value. Strictly speaking, $c(p)$ should still be multiplied by a constant to make the total probability of finding the particle equal to one. But if you do not tell people what numbers for $\Psi$ are on the vertical axes, you do not need to bother.

In doing the numerical integrations to find $\Psi(x,t)$, note that the mid point and trapezium rules of numerical integration are exponentially accurate under the given conditions, so there is probably not much motivation to try more advanced methods. The mid point rule was used.

The animations in this book used the numerical implementations daie.f, dbie.f, daide.f, and dbide.f from netlib.org for the Airy functions and their first derivatives. These offer some basic protection against underflow and overflow by splitting off an exponential for positive $\overline{x}$. It may be a good idea to check for underflow and overflow in general and use 64 bit precision. The examples here did.

For the harmonic oscillator, the larger the nominal energy is compared to the ground state energy, the more the wave packet can resemble a single point compared to the limits of motion. However, the computer program used to create the animation computed the eigenfunctions by evaluating the analytical expression given in derivation {D.12}, and explicitly evaluating the Hermite polynomials is very round-off sensitive. That limited it to a maximum of about hundred times the ground state energy when allowing for enough uncertainty to localize the wave packet. Round-off is a general problem for power series, not just for the Hermite polynomials. If you want to go to higher energies to get a smaller wave packet, you will want to use a finite difference or finite element method to find the eigenfunctions.

The plotting software used to produce the animations was a mixture of different programs. There are no doubt much simpler and better ways of doing it. In the animations presented here, first plots were created of $\Psi$ versus $x$ for a large number of closely spaced times covering the duration of the animation. These plots were converted to gifs using a mixture of personal software, netpbm, and ghostview. The gifs were then combined into a single movie using gifsicle.