Computer-Generated Motion Pictures of One-Dimensional Quantum Mechanical Transmission and Reflection Phenomena

Goldberg, A, Schey, H M, Schwartz, J L

March 1967

American Journal of Physics

Goldberg and Schey, Lawrence Radiation Laboratory. Schwartz, Science Teaching Center, MIT. Scbey seconded to the Teaching Center for this activity.

ABSTRACT

We describe the details involved in presenting the time development of one-dimensional quantum-mechanical systems in the form of computer-generated motion pictures intended for pedagogic purposes. Concentrating on reflection-transmission phenomena, we formulate the problem in terms of a Gaussian wave packet impinging on a square well or barrier and being reflected and transmitted. The wave equation is solved numerically by methods discussed in detail and photographs of the wave packet vs position at a variety of times and for a range of projectile energies are given.

I. INTRODUCTION

THE purpose of the present article is to describe the physical and mathematical details underlying the presentation of the time development of physical systems in the form of motion pictures. It is our hope that such pictures will help provide insight into the behavior of these physical systems and thereby be of some pedagogic value.

We concentrate here upon one-dimensional quantum-mechanical scattering (reflection-transmission) phenomena as described by the time-dependent Schrodinger equation. We formulate the problem in terms of a localized wave packet moving towards a potential barrier or well and being reflected and transmitted. The wave equation is solved numerically by methods to be discussed and the probability density projected on a cathode-ray tube. From the tube, photographs are made, and these in turn are processed into the successive frames of the film.

II. FORMULATION OF THE PROBLEM

The standard textbook presentation of reflection and transmission phenomena treats the problem in a time-independent fashion and requires a solution of Hψ= Eψ subject to the conditions:

See for example L. I. Schiff, Quantum Mechanics (McGraw-Hill Book Co., New York, 1949), pp. 92-95.

where the potential is short ranged and localized near x=0 and the projectile particle is assumed incident from the left. Since reflection and transmission are, in fact, time-dependent phenomena, such a formulation is, in many ways, undesirable especially from a pedagogic point of view. This fact has been noted by a number of authors, and Low, in particular, has given a complete, time-dependent formulation of scattering in three dimensions which provides a convincing justification of the standard time-independent formulation.

F. E. Low (unpublished lecture notes). See also E. Merzbacher, Quantum Mechanics (John Wiley & Sons, Inc., New York, 1961), pp. 217.

For our purposes here, namely, a motion-picture display of reflection and transmission, we clearly must adopt a time-dependent formulation in the manner of Low. Our plan, then, is to choose a one-dimensional wave packet initially centered at some position a with a spread Δx and having an average wave number k₀ with spread Δk and follow its development in the course of time as it moves into the region of a short-range potential and is reflected and transmitted.

This program, then, requires an integration of the time-dependent wave equation for a system governed by a time-independent potential V(x):

For ease of writing, we work in a system of units in which the mass m=½ and h=1. Then Eq. (1) becomes:

The standard procedure for solving this equation is in terms of the stationary states of the Hamiltonian H. If we have Hu_n= E_nu_n for the bound states and Hu_k= E(k)u_k for the continuum, then a formal solution to the problem is given by

where ψ(x,0) is the given initial state of the system. This procedure, however, has not been adopted here for several reasons. To begin with, for almost all potentials the eigenfunctions and eigenvalues have to be determined numerically. Even if we were to quantize in a large finite region (the one-dimensional analog of box normalization), in order to reduce the continuous part of the spectrum to a discrete spectrum, there would obviously be a large number of states involved in Eq. (3). Second, even with all the eigenfunctions and eigenvalues known or evaluated, the computation of the terms comprising Eq. (3) at each x and t of interest is a lengthy task. Finally, we should have to restrict the initial state of the system, if ψ(x,0), to some limited combination of eigenstates or be forced to carry out the many numerical integrations required to evaluate the coefficients a_n and a(k) given by Eq. (4) and Eq. (5) and required in Eq.( 3).

To avoid these problems we have found it convenient to deal directly with Eq. (2) making use of the eigenstate expansion Eq. (3) only for the purpose of checking the accuracy of our solution at a few isolated values of the time t. In dealing with Eq. (2), we convert it from a differential equation to a finite difference equation in both space and time, a procedure which is described in detail in the following section.

III. NUMERICAL INTEGRATION TECHNIQUES

For background and details concerning the subject of this section the reader is referred to R. D. Richtmyer, Difference Methods for Initial-Value Problems (Interscience Publishers, Inc.. New York, 1957).

In order to convert Eq. (2) to a finite difference equation, let us designate time by a superscript n and spatial position by a subscript j. Thus ψ(x,t) → ψ_jⁿ. The various x values become jε where ε is the mesh width and j=0, 1,..,J. Similarly, the time variable goes over into nδ where δ is the time step and n=0, 1,2,....

To derive the difference equation we require, let us begin by considering the spatial variation of the wavefunction ψ_j; for this purpose we suppress the time index n. We then make two Taylor's series expansions as indicated:

Upon adding these equations a little manipulation quickly yields the standard expression

This expression for the second derivative of ψ evaluated at the point x=x_j=jε then replaces ∂²ψ/∂x² in Hψ, and we have

The operation of H on ψ_j is always taken in this way in what follows.

To treat the time development of ψ, we note that the formal solution of Eq. (2) is given by

where the exponential operator is defined in terms of the Taylor's series

which, correct to terms of order δ, becomes

If we now use Eq. (6) to evaluate Hψ_j, this last equation assumes the form

One feature of this equation is especially attractive: it is an explicit differencing scheme, i.e., one which gives the wave function at time step n+l directly in terms of the function at the earlier time n. Clearly this permits a straightforward integration scheme and allows for an easy determination of ψ for all n and j. The scheme, however, has the very serious drawback of being unstable; that is to say, roundoff error and errors arising from dropping terms of order δ² and higher in the expansion of e^-iδH can grow without bounds as the time integration proceeds from step to step.

A procedure which eliminates instability involves rewriting Eq. (8) in the form

which leads to the difference equation

Though stable, this equation is implicit rather than explicit. By this we mean that ψⁿ⁺¹ is not given directly in terms of ψⁿ and certain algebraic techniques must be applied in order to obtain the wave function at time step n+1 from that at n. This is not an insurmountable problem - in fact the differencing scheme we finally settle on is implicit, and we indicate below the way it is handled. For the present we remark that Eq. (11) must be ruled out on the grounds that it is not unitary, a flaw present in Eq. (9) as well.

Unitarity is the characteristic of the original wave equation which ensures that the normalization of the wave function does not change in time. The operators e^±iδH with H Hermitian are clearly unitary. However, no approximation of these operators which involves expanding the exponential and retaining only a finite number of terms is unitary. In particular, the first-order approximations we have been discussing, 1±iδH, are not unitary. Thus, both differencing schemes, one of which leads to Eq. (9), the other to Eq. (11), involve mutilating the the time-development operators e^±iδH in such a way as to destroy an important physical characteristic. A relatively simple unitary approximation to e^±iδH is, however, provided by the Cayley form

which has the further desirable property of being correct to order δ². Let us, then, replace e^-iδH in Eq. (8) by the Cayley form to get

Carrying through the indicated algebra and using Ea. (6) to evaluate Hψ_jⁿ and Hψ_jⁿ⁺¹ yields the difference equation

where λ= 2ε²/δ. This difference equation is stable and unitary, but implicit; and we now discuss the means by which it may be solved.

Let us define

Following standard procedure we make the assumption that

where e_jⁿ and f_jⁿ are two auxiliary functions defined by this equation. Substituting Eq. (15) into Eq. (14), we find

This equation is identical with Eq. (15) exc;ept for its being written for j rather than j-1. By comparing Eqs. (15) and (16), we get

Solving these for e_jⁿ and f_jⁿ in terms of e_j-1ⁿ and f_j-1ⁿ yields

It follows from Eq. (17) that e_jⁿ is time independent, so that we may write simply e_j, there being no need for the time index. The same is not true of f_jⁿ for it depends upon Ω_jⁿ which is a time-dependent quantity.

Equations (17) and (18), which we have derived above, are recursion relations for the e and f functions. To obtain starting values for the recursion, we imagine that the physical system with which we work is situated in a large (one-dimensional) box, and that the wavefunction for the system must vanish on the walls (end points) of the box. Thus

where we assume that the box extends from 0 to L: 0 ≤ x ≤ L. Now j=0 represents x=0 and, we define j =J such that L=J_ε. Then our boundary conditions may be stated in the form

With the first of these conditions in mind, let us write Eq. (14) for the case j= 1:

By comparing this last expression with Eq. (15), we get the results

These starting values, in conjunction with the two recursion relations Eqs. (17) and (18), now provide a means of determining e₁, e₂, e₂, ... and f₁ⁿ, f₂ⁿ, f₂ⁿ ... for all n.

From the second of the two boundary conditions, ψ_Jⁿ=0, all n, and Eq. (15) written for j=J-1, we have

Thus, by recurring up on the e's and f's, we find an expression for the wavefunction at the next-to-last net point j =J-1 (it is zero at j=J). By inverting Eq. (15) we get

Now, by combining Eqs. (19) and (20), we can determine ψ_J-2ⁿ⁺¹; and this, in turn, when used in Eq. (20), provides ψ_J-3ⁿ⁺¹ etc. down to ψ₁ⁿ⁺¹ (ψ₀ⁿ⁺¹ being zero). Thus we have outlined a procedure by means of which one may obtain ψ_jⁿ for all j=0,1,2, ... ,J and all n. One sees that ultimately this procedure must begin with a knowledge of ψ_j^o for all j; i.e., of the initial wavefunction, the one conventionally written ψ(x,0), all x, for this enables us to evaluate Ω_j⁰ from which all the rest follows.

This completes the specification of the method of solution, and we now derive the criteria to be satisfied by a suitable choice of input parameters.

IV. CRITERIA FOR CHOICE OF INITIAL CONDITIONS

We treat the reflection-transmission event in a "physical" way by following the time development of a wave packet as it moves into and out of the region of a potential. This initial wave packet can, of course, be chosen in a variety of ways, the choice being dictated largely by convenience. We have elected to represent the initial state of the particle which impinges upon the potential by a Gaussian wave packet. Thus, we put

We see that this packet is centered about x=x₀ with a spread in x governed by σ₀. The factor e^ik₀x makes our initial wave function move to the right with average momentum k₀. Since we are quantizing in a large box of length L at the edges of which the wave function must vanish, it is clear that x₀ and σ₀ must be chosen so that ψ(0,0) and ψ(L,0) are essentially zero, Choosing x₀ = L/4, we find

If we choose σ₀ to be 5% of the box length L, σ₀=L/20, we find

To approximate this by zero is clearly to incur no perceptible error, and the situation is obviously even more favorable for |ψ(L,0)|. The choice σ₀=L/20 in the initial Gaussian wave packet is thus compatible with box normalization, but there are two further restrictions which must be imposed in order that box normalization not give rise to difficulties. First of all, the packet must not be allowed to travel so far that it hits the walls. We ensure this by letting the center of the packet which, as already mentioned, starts from x=L/4, move no farther to the right than x=3L/4. This, in turn, is accomplished by requiring that the average velocity of the wave packet be

where T= Nδ is the total time required by the packet to move from L/4 to 3L/4. Thus

where we have used L= J_ε and λ= 2ε²/δ. With this choice of k_- and with x₀ = L/4, we utilize approximately only the center half of the available space, omitting the outer quarters, and the packet is thus prevented from reaching the walls.

The second of the two restrictions concerns the fact that the wave packet spreads in the course of time. We must arrange matters so that the reflected and transmitted packets continue to be well enough localized so that at the end of the event they are out of the region of the potential and still far from the walls. An analysis of the spreading of a Gaussian wave packet, which is a standard problem, shows that if σ₀ is the initial spread, then after a time T the spread is given by

Actually this expression applies to a Gaussian packet moving in free space, but our actual numerical calculations show that it can provide reliable estimates of the spreading even in the presence of a potential. For negligible spreading we must therefore have 16N²ε⁴/λ ² reasonably small compared with σ₀.

E. Merzbacher, p159

There are additional limitations on the choice of input parameters. We consider first those made necessary by the replacement of the continuous variable x by the discrete index j as discussed in Sec. III. To do this we examine two extreme cases: (1) propagation with no potential at all, and, (2) propagation in a constant well set equal to the maximum value of the potential in the actual problem. We first note that since the mesh width ε is the smallest inverval in x which can be resolved by the differencing method, the largest wave number allowed is k_max =π/ε. With the initial Gaussian wave packet as given by Eq. (21), the probability distribution in momentum space is

Let k_m be the largest momentum present to any appreciable extent (e.g. with probability 10^-5 ) in this distribution. Then, clearly, we must require that k_m< k_max.

In the first case, with no potential, the eigen-states of the system are φ_k(x) = sin kx; the second derivative of φ_k is, of course,

Since k_ε is presumed to be small, this may be approximated by

Thus, in order that (24) be an accurate approximation to Eq. (23) we must require that

In the second case, where the potential is constant with the value V_max, the eigenstates are φ_k(x) =sin(k²+ V_max)^½x. Thus, an argument like the one just given results in the requirement

We turn now to the effects of the finite time step δ. Let us consider the actual system to be treated - that is, the large one-dimensional box with a potential in the center. The states of this system are denoted by χ_k(x,t) =e^-iω_ktT_k(x). Note that outside the region of the potential χ_k(x,t) is of the order e^±ikx-ω_kt and ω_k=k². If the errors introduced by the spatial mesh ε are ignored, the finite differencing in time, Eqs. (10) and (12), gives

From Eqs. (26) and (27) there result two important consequences. First, the over-all normalization of each mode is constant in time; so those modes which were negligible at t=0 remain negligible. This, of course, is a consequence of the unitarity of approximation (12). Second, at every time step a phase error φ is introduced into each mode of the system; by comparing Eqs. (26) and (27) we find

The complete wavefunction ψ(x,t) is a coherent sum of the eigenstates with average energy ω₀=k₀². Thus, in the probability density

it is only the relative phase error which will enter, and to minimize it we require that

In summary then, the following conditions must be satisfied in constructing a one-dimensional scattering event:

Of course, the last four conditions, being inequalities, do not tell us precisely what values of ε and δ are needed to give accurate results. Therefore, we have constructed a scattering event in the above manner with a square-well potential utilizing the least favorable values of well depth and momentum k₀ of the set we have considered and choosing several values of δ and ε (or alternatively λ and ε). We have also analyzed the scattering event by the standard expansion in eigenstates and used this exact solution to check the accuracy of our numerical integration procedures. In this way we have chosen the quantities ε and λ to be used in the actual computation for film production. The specific quantitative results are contained in the next section.

As we have remarked, the calculation of the solution in terms of an eigenstate expansion is much too time consuming to be used in the production of a film. The eigenstate expansion, evaluated at a few time steps, was used only as a check to determine the accuracy of the much more rapid numerical-integration method.

V. INPUT PARAMETERS

There is no loss in generality in choosing L = J_ε=1, that is to say, in normalizing in an interval of unit length. Therefore, from the above discussion we have σ₀=0.05, and the initial position of the center of the wave packet is x₀ = ¼. Trial computation and comparison with the exact result for a square-well potential [Eqs. 3), (4) and (5)] indicates that if we choose

the spatial differencing procedure is quite accurate. Since the distribution of momenta in the packet is Gaussian and sharply peaked, the largest effective momentum k_m is not too different from k₀, and hence criterion 3 is well satisfied. Apart from an exception to be noted later, the limitations on computer speed and memory make it convenient to restrict the number of mesh points to, at most, 2000; and thus we can treat momenta k₀≤100π. We have then made computations at three values of k₀, choosing the mesh interval to satisfy Eq. (28).

Thus:

This gives a range of energy which should satisfactorily demonstrate the energy dependence of one-dimensional scattering phenomena.

We have also chosen a value of δ (or rather λ) by selecting a series of values and comparing the results of the numerical integration with the exact results. We note from Eq. (22) that the distribution of momenta is very strongly peaked about k= k₀. Since the momentum resolution of the packet is Δk/k₀∼1/σ₀k₀, the most poorly resolved distribution is that with k₀ equal to its smallest value, 50π. In this worst case, the density of momenta falls to 10^-5 of its maximum values for k_m=1.5k₀. Thus, if we choose δ so that many time steps are contained in one period of the dominant mode, k=k₀, all other modes in the packet will also be treated accurately. The value λ=1 gives a time step δ=2ε², two hundred and fifty times smaller than the period of the dominant mode. That this is adequate is verified by comparison of the results of the integration with the exact solution. Criterion 1 then determines the number of time steps required to complete a scattering event. For the above momenta, between 800 and 1600 steps are needed, a rather moderate number of computations.

We may also note that with the above parameters the width of the packet at the end of the scattering event is always less than 1.65 σ₀, satisfying criterion 2.

Finally, we must choose a potential. For this film and for the purpose of checking the accuracy of our integration method, we have used square wells and square barriers with depth (or height) V₀= ±2(50π)². In all cases, criterion 4 is satisfied insofar as 0.004 is small as compared with unity. We have chosen the half-width of the potential to be 0.032, guaranteeing that there is essentially no spatial overlap between the wave packet and the potential both at the beginning and at the end of the scattering event.

We have made a quantitative comparison between the results of the numerical integration and the exact solution for the problem of transmission and reflection by a square well, with the parameters k₀= 50π, ε= 10^-3, λ=1. The two calculations give results for the reflection and transmission coefficients which agree to within 2%, while the positions of the peaks and valleys in the distribution functions |ψ(x,t)|² (see Figs. 1-6) are virtually identical. We, therefore, feel that films generated in this manner are accurate, both qualitatively and quantitively, and can be pedagogically useful.

VI. RESULTS AND DISCUSSION

The film we have produced by the methods discussed here consist of six scattering events each lasting about 1 min. The entire film, including titles, has a running time of about 10 min. Examples of what is seen in the film are shown in Figs. 1-6. The first three sets of diagrams represent reflection and transmission by a potential well, and the second set of three by a barrier. For the most part these illustrations speak for themselves. We should emphasize, however, that the rapid oscillations appearing when the packet is near the potential are accurate and do not result from approximations or machine error.

Figure 5, in which the particle energy is equal to the barrier height, requires some special comment, for it represents a type of resonance. A portion of the wave packet is captured by the barrier and remains trapped for a period which is long as compared with the usual time of transmission through the barrier. The captured part of the packet is seen in the diagrams to bounce to and fro between the barrier walls, a small amount of probability escaping at each collision. Eventually, by this mechanism, the entire packet escapes. To generate a film of this event, it was necessary to alter some of the input parameters of Sec. V. Specifically, we chose J = 5000 and ε=0.7071×10^-3 so that L=Jε ∼ 3.5 rather than 1 as in all other cases. (Note that the large number of net points here would not be feasible in all cases since it increases computation time by a factor of about five.) Further, only the center portion of the region between x=0 and x=L was photographed, enabling the transmitted and reflected parts of the packet to leave the field of view so that the nonphysical reflection from the walls (which takes place too far from the region of the barrier to influence what occurs there) will not appear in the film to interfere with a clear presentation of the escape of the trapped part of the packet.

The reader will note occasional breaks or discontinuities ar various points in some of these pictures. This effect is due to an evidently inherent malfunction of the equipment involved in rendering the machine calculations into graphical form. At the present time there appears to be no simple way to avoid this problem; but, fortunately, the discontinuities, while annoying, by no means destroy the effect of the film. This situation emphasizes the fact that the use of computers to illustrate time development in physical systems by motion pictures is still in a preliminary, if no longer rudimentary, stage. Our feeling, nonetheless, is that despite their shortcomings these films already have considerable merit as a pedagogical tool and, as methods are refined and scope broadened, they will play an increasingly important role in science teaching both at the college and graduate levels.

ACKNOWLEDGMENTS

The authors are deeply grateful to Raymond L. Tarp, who expertly programmed all computations and plotting routines involved in this project. We wish also to thank Dr. Cecil E. Leith for an important conversation. Computations for this work were done on the Livermore IBM-7094 and CDC-3600.

COMPUTER MOVIES