Using an iterative equation solver, and the quantum-mechanical foundations of the problem

This article was prompted by learning how to solve transcendental equations numerically by iteration, using the Hewlett-Packard HP-48G calculator. A good example was necessary, and I chose the one-dimensional particle in a square well potential. This is an interesting problem and a good introduction to quantum mechanics, so the first part of this article is devoted to a discussion of it. Readers interested only in how to solve transcendental equations with the HP-48G can skip to that section immediately.

The "particle" that we will use as an example is the electron. The word "particle" here only means that it will be treated as a mass point in studying its motion, with a position x that is a function of time. For simplicity, we will study only one-dimensional motion, but the results are more general, as is often the case. In Newtonian mechanics, we would write the equation of motion mx" = f, where x" is the second derivative of x with respect to t, m is the mass of the particle, and f the force acting on it. Although this method gives good results in many cases, in the present problem it fails seriously to describe what is observed. The behavior of very small particles is properly described by quantum mechanics, and we shall see that this gives the correct result in the present case.

The electron has a very small mass, m = 9.109 x 10^{-31} kg, which indicates that it does not have strong interactions. In fact, the electron interacts with the electromagnetic and "weak" fields, which are closely connected at a fundamental level. Its electric charge is -e, where e = 1.602 x 10^{-19} C, large enough in comparison with its mass to make electromagnetic forces very effective in moving it around. The gravitational force on its mass, however, is negligibly small. Its interaction with the weak field is responsible for its creation in beta decay, but is negligible in the present problem. The electron also has an angular momentum and a magnetic dipole moment, which are also not important in the present case. The structure of the electron has been investigated since its discovery, but it cannot be described in macroscopic terms. Dirac's relativistic theory of electrons gives the best description of the electron, showing the relation of the electron to the positron, its *antiparticle* with positive charge, and explaining its angular momentum (spin) and magnetic moment. These properties are refined by quantum electrodynamics and the theory of the electroweak interaction, so nearly every detail is known and explained, except for the lack of a macroscopic picture. About the best that can be done is that the electron is a charged point in constant jittering motion at the speed of light within the region where its charge density is appreciable ("Zitterbewegung").

Our problem can be explained physically if we conceive of the electric double layer, two plane sheets of uniform opposite charge density very close together. The field outside the sheets is zero, but between them the electric field is very strong. If σ is the charge density, and s the separation of the sheets, then the electrostatic potential difference between the sheets is Δφ = 4πσs, higher on the positive side and lower on the negative. As in the case of defining the dipole moment, we can set σs = p, the *strength* of the double layer, and assume that s is smaller than distances of interest in the problem. The mechanical potential energy of an electron will change by V = -eΔφ on crossing a double layer. Let us now consider two double layers, one at x = -a and the other at x = +a, with the negative sides outward. If an electron approaches from either direction, it is accelerated by a sudden jerk as it crosses a layer into the region between them, its kinetic energy increasing by V. When it leaves, it is decelerated by the layer, losing kinetic energy V. The region between the two layers is called a *square well potential*. This is not just a mathematical abstraction: such double layers keep electrons inside metals.

We can finally state the problem. We have a particle of mass m that moves in one dimension, x, under the action of no forces, except for those caused by a square well potential of depth V between x = -a and x = a. We want to find solutions in which the electron is bound to the square well, and cannot escape, because its kinetic energy inside the well is less than V. It is convenient to make the potential energy inside the well zero, and V outside, so that all energies in the problem will be positive, and measured with respect to the floor of the well. This situation is shown in the figure at the right. W is total energy, E is the total energy of the particle inside the well, which is also equal to its kinetic energy. V is the height of the well. If E > V, then the kinetic energy would be positive outside the well, and the particle could move freely, only accelerating while crossing the well. For E < V, the particle can only move (classically) between -a and +a, where its kinetic energy is positive. Classically, E can have any value from 0 up to V with the particle confined to the well, but this is not observed. What is found is that only a small number of energies E are possible, sometimes only one or two, and perhaps none at all, so the particle cannot be confined to the well.

The quantities of interest here are the position, x, and the momentum, p = mv. They have the property that xp - px = ih', where h' is Planck's constant divided by 2π, h' = h/2π, called h-bar familiarly. This extraordinary relation is called the *commutator*, expressing everything that is different about quantum mechanics. We chose p instead of v so that this relation would not contain m. The nonzero commutator means that x and p cannot be ordinary numbers, since ordinary numbers commute. Also, the imaginary unit i = √-1 appears, showing that x and p must be complex quantities, in general. Since x and p are now esoteric quantities, we must also devise some way of extracting the usual position and momentum, which are ordinary numbers, from them.

Lack of commutation is a familiar property of square matrices, which suggests that x and p should be represented by matrices, which would represent *operators*. Heisenberg followed up this line, and determined that the matrices are infinite-dimensional, and the usual mechanical quantities are diagonal elements of these matrices. This approach succeeds, but the computational difficulties are great. Nevertheless, Heisenberg was able to show that ΔxΔp = h/2 was a consequence of the commutation relation, where a Δ was a root-mean-square time average of a quantity in the motion. This *Heisenberg Uncertainty Principle* can be used to explain qualitatively many quantum effects. The word "uncertainty" is rather unfortunate, since there is really no uncertainty, just a different way of expressing motion.

Another route was followed by Schrödinger, who noted that (d/dx)xf(x) = x(d/dx)f(x) + f(x), so that x(d/dx) - (d/dx)x = -1 when these differential operators operated on a function of x. If we let the operator p = -ih'(d/dx), then we have xp - px = ih', as required. f(x) is called the *wave function*, which must be a smooth function of x (the function and its first derivative must be continuous) and normalizable, or ∫(-∞,+∞)f*(x)f(x)dx = 1. The average value of any operator O in a motion described by f(x) is o' = ∫(-∞,+∞)f*(x)Of(x)dx.

In general, then, f(x) is a complex function that completely describes the motion. We are leaving out an important factor here, time dependence. In general, the wave function is f(x,t) and its time dependence is determined by the equation ih'(∂f/∂t) = Ef, where E is the total energy. If E is constant, then the time part of f(x,t) is just exp(iEt/h'), so f(x,t) splits into an x factor, f(x), and a time factor. All the states that we will consider here have constant energy E, so they can be described by a wave function f(x) that is a function of x only. Because of the absence of t, these states are called "stationary" where the stationary is meant in a statistical sense, that averages are independent of time (since the exponential time factors multiply to 1). This does not mean that the particle is not moving! The solution of time-dependent problems is so difficult that they are cast into the form of time-independent problems when possible, as in scattering theory. This is the reason for the ubiquity of the Schrödinger wave function, which can be found with the powerful methods for the solution of differential equations.

The energy operator, or *Hamiltonian*, is H(x,p) = p^{2}/2m + V(x), where V(x) is the potential energy as a function of x. We must have H(x,p)f(x) = Ef(x) for a stationary state. Expressing H as an operator by means of the p operator, we obtain the famous *Schrödinger Equation*: (h'^{2}/2m)f"(x) - V(x)f(x) = -Ef(x), or f" + [(2m/h'^{2})(E - V(x))]f = 0. Note that p^{2} must be interpreted as p*p. In the current problem, the wave function outside the well is f" - [(2m/h'^{2})(V - E)]f = 0, and inside the well it will be f" + [(2m/h'^{2})E]f = 0. Both of these equations are easy to solve, since they are second order equations with constant coefficients, and have solutions that are exponentials. If we set α^{2} = 2mE/h'^{2} and β^{2} = 2m(V - E)/h'^{2}, the solutions outside can be expressed as linear combinations of e^{±βx}, while those inside are linear combinations of sin(αx) and cos(αx). The oscillatory wave functions inside the well, where the kinetic energy is positive, have exponential tails outside the well, where the kinetic energy is negative. In a similar problem, the well becomes a wall of height V, which particles with kinetic energy E < V cannot penetrate classically. However, it is found that some particles do get through anyway. This "tunnelling" is a quantum-mechanical effect that is often observed, and is quite real.

The symmetry of the problem in x means that only e^{βx} is valid for x < 0, and only e^{-βx} for x > 0, or else the wave function would grow without limit and not be normalizable. Inside the well, both the cosine and sine solutions exist, but f(x) will include only one or the other exclusively. At x = ±a, the inside and outside functions must join smoothly, with f and df/dx continuous. The details are quite boring, so only the result will be quoted here. The equations are: for the cosine solution, α tan(αa) = β, and for the sine solution, α cot(αa) = -β.. These are the transcendental equations we shall solve to find the energy E, since both α and β are functions of E. I'll come back to them later, after some further observations on the problem.

The *infinite square well*, in which V >> E, is particularly simple to solve. Now all we have are the interior solutions, which must vanish for x = ±a. The derivative f' need not vanish there, since it is equal to the slope of an exponential that approaches zero as V becomes large. Therefore, sin(αa) = 0 and cos(αa) = 0, which together say that αa = &pi/2, π, 3π/2, ..., nπ/2, cosine and sine functions alternating. The lowest state, n = 1, is a single cosine arch. The next state, n = 2, is a single wavelength of the sine, with a zero, or node, at x = 0. State n has n - 1 nodes. The number n is a *quantum number*, and this is typically the way one arises from a problem, as a natural consequence.

The definition of α gives the energy. (αa)^{2} = (nπ/2)^{2} = 2mE/h'^{2}, so E = (πnh')^{2}/8ma^{2}. The energy is proportional to n^{2}, and inversely proportional to ma^{2}. This gives the interesting result that large energy is required to confine a light particle to small dimensions. An electron requires lots of room, so, in a classical sense, it is "big," while the much more massive proton, is "small." In a hydrogen atom, most of the room is occupied by the electron. This is true of all matter, whose apparent size is mainly space occupied by electrons. When you look at a macroscopic body, you are looking mainly at electrons.

The quantity E' = h'^{2}/2ma^{2} is a good reference energy to reflect these things. Let's suppose we have an electron restricted to a distance 2a = 0.1 nm, about the size of a hydrogen atom. The reference energy is 2.4417 x 10^{-18} J, or 15.24 eV, not far from the ionization potential of hydrogen (13.5 eV). For the infinite square well, E = (nπ/2)^{2} E'. The lowest energy is not zero, but (π/2)^{2}E' = 2.47E'. This is called the *zero-point energy*, and is quite real.

In the uncertainty relation, let the rms value of p be Δp, and 2a the rms value of x. Then, 2ap = h/2, or p = h/4a, so the energy is E = h^{2}/32ma^{2} = (2π)^{2}h'^{2}/32ma^{2} = π^{2}h'^{2}/8ma^{2}, exactly the result that we found from the Schrödinger equation. If we assume that for a quantum number n that ΔxΔp = nh/2 (greater uncertainty in the higher states), then we get the energy formula for the square well. The zero-point energy is a result of the uncertainty relation (really a result of the noncommutation of x and p, not of any gambling).

Now we return to the transcendental equations for finding the energy E for a finite one-dimensional square well. In terms of the quantities α^{2} = 2mE/h'^{2} and β^{2} = 2m(V - E)/h'^{2}, the equations are: for the cosine solution, α tan(αa) = β, and for the sine solution, α cot(αa) = -β. We shall see that cosine and sine solutions alternate, and that the number of nodes is n - 1, where n is the serial number of the levels counting from the bottom.

It is convenient to define a new variable x = αa = √(2ma^{2}E/h'^{2}), so that when we have determined x, then E = x^{2}(h'^{2}/2ma^{2}). Then (βa)^{2} = (2ma^{2}V/h'^{2}) - x^{2}. Let r^{2} = V/(h'^{2}/2ma^{2}) be the well depth in units of h'^{2}/2ma^{2}. Then, βa = √(r^{2} - x^{2}). This allows us to express the transcendental equations in terms of the variable x and the parameter r as: x tan x = √(r^{2} - x^{2}) and x cot x = -√(r^{2} - x^{2}). These equations are now ready for solution.

First, create variables X and R, with any values in them. The SOLVE routine will accept an equation, as written above, or an expression, when one term is transposed so that the expression gives the difference between them, or a program that calculates the difference. The easiest seems to be the expression. Create, name and store the two expressions. The Equation Writer is a suitable environment, or they can be created in the command line between the identifiers ' '. One is 'x*tan(x) - √(r^{2} - x^{2})', and the other is 'x/tan(x) + √(r^{2} - x^{2})'. Name the first F1, the second F2, and store them. Set the calculator to RAD mode by pressing ← RAD if necessary.

Now go to → SOLVE. The CHOOS box will have the cursor at "Solve equation" by default. Accept this by pressing the menu key OK. With the prompt on the EQ: line, press CHOOS, select F1 from the list, and press OK. The next two lines should be labelled R: and X:, corresponding to the variables in the equation. Variables of those names should also be available for the program. Move the cursor to R:, key in 1 and press OK. This value will be seen after R:. Now move the cursor to X:. The program will solve for the variable on which the cursor lies when you press SOLVE. You can input a starting trial value for the unknown variable. This value should be near the root to be found, so the iterative procedure moves toward the proper root. In this case, a value of 0 will do. Now press SOLVE, and after a while the value 0.739085133215 will appear, to the full accuracy of the calculator.

CANCEL will now take you to the stack, where X will be displayed in level 1, labelled and in the current display format. The energy E, in terms of V, is E = (x/r)^{2}V. In general, V = r^{2}h'^{2}/2ma^{2} in joule. Here, we have E = 0.5462V (to see all 12 places, press ← EDIT, and CANCEL to get back). This is the only bound state for a well of this depth. For an electron in an 0.1 nm well, V = 15.24 eV and the energy level is at 8.32 eV.

A graphical solution of the problem is easy, and is shown in Schiff (see References). The values of X corresponding to the states are the intersections of a circle of radius R with the curves X tan(X) and -X/tan(X). These curves cross the X-axis at 0, and at values of X such that tan(X) = 0 or tan(X) = ∞. These values are X = 0, π/2, π, 3π/2, ..., where the S-shaped tangent curves cross the axis. Until R reaches these values, there can be no intersection with the corresponding tangent curve, and no solution. For R = 1, the only intersection can be with the tangent curve coming from X = 0, and so there is just one state. When R reaches π/2 = 1.5708, the circle just intersects on the X-axis, and there is a state at E = V, in addition to the one from the intersection with the first tangent curve. Now there are two stable states in the well. At R = π, a third state appears. These are alternately cosine and sine wave functions. As the well becomes deeper and deeper, more states are added.

The reader is encouraged to make a plot up to X = 10, labelling each curve as belonging to a sine wavefunction or a cosine. Then, drawing circles with a compass at intervals up to R = 10, estimate the values of X for each state, and use them as first guesses in the SOLVE routine. After finding the accurate values of X, plot curves of energy versus well depth in eV. the HP-48G makes this very easy. If a manual iterative process is used instead, calculations with trig tables and logarithms are exceedingly tedious, and it is not much better if a calculator is available, unless it is one like the HP-48G with a SOLVE facility.

This theory also applies in three dimensions for a spherical potential well of radius a for states of zero angular momentum, since the radial problem is then a one-dimensional problem of the same nature that can be solved the same way (except for r = 0 to ∞). For details, see Schiff.

L. I. Schiff, *Quantum Mechanics*, 3rd ed. (New York: McGraw-Hill, 1968). pp. 37-42.

*HP 48G Series User's Guide* (Corvallis, OR: Hewlett-Packard Co., 1993). Part No. 00048-90126. Chapter 18.

Return to Physics Index

Composed by J. B. Calvert

Created 5 May 2003

Last revised 6 May 2003