A.27 De­tails of the an­i­ma­tions

This note ex­plains how the wave packet an­i­ma­tions of chap­ter 7.11 and 7.12 were ob­tained. If you want a bet­ter un­der­stand­ing of un­steady so­lu­tions of the Schrö­din­ger equa­tion and their bound­ary con­di­tions, this is a good place to start. In fact, de­riv­ing such so­lu­tions is a pop­u­lar item in quan­tum me­chan­ics books for physi­cists.

First con­sider the wave packet of the par­ti­cle in free space, as shown in chap­ter 7.11.1. An en­ergy eigen­func­tion with en­ergy $E$ in free space takes the gen­eral 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 mo­men­tum of the par­ti­cle and $C_{\rm {f}}$ and $C_{\rm {b}}$ are con­stants.

To study a sin­gle wave packet com­ing in from the far left, the co­ef­fi­cient $C_{\rm {b}}$ has to be set to zero. The rea­son was worked out in chap­ter 7.10: com­bi­na­tions of ex­po­nen­tials of the form $C_{\rm {b}}e^{-{\rm i}{p}x/\hbar}$ pro­duce wave pack­ets that prop­a­gate back­wards in $x$, from right to left. There­fore, a nonzero value for $C_{\rm {b}}$ would add an un­wanted sec­ond wave packet com­ing in from the far right.

Fig­ure A.9: Ex­am­ple en­ergy eigen­func­tion for the par­ti­cle in free space.

With only the co­ef­fi­cient $C_{\rm {f}}$ of the for­ward mov­ing part left, you may as well scale the eigen­func­tion so that $C_{\rm {f}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1, sim­pli­fy­ing it to

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

A typ­i­cal ex­am­ple is shown in fig­ure A.9. Plus and mi­nus the mag­ni­tude of the eigen­func­tion are shown in black, and the real part is shown in red. This wave func­tion is an eigen­func­tion of lin­ear mo­men­tum, with $p$ the lin­ear mo­men­tum.

To pro­duce a co­her­ent wave packet, eigen­func­tions with some­what dif­fer­ent en­er­gies $E$ have to be com­bined to­gether. Since the mo­men­tum is given by $p$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\sqrt{2mE}$, dif­fer­ent en­ergy means dif­fer­ent mo­men­tum $p$; there­fore the wave packet can be writ­ten as

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

where $c(p)$ is some func­tion that is only nonzero in a rel­a­tively nar­row range of mo­menta $p$ around the nom­i­nal mo­men­tum. Ex­cept for that ba­sic re­quire­ment, the choice of the func­tion $c(p)$ is quite ar­bi­trary. Choose some suit­able func­tion $c(p)$, then use a com­puter to nu­mer­i­cally in­te­grate the above in­te­gral at a large num­ber of plot points and times. Dump the re­sults into your fa­vorite an­i­ma­tion soft­ware and bingo, out comes the movie.

Fig­ure A.10: Ex­am­ple en­ergy eigen­func­tion for a par­ti­cle en­ter­ing a con­stant ac­cel­er­at­ing force field.

Next con­sider the an­i­ma­tion of chap­ter 7.11.2, where the par­ti­cle ac­cel­er­ates along a down­ward po­ten­tial en­ergy ramp start­ing from point A. A typ­i­cal en­ergy eigen­func­tion is shown in fig­ure A.10. Since to the left of point A, the po­ten­tial en­ergy is still zero, in that re­gion the en­ergy eigen­func­tion is still of the form

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

where $p_{\rm {c}}^{\rm {l}}$ is the mo­men­tum that a clas­si­cal par­ti­cle of en­ergy $E$ would have in the left re­gion. (Quan­tum me­chan­ics looks at the com­plete wave func­tion, not just a sin­gle point of it, and would say that the mo­men­tum is un­cer­tain.)

In this case, it can no longer be ar­gued that the co­ef­fi­cient $C^{\rm {l}}_{\rm {b}}$ must be zero to avoid a packet en­ter­ing from the far right. Af­ter all, the $C^{\rm {l}}_{\rm {b}}e^{-{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$ term does not ex­tend to the far right any­more. To the right of point $A$, the po­ten­tial changes lin­early with po­si­tion, and the ex­po­nen­tials are no longer valid.

In fact, it is known that the so­lu­tion of the Hamil­ton­ian eigen­value prob­lem in a re­gion with a lin­early vary­ing po­ten­tial is a com­bi­na­tion of two weird func­tions Ai and Bi that are called the Airy func­tions. The bad news is that if you are in­ter­ested in learn­ing more about their prop­er­ties, you will need an ad­vanced math­e­mat­i­cal hand­book like [1] or at least look at ad­den­dum {A.29}. The good news is that free soft­ware to eval­u­ate these func­tions and their first de­riv­a­tives is read­ily avail­able on the web. The gen­eral so­lu­tion for a lin­early vary­ing po­ten­tial is of the form

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

Note that $(V-E)$$\raisebox{.5pt}{$/$}$​$V'$ is the $x$-​po­si­tion mea­sured from the point where $V$ $\vphantom0\raisebox{1.5pt}{$=$}$ $E$. Also note that the cube root is neg­a­tive, so that $\overline{x}$ is.

It may be de­duced from the ap­prox­i­mate analy­sis of ad­den­dum {A.28} that to pre­vent a sec­ond wave packet com­ing in from the far right, Ai and Bi must ap­pear to­gether in the com­bi­na­tion ${\rm {Bi}}+{\rm i}{\rm {Ai}}$ as shown in fig­ure A.10. The fact that no sec­ond packet comes in from the far right in the an­i­ma­tion can be taken as an ex­per­i­men­tal con­fir­ma­tion of that re­sult, so there seems lit­tle jus­ti­fi­ca­tion to go over the messy ar­gu­ment.

To com­plete the de­ter­mi­na­tion of the eigen­func­tion for a given value of $E$, the con­stants $C^{\rm {l}}_{\rm {f}}$, $C^{\rm {l}}_{\rm {b}}$ and $C^{\rm {r}}$ must still be de­ter­mined. That goes as fol­lows. For now, as­sume that $C^{\rm {r}}$ has the pro­vi­sional value $c^{\rm {r}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1. Then pro­vi­sional val­ues $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$ for the other two con­stants may be found from the re­quire­ments that the left and right re­gions give the same val­ues for $\psi_E$ and ${\rm d}\psi_E$$\raisebox{.5pt}{$/$}$​${\rm d}{x}$ at the point A in fig­ure A.10 where they meet:

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

That is equiv­a­lent to two equa­tions for the two con­stants $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$, since every­thing else can be eval­u­ated, us­ing the men­tioned soft­ware. So $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$ can be found from solv­ing these two equa­tions.

As the fi­nal step, it is de­sir­able to nor­mal­ize the eigen­func­tion $\psi_E$ so that $C^{\rm {l}}_{\rm {f}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1. To do so, the en­tire pro­vi­sional eigen­func­tion can be di­vided by $c^{\rm {l}}_{\rm {f}}$, giv­ing $C^{\rm {l}}_{\rm {b}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ $c^{\rm {l}}_{\rm {b}}$$\raisebox{.5pt}{$/$}$​$c^{\rm {l}}_{\rm {f}}$ and $C^{\rm {r}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ $c^{\rm {r}}$$\raisebox{.5pt}{$/$}$​$c^{\rm {l}}_{\rm {f}}$. The en­ergy eigen­func­tion has now been found. And since $C^{\rm {l}}_{\rm {f}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1, the $e^{{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$ term is ex­actly the same as the free space en­ergy eigen­func­tion of the first ex­am­ple. That means that if the eigen­func­tions $\psi_E$ are com­bined into a wave packet in the same way as in the free space case, (A.206) with $p$ re­placed by $p_{\rm {c}}^{\rm {l}}$, the $e^{{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$ terms pro­duce the ex­act same wave packet com­ing 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 pro­duce a re­flected wave packet that re­turns to­ward the far left. Note that $e^{-{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$ is the com­plex con­ju­gate of $e^{{\rm i}{p}_{\rm {c}}^{\rm {l}}x/\hbar}$, and it can be seen from the un­steady Schrö­din­ger equa­tion that if the com­plex con­ju­gate of a wave func­tion is taken, it pro­duces a re­ver­sal of time. Wave pack­ets com­ing in from the far left at large neg­a­tive times be­come wave pack­ets leav­ing to­ward the far left at large pos­i­tive times. How­ever, the con­stant $C^{\rm {l}}_{\rm {b}}$ turns out to be very small in this case, so there is lit­tle re­flec­tion.

Fig­ure A.11: Ex­am­ple en­ergy eigen­func­tion for a par­ti­cle en­ter­ing a con­stant de­cel­er­at­ing force field.

Next con­sider the an­i­ma­tion of chap­ter 7.11.3, where the par­ti­cle is turned back by an up­ward po­ten­tial en­ergy ramp. A typ­i­cal en­ergy eigen­func­tion for this case is shown in fig­ure A.11. Un­like in the pre­vi­ous ex­am­ple, where the ar­gu­ment $\overline{x}$ of the Airy func­tions was neg­a­tive at the far right, here it is pos­i­tive. Ta­ble books that cover the Airy func­tions will tell you that the Airy func­tion Bi blows up very strongly with in­creas­ing pos­i­tive ar­gu­ment $\overline{x}$. There­fore, if the so­lu­tion in the right hand re­gion would in­volve any amount of Bi, it would lo­cate the par­ti­cle at in­fi­nite $x$ for all times. For a par­ti­cle not at in­fin­ity, the so­lu­tion in the right hand re­gion can only in­volve the Airy func­tion Ai. That func­tion de­cays rapidly with pos­i­tive ar­gu­ment $\overline{x}$, as seen in fig­ure A.11.

The fur­ther de­ter­mi­na­tion of the en­ergy eigen­func­tions pro­ceeds along the same lines as in the pre­vi­ous ex­am­ple: give $C^{\rm {r}}$ a pro­vi­sional value $c^{\rm {r}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1, then com­pute $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$ from the re­quire­ments that the left and right re­gions pro­duce the same val­ues for $\psi$ and ${\rm d}\psi$$\raisebox{.5pt}{$/$}$​${\rm d}{x}$ at the point A where they meet. Fi­nally di­vide the eigen­func­tion by $c^{\rm {l}}_{\rm {f}}$. The big dif­fer­ence is that now $C^{\rm {l}}_{\rm {b}}$ is no longer small; $C^{\rm {l}}_{\rm {b}}$ turns out to be of unit mag­ni­tude just like $C^{\rm {l}}_{\rm {f}}$. It means that the in­com­ing wave packet is re­flected back com­pletely.

Fig­ure A.12: Ex­am­ple en­ergy eigen­func­tion for the har­monic os­cil­la­tor.

For the har­monic os­cil­la­tor of chap­ter 7.11.4, the analy­sis is some­what dif­fer­ent. In par­tic­u­lar, chap­ter 4.1.2 showed that the en­ergy lev­els of the one-di­men­sion­al har­monic os­cil­la­tor are dis­crete,

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

so that un­like the mo­tions just dis­cussed, the so­lu­tion of the Schrö­din­ger equa­tion is a sum, rather than the in­te­gral (A.206),

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

How­ever, for large $n$ the dif­fer­ence be­tween sum­ma­tion and in­te­gra­tion is small.

Also, while the en­ergy eigen­func­tions $h_n(x)$ are not ex­po­nen­tials as for the free par­ti­cle, for large $n$ they can be pair­wise com­bined to ap­prox­i­mate such ex­po­nen­tials. For ex­am­ple, eigen­func­tion $h_{50}$, shown in fig­ure A.12, be­haves near the cen­ter point much like a co­sine if you scale it prop­erly. Sim­i­larly, $h_{51}$ be­haves much like a sine. A co­sine plus ${\rm i}$ times a sine gives an ex­po­nen­tial, ac­cord­ing to the Euler for­mula (2.5). Cre­ate sim­i­lar ex­po­nen­tial com­bi­na­tions of eigen­func­tions with even and odd val­ues of $n$ for a range of $n$ val­ues, and there are the ap­prox­i­mate ex­po­nen­tials that al­low you to cre­ate a wave packet that is at the cen­ter point at time $t$ $\vphantom0\raisebox{1.5pt}{$=$}$ 0. In the an­i­ma­tion, the range of $n$ val­ues was cen­tered around $n$ $\vphantom0\raisebox{1.5pt}{$=$}$ 50, mak­ing the nom­i­nal en­ergy hun­dred times the ground state en­ergy. The ex­po­nen­tials de­gen­er­ate over time, since their com­po­nent eigen­func­tions have slightly dif­fer­ent en­ergy, hence time evo­lu­tion. That ex­plains why af­ter some time, the wave packet can re­turn to the cen­ter point go­ing the other way.

Fig­ure A.13: Ex­am­ple en­ergy eigen­func­tion for a par­ti­cle en­coun­ter­ing a brief ac­cel­er­at­ing force.

For the par­ti­cle of chap­ter 7.12.1 that en­coun­ters a brief ac­cel­er­at­ing force, an ex­am­ple eigen­func­tion looks like fig­ure A.13. In this case, the so­lu­tion in the far right re­gion is sim­i­lar to the one in the far left re­gion. How­ever, there can­not be a term of the form $e^{-{{\rm i}}p_{\rm {c}}^{\rm {r}}x/\hbar}$ in the far right re­gion, be­cause when the eigen­func­tions are com­bined, it would pro­duce an un­wanted wave packet com­ing in from the far right. In the mid­dle re­gion of lin­early vary­ing po­ten­tial, the wave func­tion is again a com­bi­na­tion of the two Airy func­tions. The way to find the con­stants now has an ad­di­tional step. First give the con­stant $C^{\rm {r}}$ of the far right ex­po­nen­tial the pro­vi­sional value $c^{\rm {r}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1 and from that, com­pute pro­vi­sional val­ues $c^{\rm {m}}_{\rm {A}}$ and $c^{\rm {m}}_{\rm {B}}$ by de­mand­ing that the Airy func­tions give the same val­ues for $\psi$ and ${\rm d}\psi$$\raisebox{.5pt}{$/$}$​${\rm d}{x}$ as the far-right ex­po­nen­tial at point B, where they meet. Next com­pute pro­vi­sional val­ues $c^{\rm {l}}_{\rm {f}}$ and $c^{\rm {l}}_{\rm {b}}$ by de­mand­ing that the far-left ex­po­nen­tials give the same val­ues for $\psi$ and ${\rm d}\psi$$\raisebox{.5pt}{$/$}$​${\rm d}{x}$ as the Airy func­tions at point A, where they meet. Fi­nally, di­vide all the con­stants by $c^{\rm {l}}_{\rm {f}}$ to make $C^{\rm {l}}_{\rm {f}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ 1.

Fig­ure A.14: Ex­am­ple en­ergy eigen­func­tion for tun­nel­ing through a bar­rier.

For the tun­nel­ing par­ti­cle of chap­ter 7.12.2, an ex­am­ple eigen­func­tion is as shown in fig­ure A.14. In this case, the so­lu­tion in the mid­dle part is not a com­bi­na­tion of Airy func­tions, but of real ex­po­nen­tials. It is es­sen­tially the same so­lu­tion as in the left and right parts, but in the mid­dle re­gion the po­ten­tial en­ergy is greater than the to­tal en­ergy, mak­ing $p_{\rm {c}}^{\rm {m}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\sqrt{2m(E-V_{\rm {m}})}$ an imag­i­nary num­ber. There­fore the ar­gu­ments of the ex­po­nen­tials be­come real when writ­ten in terms of the ab­solute value of the mo­men­tum $\vert p_{\rm {c}}^{\rm {m}}\vert$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\sqrt{2m(V_{\rm {m}}-E)}$. The rest of the analy­sis is sim­i­lar to that of the pre­vi­ous ex­am­ple.

Fig­ure A.15: Tun­nel­ing through a delta func­tion bar­rier.

For the par­ti­cle tun­nel­ing through the delta func­tion po­ten­tial in chap­ter 7.12.2, an ex­am­ple en­ergy eigen­func­tion is shown in fig­ure A.15. The po­ten­tial en­ergy in this case is $V$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\nu\delta(x-x_{\rm {A}})$, where $\delta(x-x_{\rm {A}})$ is a spike at point A that in­te­grates to one and the strength $\nu$ is a cho­sen con­stant. In the ex­am­ple, $\nu$ was cho­sen to be $\sqrt{2\hbar^2E_{\rm {nom}}/m}$ with $E_{\rm {nom}}$ the nom­i­nal en­ergy. For that strength, half the wave packet will pass through.

For a delta func­tion po­ten­tial, a mod­i­fi­ca­tion must be made in the analy­sis as used so far. As fig­ure A.15 il­lus­trates, there are kinks in the en­ergy eigen­func­tion at the lo­ca­tion A of the delta func­tion. The left and right ex­pres­sions for the eigen­func­tion do not pre­dict the same value for its de­riv­a­tive ${\rm d}\psi$$\raisebox{.5pt}{$/$}$​${\rm d}{x}$ at point A. To find the dif­fer­ence, in­te­grate the Hamil­ton­ian eigen­value prob­lem from a point a very short dis­tance $\varepsilon$ be­fore point A to a point the same very short dis­tance be­hind it:

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

The in­te­gral in the right hand side is zero be­cause of the van­ish­ingly small in­ter­val of in­te­gra­tion. But the delta func­tion spike in the left hand side in­te­grates to one re­gard­less of the small in­te­gra­tion range, so

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

For van­ish­ingly small $\varepsilon$, ${\rm d}\psi$$\raisebox{.5pt}{$/$}$​${\rm d}{x}$ at $x_{\rm {A}}+\varepsilon$ be­comes what the right hand part of the eigen­func­tion gives for ${\rm d}\psi$$\raisebox{.5pt}{$/$}$​${\rm d}{x}$ at $x_{\rm {A}}$, while ${\rm d}\psi$$\raisebox{.5pt}{$/$}$​${\rm d}{x}$ at $x_{\rm {A}}-\varepsilon$ be­comes what the left hand part gives for it. As seen from the above equa­tion, the dif­fer­ence is not zero, but $2m\nu\psi(x_{\rm {A}})$$\raisebox{.5pt}{$/$}$​$\hbar^2$.

So the cor­rect equa­tions for the pro­vi­sional con­stants 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*}$$

Com­pared to the analy­sis as used pre­vi­ously, the dif­fer­ence is the fi­nal term in the sec­ond equa­tion that is added by the delta func­tion.

The re­main­der of this note gives some tech­ni­cal de­tails for if you are ac­tu­ally plan­ning to do your own an­i­ma­tions. It is a good idea to as­sume that the units of mass, length, and time are cho­sen such that $\hbar$ and the nom­i­nal en­ergy are one, while the mass of the par­ti­cle is one-half. That avoids hav­ing to guessti­mate suit­able val­ues for all sorts of very small num­bers. The Hamil­ton­ian eigen­value prob­lem then sim­pli­fies to

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

where the val­ues of $E$ of in­ter­est clus­ter around 1. The nom­i­nal mo­men­tum will be one too. In those units, the length of the plot­ted range was one hun­dred in all but the har­monic os­cil­la­tor case.

It should be noted that to se­lect a good func­tion $c(p)$ in (A.206) is some­what of an art. The sim­plest idea would be to choose $c(p)$ equal to one in some lim­ited range around the nom­i­nal mo­men­tum, and zero else­where, 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 rel­a­tive de­vi­a­tion from the nom­i­nal mo­men­tum be­low which $c(p)$ is nonzero. How­ever, it is know from Fourier analy­sis that the lo­ca­tions where $c(p)$ jumps from one to zero lead to lengthy wave pack­ets when viewed in phys­i­cal space. {D.44}. Func­tions $c(p)$ that do lead to nice com­pact wave pack­ets 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 es­sen­tially the func­tion $c(p)$ used in this study. The typ­i­cal width of the mo­men­tum range was cho­sen to be $r$ $\vphantom0\raisebox{1.5pt}{$=$}$ 0.15, or 15%, by trial and er­ror. How­ever, it is nice if $c(p)$ be­comes not just very small, but ex­actly zero be­yond some point, for one be­cause it cuts down on the num­ber of en­ergy eigen­func­tions that have to be eval­u­ated nu­mer­i­cally. Also, it is nice not to have to worry about the pos­si­bil­ity of $p$ be­ing neg­a­tive in writ­ing en­ergy eigen­func­tions. There­fore, the fi­nal func­tion used was

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

The ac­tual dif­fer­ence in nu­mer­i­cal val­ues is small, but it does make $c(p)$ ex­actly zero for neg­a­tive mo­menta and those greater than twice the nom­i­nal value. Strictly speak­ing, $c(p)$ should still be mul­ti­plied by a con­stant to make the to­tal prob­a­bil­ity of find­ing the par­ti­cle equal to one. But if you do not tell peo­ple what num­bers for $\Psi$ are on the ver­ti­cal axes, you do not need to bother.

In do­ing the nu­mer­i­cal in­te­gra­tions to find $\Psi(x,t)$, note that the mid point and trapez­ium rules of nu­mer­i­cal in­te­gra­tion are ex­po­nen­tially ac­cu­rate un­der the given con­di­tions, so there is prob­a­bly not much mo­ti­va­tion to try more ad­vanced meth­ods. The mid point rule was used.

The an­i­ma­tions in this book used the nu­mer­i­cal im­ple­men­ta­tions daie.f, dbie.f, daide.f, and dbide.f from netlib.org for the Airy func­tions and their first de­riv­a­tives. These of­fer some ba­sic pro­tec­tion against un­der­flow and over­flow by split­ting off an ex­po­nen­tial for pos­i­tive $\overline{x}$. It may be a good idea to check for un­der­flow and over­flow in gen­eral and use 64 bit pre­ci­sion. The ex­am­ples here did.

For the har­monic os­cil­la­tor, the larger the nom­i­nal en­ergy is com­pared to the ground state en­ergy, the more the wave packet can re­sem­ble a sin­gle point com­pared to the lim­its of mo­tion. How­ever, the com­puter pro­gram used to cre­ate the an­i­ma­tion com­puted the eigen­func­tions by eval­u­at­ing the an­a­lyt­i­cal ex­pres­sion given in de­riva­tion {D.12}, and ex­plic­itly eval­u­at­ing the Her­mite poly­no­mi­als is very round-off sen­si­tive. That lim­ited it to a max­i­mum of about hun­dred times the ground state en­ergy when al­low­ing for enough un­cer­tainty to lo­cal­ize the wave packet. Round-off is a gen­eral prob­lem for power se­ries, not just for the Her­mite poly­no­mi­als. If you want to go to higher en­er­gies to get a smaller wave packet, you will want to use a fi­nite dif­fer­ence or fi­nite el­e­ment method to find the eigen­func­tions.

The plot­ting soft­ware used to pro­duce the an­i­ma­tions was a mix­ture of dif­fer­ent pro­grams. There are no doubt much sim­pler and bet­ter ways of do­ing it. In the an­i­ma­tions pre­sented here, first plots were cre­ated of $\Psi$ ver­sus $x$ for a large num­ber of closely spaced times cov­er­ing the du­ra­tion of the an­i­ma­tion. These plots were con­verted to gifs us­ing a mix­ture of per­sonal soft­ware, netpbm, and ghostview. The gifs were then com­bined into a sin­gle movie us­ing gif­si­cle.