On the connection between bound and scattering states of finite square-well potentials: a unified approach

We start by writing the general solution of the time-independent Schrödinger equation for a symmetric potential well with depth V_o and width 2a, as shown in figure 1. For any coordinate interval where the potential has a constant value V₁, the general solution for the time-independent Schrödinger equation is:

$$\begin{equation}{\Psi}\left(x\right)=A\enspace {\text{e}}^{\text{i}{k}_{1}x}+B\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{1}x},\qquad {k}_{1}=\frac{\sqrt{2m\left(E-{V}_{1}\right)}}{\hslash }.\end{equation} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

It is a standard textbook problem to calculate the transmission coefficient for a particle with an energy above the well [1, 2, 4, 7]. In this case, the general solution consists of three parts that can be written in the same form as equation (1), without assuming anything about the relative magnitude of the energy E compared to the value of the potential over each coordinate interval:

$$\begin{equation}{\Psi}\left(x\right)=\begin{cases}A\enspace {\text{e}}^{\text{i}{k}_{1}x}+B\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{1}x},\qquad \quad x{< }-a\hfill \\ C\enspace {\text{e}}^{\text{i}{k}_{2}x}+D\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{2}x},\quad \qquad -a{< }x{< }a\hfill \\ F\enspace {\text{e}}^{\text{i}{k}_{1}x}+\left(G=0\right){\mathrm{e}}^{-\mathrm{i}{k}_{1}x},\quad x{ >}a\hfill \end{cases}\end{equation} \tag{ 2 }$$

where

$$\begin{equation}{k}_{1}=\frac{\sqrt{2mE}}{\hslash }\enspace \mathrm{R}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}\enspace \mathrm{I}\enspace \mathrm{a}\mathrm{n}\mathrm{d}\enspace \mathrm{I}\mathrm{I}\mathrm{I}\end{equation} \tag{ 3 }$$

$$\begin{equation}{k}_{2}=\frac{\sqrt{2m\left(E-{V}_{\mathrm{o}}\right)}}{\hslash }\enspace \left(\mathrm{Re}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{n}\enspace \mathrm{I}\mathrm{I}\right)\end{equation} \tag{ 4 }$$

Since V_o < 0, k₂ will always be real if E is above the bottom of the well. The coefficients A, B, ..., G allow us to construct a variety of solutions subject to the conditions of continuity in Ψ(x) and Ψ'(x) at x = ±a, which provides four equations that relate these coefficients. In keeping with the fact that the Schrödinger equation is a 2nd order ordinary differential equation, one must also impose two boundary conditions at x = ±∞. Imposing the boundary conditions completely defines the coefficients. We chose to examine solutions where G = 0, a choice that ensures that the particle has positive momentum in region III [27]. In many publications and textbooks, this choice is usually described as necessary to ensure that all particles are incident from the left-hand side of the well. While dynamical language such as 'incident from the left' or 'traveling to the right' is commonly used when describing particles in scattering problems, one must remember that we are dealing with the time-independent Schrödinger equation. The wave function defined by equation (2) is a 'steady-state' solution for a particle as it could exist in a time-independent potential. Upon restoring the time-dependence of the wave function, solutions with G = 0 correspond to right-moving plane waves in region III.

Applying the continuity conditions for the wave functions and their derivatives at the boundaries yields the following four equations that relate the coefficients A, B, C, D, and F:

$$\begin{equation}{\Psi}\left(x=-a\right){\Rightarrow}A\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{1}a}+B\enspace {\text{e}}^{\text{i}{k}_{1}a}=C\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{2}a}+D\enspace {\text{e}}^{\text{i}{k}_{2}a}\end{equation} \tag{ 5 }$$

$$\begin{equation}{{\Psi}}^{\prime }\left(x=-a\right){\Rightarrow}A\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{1}a}-B\enspace {\text{e}}^{\text{i}{k}_{1}a}=\left(\frac{{k}_{2}}{{k}_{1}}\right)\left(C\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{2}a}-D\enspace {\text{e}}^{\text{i}{k}_{2}a}\right)\end{equation} \tag{ 6 }$$

$$\begin{equation}{\Psi}\left(x=a\right){\Rightarrow}C\enspace {\text{e}}^{\text{i}{k}_{2}a}+D\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{2}a}=F\enspace {\text{e}}^{\text{i}{k}_{1}a}\end{equation} \tag{ 7 }$$

$$\begin{equation}{{\Psi}}^{\prime }\left(x=a\right){\Rightarrow}C\enspace {\text{e}}^{\text{i}{k}_{2}a}-D\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{2}a}=\left(\frac{{k}_{1}}{{k}_{2}}\right)F\enspace {\text{e}}^{\text{i}{k}_{1}a}\end{equation} \tag{ 8 }$$

When E > 0, it is common to explore the effect of the potential on the particle by its transmission through the potential and its reflection off of the potential, which are quantitatively defined in terms of transmission and reflection coefficients. The transmission coefficient is defined as the ratio of the probability flux (or probability current) of the transmitted wave to the incident wave. Analogously, the reflection coefficient is defined as the ratio of the probability flux of the reflected wave to that of the incident wave. In one dimension, the probability flux for a plane wave, such as that of the incident component of our particle's wave function, is [4, 26]

$$\begin{equation}j=\frac{\hslash {k}_{1}}{m}{\left\vert A\right\vert }^{2}\end{equation} \tag{ 9 }$$

In our case, the transmission coefficient is then given by [28],

$$\begin{equation}T\left(E\right)=\frac{{k}_{1}}{{k}_{1}}{\left\vert \frac{F}{A}\right\vert }^{2}={\left\vert \frac{F}{A}\right\vert }^{2}\end{equation} \tag{ 10 }$$

When E is greater than zero, the wave function cannot be normalized, and it is defined up to an arbitrary constant. One often sets A = 1, but we prefer not to do so in the spirit of keeping this example as general as possible, later allowing us to discuss the role of this coefficient explicitly, so we will let A take on any value except for zero. Solving for T(E) with this condition on A yields the expression:

$$\begin{equation}T\left(E\right)={\left\vert \frac{F}{A}\right\vert }^{2}={\left[1+{\left(\frac{{k}_{1}^{2}-{k}_{2}^{2}}{2{k}_{1}{k}_{2}}\right)}^{2}{\mathrm{sin}}^{2}\left(a{k}_{2}\right)\right]}^{-1}\end{equation} \tag{ 11 }$$

Restricting ourselves to the case that V_o is negative in region II and substituting equations (3) and (4) into equation (11) yields,

$$\begin{equation}T\left(E\right)={\left\vert \frac{F}{A}\right\vert }^{2}={\left[1+\frac{\left\vert {V}_{o}\right\vert }{4E\left(1+\frac{E}{\left\vert {V}_{\mathrm{o}}\right\vert }\right)}\cdot {\mathrm{sin}}^{2}\left(\frac{a}{\hslash }\sqrt{2m\left(\left\vert {V}_{\mathrm{o}}\right\vert +E\right)}\right)\right]}^{-1}\end{equation} \tag{ 12 }$$

Let us first examine this equation with a concrete example to get a sense of the dimensionality of this type of problem before moving onto to more generalized cases. Take an electron incident on a square-well potential of depth V_o = −9 eV and width 2a = 1 nm. There is nothing in the mathematical derivation of equation (12) that would limit it to the case where E > 0. Figure 2(a) shows a plot of equation (12) as a function of the electron's energy from 0 eV to 10 eV. The standard oscillatory behavior of the transmission into region III is observed, reaching T(E) = 1 for the energies corresponding to scattering resonances, where the particle has a 100% chance of making it into region III [1, 2, 23–26]. Figure 2(b) shows a plot of T(E) for E < 0, which has poles at a set of discrete energy values, with the spacing between these values increasing toward higher energies. For our case, the energy values where T(E) diverges are: −8.706 eV, −7.802 eV, −6.371 eV, −4.403 eV, and −2.020 eV. What do the energies where T(E) goes to infinity represent? It can be shown that each of these poles occurs at a bound state energy eigenvalue [13]. Some students do recognize that these values should be the energy eigenvalues, mostly due to the increased spacing at higher energies which they recognize from energy spectrum of the infinite square well, however they typically cannot correctly explain why this is the case. So why do these poles represent bound states when the problem was initially solved assuming that the particle was above the well, and why do the poles take on values greater than 1 for E < 0?

Zoom In Zoom Out Reset image size

Although T(E) is typically interpreted as the 'transmission probability' for energies above the well and has an upper limit of unity due to the conservation of probability flux [4], it does not lose significance when E < 0, however its interpretation requires explanation: in this example, T(E) is technically just a function that defines the ratio between the amplitudes of two solutions to the Schrödinger equation. When V₀ < E < 0, the momenta related to k₁ in the outer regions (I and III) become imaginary, i.e., k₁ → i|k₁|, and so the wave functions in these regions are real exponentials. For the overall solution to remain well-behaved in regions I and III, it must decay exponentially to the left and right of the well, respectively. This occurs when A = G = 0 and F = B ≠ 0. Our boundary condition that G = 0 takes care of region III, but A was assumed to be a non-zero value, and so the component of the solution associated with A will diverge as x → −∞. However, if A were to become zero for V₀ < E < 0, the function T(E) would diverge, and this makes T(E < 0) a handy quantity since the divergences indicate where A has vanished.

To understand how A can become zero, but only at certain energies, let us revisit equation (12). Upon inspection, the function T(E) will go to infinity when the term inside the brackets becomes zero. If we define everything inside of the brackets as β , then equation (12) becomes:

$$\begin{equation}T\left(E\right)={\left\vert \frac{F}{A}\right\vert }^{2}=\frac{1}{\boldsymbol{\beta }}\end{equation} \tag{ 13 }$$

Upon rearranging equation (13), we see that $A=\sqrt{\beta }\left\vert F\right\vert$, and under the condition that β = 0 and F is finite, A vanishes. We never preemptively imposed the physical requirement that A = 0, but this condition is naturally satisfied when β = 0. Consequently, any energy where β = 0 is the energy of an eigenstate, which allows one to go back to equation (2) and obtain the wave function for that eigenstate. It should be noted that the reflection coefficient for this potential would be defined as R(E) = |B/A| [1, 2, 2–5] and while R(E) = 1 − T(E) for E > 0, R(E) diverges at the same energies as T(E) for E < 0.

Using the transmission coefficient to determine the eigenvalues of the square-well demonstrates the difference between solving the Schrödinger equation and picking relevant solutions. Equation (11) was derived under the mathematical restriction that the solutions to equation (2) and their derivatives were continuous across the boundaries. This yielded a continuum of mathematically valid solutions for electron energies less than zero. We never imposed the boundary condition that the wave function remained finite, and therefore most of these solutions are physically meaningless. However, for a specific subset of these solutions, the physical requirement of finiteness is naturally satisfied. In the conventional method, one picks out this subset of solutions when the physical condition that A = 0 and G = 0 is imposed. Essentially, one preemptively limits the solutions to only those that correspond to β = 0, before solving the problem for the bound state energies.

We have established that the poles in T(E) represent bound states, which means that the energy values where these occur can be used in equation (2) to obtain their associated wave functions, but how precisely does one need to know these energies? It turns out that a small numerical error in the energy eigenvalue (even if only in the precision of the value) can still result in solutions to equation (2) that are not well-behaved for all values of x. Figure 3 shows a plot of equation (2) for energies around the pole associated with the lowest bound state (E ∼ −8.706 eV). A divergence exists in region I (x < −0.5 nm) that becomes unacceptably close to the barrier of the potential well when the energy is off from the 'exact' pole by just 0.001 eV or 0.01%. This subtlety is not typically addressed but is critical in making appropriate numerical approximations [13, 14, 23–26]. The function in region I is a superposition of two solutions with coefficients A and B. This divergence originates from the solution associated with A, which is usually set to zero in textbooks when dealing with bound states. However, as we mentioned above, it is unnecessary to set A to zero to get a 'well-behaved' wave function. The amplitude of the divergent component goes to zero naturally, so long as one uses a 'reasonably' precise numerical value for the bound state energy. The fact that only at infinitely precise energies will true bound states exist for static problems is an essential point in quantum mechanics. When performing numerical evaluations, one must be mindful of the precision used and the consequences of that decision.

Zoom In Zoom Out Reset image size

Next, we will follow the methodology introduced by Walker and Gilmore [23–25], but also expand upon it using the treatment described above, which allows us to quickly probe the energies of the bound and scattering states for wells of various depths and widths by simply using equation (12), without any rederivation. Moreover, we can probe the solutions of the function defined by equation (2) and gain useful insight on what happens to these solutions as one approaches the energies corresponding to the poles in T(E). In the following, the term 'wave function' will only be used to describe a solution to equation (2) that also behaves acceptably, i.e., there is no divergence close to the edge of the well (figure 3).

With these considerations in mind, let us investigate T(E) and the associated solutions to equation (2) for a well of fixed width, but increasing depth. In figure 4(a), T(E) is plotted for four different well depths where the number of bound states, as seen from the number of poles, goes from three to four. The green arrow tracks a scattering resonance (for E > 0), which moves closer to the top of the well as the depth of the well is increased. If the plot of T(E) is monitored as the well depth is continuously lowered, it is clear that every new bound state is correlated to a scattering resonance that was 'pulled' into the well (see supplementary material for the movie) [23]. In figure 4(b), the absolute value of equation (2) evaluated at the energy corresponding to the transmission resonance indicated by the arrow for a well depth of −3.32 eV is shown on the right. On the left of figure 4(b), the same is shown for the highest bound state energy at a well depth of ∼−3.39 eV. The part of the solution representing the electron in region II has a similar shape in both cases, except that the bound state wave function has four peaks inside the well. In contrast, for the scattering resonance, the peak of the last oscillation exactly coincides with the boundary of the well.

Zoom In Zoom Out Reset image size

We have just shown that each new bound state that appears when the well's depth is increased has a one-to-one correspondence to a scattering resonance that has been pulled into the well. Another way to think of this is that the well width and particle mass restrict the value of the wavenumbers that are needed to meet the conditions for a scattering resonance or a bound state. Most of the energies where these conditions are met exist as scattering resonances. However, once the well is deepened, the energies where the required wavenumbers exist get closer to the top of the well until they finally exist within the well, at which point the scattering resonance becomes a bound state, which was also pointed out by Gilmore using a different representation [23]. Overall, once the students associate a pole in T(E) as a bound state, and a peak for E > 0 as a resonance, they can directly observe that each pole (bound state) is always a result of a transmission resonance peak being pulled below the well. They can then probe the wave functions above and below the well at the point of crossover and observe, as Gilmore puts it, that scattering resonances near the top of the well are 'almost bound states' [23], which the presented method makes very clear.

Next, we will show that the calculations above can be significantly simplified using a well-known matrix formalism based on the so-called 'transfer matrices'. In essence, transfer matrices generalize the concept of the transmission coefficient, which will allow us to tackle more complex arbitrary piecewise 1D square-well potentials.

Let us consider a step potential located at the origin,

$$\begin{equation}V\left(x\right)=\begin{cases}_{1},\quad x{< }0\quad \hfill \\ {V}_{2},\quad x{ >}0\quad \hfill \end{cases}\end{equation} \tag{ 14 }$$

To the left and right of the step, depicted in figure 5, the general solution to the Schrödinger equation has the usual form given in equation (1). After imposing continuity of Ψ(x) and Ψ' (x) at x = 0, we obtain two equations that relate the coefficients A, B, C, D:

$$\begin{equation}\begin{gathered}{{\Psi}}_{\mathrm{L}}\left(0\right)={{\Psi}}_{\mathrm{R}}\left(0\right){\Rightarrow}A+B=C+D\\ {{\Psi}}_{\mathrm{L}}^{\prime }\left(0\right)={{\Psi}}_{\mathrm{R}}^{\prime }\left(0\right){\Rightarrow}\mathrm{i}{k}_{1}\left(A-B\right)=\mathrm{i}{k}_{2}\left(C-D\right)\end{gathered}\end{equation} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

These equations can then be compactly written as a matrix equation:

$$\begin{equation}\left(\begin{matrix}\hfill A\hfill \\ \hfill B\hfill \end{matrix}\right)=\frac{1}{2}\left(\begin{matrix}\hfill 1+\frac{{k}_{2}}{{k}_{1}}\hfill & \hfill 1-\frac{{k}_{2}}{{k}_{1}}\hfill \\ \hfill 1-\frac{{k}_{2}}{{k}_{1}}\hfill & \hfill 1+\frac{{k}_{2}}{{k}_{1}}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill C\hfill \\ \hfill D\hfill \end{matrix}\right)={\mathbb{N}}_{12}\left(\begin{matrix}\hfill C\hfill \\ \hfill D\hfill \end{matrix}\right)\end{equation} \tag{ 16 }$$

The matrix $\mathrm{N}$ ₁₂ is formally known as a transfer matrix, which relates the rightward and leftward traveling waves on either side of the step. The transfer matrix is the fundamental solution to the Schrödinger equation for the potential given by equation (14); it neatly stores the matching conditions that relate the four coefficients A, B, C, D, and that must be satisfied by any solution regardless of whether the energy is above or below V₁. While the transfer matrix encodes the step potential information, it is ignorant of additional physical restrictions one might impose, such as boundary conditions at x = ±∞ or normalizability, which is typically achieved by imposing conditions on individual coefficients (i.e., setting them to unity or zero).

If the step potential is located away from the origin, we can shift the step back to the origin by shifting our x coordinate, x → x + a. Under this change of coordinates, the left and right solutions transform as:

$$\begin{equation}\begin{gathered}{{\Psi}}_{\mathrm{L}}\left(x\right)\to {{\Psi}}_{\mathrm{L}}\left(x+a\right)=A\enspace {\text{e}}^{\text{i}{k}_{1}a}\enspace {\text{e}}^{\text{i}{k}_{1}x}+B\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{1}a}\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{1}x}\\ {{\Psi}}_{\mathrm{R}}\left(x\right)\to {{\Psi}}_{\mathrm{R}}\left(x+a\right)=C\enspace {\text{e}}^{\text{i}{k}_{2}a}\enspace {\text{e}}^{\text{i}{k}_{2}x}+D\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{2}a}\enspace {e}^{-i{k}_{2}x}\end{gathered}\end{equation} \tag{ 17 }$$

We see that moving the location of the potential merely adds relative phases (in blue) to the coefficients of the solution. Thus, for a step potential located at x = a (and which we wish to shift back to the origin), the relation between coefficients is:

$$\begin{equation}\left(\begin{matrix}\hfill A\enspace {\text{e}}^{\text{i}{k}_{1}a}\hfill \\ \hfill B\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{1}a}\hfill \end{matrix}\right)={\mathbb{N}}_{12}\left(\begin{matrix}\hfill C\enspace {\text{e}}^{\text{i}{k}_{2}a}\hfill \\ \hfill D\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{2}a}\hfill \end{matrix}\right),\end{equation} \tag{ 18 }$$

We can absorb the phases into the transfer matrix, and instead write:

$$\begin{equation}\left(\begin{matrix}\hfill A\hfill \\ \hfill B\hfill \end{matrix}\right)={\mathbb{N}}_{12}\left(a\right)\left(\begin{matrix}\hfill C\hfill \\ \hfill D\hfill \end{matrix}\right)\end{equation} \tag{ 19 }$$

where the transfer matrix $\mathrm{N}$ ₁₂(a) describes a potential step located at x = a:

$$\begin{equation}{\mathbb{N}}_{12}\left(a\right)=\frac{1}{2}\left(\begin{matrix}\hfill \left(1+\frac{{k}_{2}}{{k}_{1}}\right)\mathrm{exp}\enspace \left[\mathrm{i}\left({k}_{2}-{k}_{1}\right)a\right]\hfill & \hfill \left(1-\frac{{k}_{2}}{{k}_{1}}\right)\mathrm{exp}\enspace \left[-\mathrm{i}\left({k}_{1}+{k}_{2}\right)a\right]\hfill \\ \hfill \left(1-\frac{{k}_{2}}{{k}_{1}}\right)\mathrm{exp}\enspace \left[\mathrm{i}\left({k}_{1}+{k}_{2}\right)a\right]\hfill & \hfill \left(1+\frac{{k}_{2}}{{k}_{1}}\right)\mathrm{exp}\enspace \left[\mathrm{i}\left({k}_{1}-{k}_{2}\right)a\right]\hfill \end{matrix}\right).\end{equation} \tag{ 20 }$$

With the transfer matrix $\mathrm{N}$ ₁₂(a) in hand, it is easy to build solutions to more complicated piecewise constant potentials. Consider a finite well potential that is not necessarily symmetric, noting that this is just two step potentials put together (figure 6). Initially, we will focus on attaining information on the energies of the well, particularly for the bound states, and then solve for the complete wave functions to which they correspond.

Zoom In Zoom Out Reset image size

Outside the well, the solutions to the Schrödinger equation have the form:

$$\begin{equation}\begin{gathered}{{\Psi}}_{\mathrm{L}}\left(x\right)=A\enspace {\text{e}}^{\text{i}{k}_{1}x}+B\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{1}x},\begin{matrix}\hfill \hfill & \hfill x{< }-a\hfill \end{matrix},\\ {{\Psi}}_{\mathrm{R}}\left(x\right)=F\enspace {\text{e}}^{\text{i}{k}_{3}x}+G\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{3}x},\begin{matrix}\hfill \hfill & \hfill x{ >}a\hfill \end{matrix}\cdot \end{gathered}\end{equation} \tag{ 21 }$$

The coefficients A, B, ..., F are related by transfer matrices at each step of the finite square-well,

$$\begin{equation}\left(\begin{matrix}\hfill A\hfill \\ \hfill B\hfill \end{matrix}\right)={\mathbb{N}}_{12}\left(-a\right)\left(\begin{matrix}\hfill C\hfill \\ \hfill D\hfill \end{matrix}\right),\begin{matrix}\hfill \hfill \end{matrix}\begin{matrix}\hfill \hfill \end{matrix}\left(\begin{matrix}\hfill C\hfill \\ \hfill D\hfill \end{matrix}\right)={\mathbb{N}}_{23}\left(a\right)\left(\begin{matrix}\hfill F\hfill \\ \hfill G\hfill \end{matrix}\right)\end{equation} \tag{ 22 }$$

Thus, the coefficients A, B, G, and F can be directly determined without considering the coefficients C and D, by taking the product of two step-potential transfer matrices at the respective location of the steps as [14, 23–26],

$$\begin{equation}\left(\begin{matrix}\hfill A\hfill \\ \hfill B\hfill \end{matrix}\right)={\mathbb{N}}_{123}\left(\begin{matrix}\hfill F\hfill \\ \hfill G\hfill \end{matrix}\right),\begin{matrix}\hfill \hfill & \hfill \qquad {\mathbb{N}}_{123}\equiv {\mathbb{N}}_{12}\left(-a\right){\mathbb{N}}_{23}\left(a\right)\hfill \end{matrix}\end{equation} \tag{ 23 }$$

And just like that, we have obtained the matching conditions for the finite well, where the width of the well is automatically encoded via the position of neighboring steps when they are multiplied together properly [30].

In our new formalism, setting G = 0, the bound states correspond to the zeros of the top-left element of the overall transfer matrix,

$$\begin{equation}{\left({\mathbb{N}}_{123}\right)}_{11}=0\end{equation} \tag{ 24 }$$

T(E) would be given in general by: [28]

$$\begin{equation}T\left(E\right)=\frac{{k}_{3}}{{k}_{1}}{\left\vert \frac{1}{{\left({\mathbb{N}}_{123}\right)}_{11}}\right\vert }^{2}\end{equation} \tag{ 25 }$$

Plotting this and looking for poles at negative energies yields the energy eigenvalues of the bound states of the well. If k₃ = k₁, it can be shown that the left side of equation (24) reduces to √ β from equations (13), and (24) is then equivalent to the case of β = 0.

Moreover, after a bit of algebra, one finds that equation (25) is also equivalent to:

$$\begin{equation}-\mathrm{i}\mathrm{tan}\enspace \left(2a{k}_{2}\right)=\frac{{k}_{2}\left({k}_{1}+{k}_{3}\right)}{{k}_{2}^{2}+{k}_{1}{k}_{3}}.\end{equation} \tag{ 26 }$$

In our notation, this is the general transcendental equation satisfied by the bound state energies of a finite square-well potential, cf reference [31]. Overall, the condition given by equation (24) leads to equation (26) and is the condition that defines the poles in T(E), as defined by equation (25). This mathematical fact reiterates that the poles of T(E) are the energy eigenvalues of the bound states.

3.2.1. Example: constructing a wave function using the transfer matrix method

The elegance of using the transfer matrix method comes from the fact that we do not need to address solutions inside of the well to determine the bound state energies of the well. The power of the transfer matrix becomes even more apparent when one considers that the coefficients in each region can also be determined without directly addressing the solutions inside the well. Below, we will provide a detailed example of using the transfer matrix method to determine the bound state energies of a two-step asymmetric potential well and compute the associated wave functions. This method can then be directly applied to a potential with more than two steps.

Let us consider the uneven potential well shown in figure 6. With the concrete example of section 2 in hand, we are ready to consider a more generalized approach. We will define the dimensionless quantity x* = (x/a). With this convention, the well has length 2, and the potential goes from x* = −1 to x* = 1. Next, we define a dimensionless potential (V*) and energy (E*):

$$\begin{equation}{V}^{{\ast}}=V\left(\frac{2m{a}^{2}}{{\hslash }^{2}}\right)\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad {E}^{{\ast}}=E\left(\frac{2m{a}^{2}}{{\hslash }^{2}}\right)\end{equation} \tag{ 27 }$$

In this convention, let the values of the potential in each region be such that V*₁ = 1, V*₂ = −2, and V*₃ = 0. We treat the problem as before, but making the substitutions V → V*, E → E*, and x → x*. The first step will be to define our initial conditions on the solutions. In this method, we will be using matrix multiplication, and it is more convenient to define both coefficients in region III, and then use the transfer matrix to determine the coefficients in region II and region I.

Let us define G = 0 as before, which allows us to ensure a positive momentum in region III for E greater than any of the potential pieces. We will also define F = 1, which is an arbitrary choice to make the matrix multiplication more straight-forward, and one that can be modified to normalize the wave function at the end. With our boundary conditions set, we will use equation (23) to compute A and B from F and G. Plotting │1/A│² and looking for poles; one finds a single pole at E* = −1.136 56 (figure 7(a)). This determines the numerical values of k₁, k₂, and k₃, and thus, the numerical values of the transfer matrix at the bound state energy. With this in hand, one can use the transfer matrix at each step to find the coefficients in each region. In our case:

$$\begin{equation}\left(\begin{matrix}\hfill C\hfill \\ \hfill D\hfill \end{matrix}\right)={\mathbb{N}}_{23}\left(a\right)\left(\begin{matrix}\hfill F=1\hfill \\ \hfill G=0\hfill \end{matrix}\right){\Rightarrow}\left(\begin{matrix}\hfill A\hfill \\ \hfill B\hfill \end{matrix}\right)={\mathbb{N}}_{12}\left(-a\right)\left(\begin{matrix}\hfill C\hfill \\ \hfill D\hfill \end{matrix}\right)\end{equation} \tag{ 28 }$$

With all the values of the coefficients known, the total wave function is also known and can be normalized and plotted in the usual way (figure 7(b)).

Zoom In Zoom Out Reset image size

Although Gilmore addressed the construction of wave functions using transfer matrices, the method presented here to generate coefficients and compute wave functions echoes the straight-forward process typically outlined in undergraduate textbooks, and without the introduction of, e.g., different sets of solutions to the Schrödinger equation or hyperbolic trigonometric functions [23]. As we will see, by translating this process into an intuitive code, the wave function of almost any finite piecewise square-well potential, including approximated potentials, can be easily and quickly computed and plotted using the same code by merely updating the parameters of the potential.

Consider the case of many step potentials of varying depths connected to each other, or a so-called 'upside-down cityscape' (figure 8(a)). The procedure outlined above can be extended to determine the energy spectrum and wave functions for this potential, where again we will use the dimensionless variables V*, E*, and x*. One simply set the coefficients in the last region to 0 and 1 for the positive and negative momentum solutions, respectively, and then multiply these coefficients by each step consecutively from right to left. It is possible to simply add a new matrix multiplication line for each new step in the two-step example, but this can get messy as the number of steps increases. To keep things compact, a loop can be created in Mathematica^TM where the values of each potential and step location are stored and called on as the loop cycles through the steps from right to left [23–26], From this, A can be computed, T(E*) = │F/A│² can be plotted for E^* < 6, and the energy of the poles can be identified as shown in figure 8(b). Once the bound state energies are known, they can be plugged back into the transfer matrix. The coefficients for the wave function associated with each pole can then be systematically computed, and the wave function can be plotted (figure 8(a)). As was the case for two-step potential well, the solutions inside the well did not need to be directly addressed in this procedure.

Zoom In Zoom Out Reset image size

The power of this approach is further revealed in the ability to discretize curved potentials as a series of steps and then apply the above procedure to quickly calculate a numeric approximation of the energy spectrum and wave functions for the potential. Some examples of known curved potentials would be the harmonic oscillator and the Morse potential. These potentials go to infinity at one or both ends, but it has been shown that they can be well-approximated by truncating them and observing, e.g., the lowest states of the harmonic potential or potentials with a Gaussian shape [14, 23].

We have chosen the Pöschl–Teller potential to further illustrate the benefits of the transmission coefficient approach because it is a finite potential commonly seen in chemistry and classical mechanics whose solutions exist in closed form, thus facilitating easy validation of our method [32]. In its symmetric form, the Pöschl–Teller potential is:

$$\begin{equation}V\left(x\right)=-\frac{{\hslash }^{2}}{2m{a}^{2}}\left[l\left(l+1\right){\text{sech}}^{2}\left(\frac{x}{a}\right)\right],\end{equation} \tag{ 29 }$$

where l is a real parameter that is greater than zero. For integer values of l, the bound state wave functions of this potential are the Legendre functions:

$$\begin{equation}{P}_{l}^{\mu }\left(\mathrm{tanh}\left(\frac{x}{a}\right)\right),\end{equation} \tag{ 30 }$$

where l = 1, 2, 3, ... and μ = 1, 2, ... l − 1, l. Note: the general solution to the Pöschl–Teller potential for generic l and μ is a linear combination of Legendre functions (${P}_{l}^{\mu }\left(\mathrm{tanh}\enspace \left(x/a\right)\right)$ and ${Q}_{l}^{\mu }\left(\mathrm{tanh}\enspace \left(x/a\right)\right)$). When l is an integer and μ takes certain discrete values, the ${P}_{l}^{\mu }$ solutions give normalizable bound states. The energy eigenvalues of the symmetric Pöschl–Teller potential have the form:

$$\begin{equation}E=-\frac{{\hslash }^{2}{\mu }^{2}}{2m{a}^{2}}.\end{equation} \tag{ 31 }$$

Choosing l = 1, the potential was discretized as a series of square steps starting with two steps and increasing to 20 steps. The discretized potentials and the corresponding wave functions are shown in figures 9(a)–(c). The analytic expression for the wave function of this potential is:

$$\begin{equation}{\Psi}\left(x\right)\propto {P}_{1}^{1}\left(\mathrm{tanh}\left(\frac{x}{a}\right)\right)\propto \enspace \text{sech}\left(\frac{x}{a}\right)\end{equation} \tag{ 32 }$$

Zoom In Zoom Out Reset image size

As can be seen from figures 9(a)–(c), the wave functions computed from the discretized potential are in excellent qualitative agreement with the analytic expression [32], which further demonstrates this method's utility. The value of the dimensionless ground state energy (after dividing E₀ by ħ²/2ma²) should be −1, according to equation (31). Figure 9(d) shows the dimensionless ground state energy as a function of the number of steps used to approximate the potential, and after six steps, there is a clear trend toward −1. The fractional uncertainty of the estimated ground state energy compared to the theoretical value is shown in figure 9(e) and is ∼0.02 by as few as twelve steps.

It is well-known that the number of bound states increases with the depth of a potential, and one might expect that anomalous bound states may arise depending on the level of discretization. However, only one bound state is found in all cases, regardless of the number of steps used. It turns out to be true in general that the number of bound states remains accurate regardless of the number of steps used; [23] however, the value for the energy of that state becomes more accurate as the number of steps used increases. The potential goes to zero asymptotically. The width of the plotted area was chosen such that the steps on the outside were at zero and sufficiently followed the actual potential. Analogous to increasing the number of steps, in this example, the energy eigenvalue does not become significantly more accurate by extending the width beyond the range of ±4 nm.

Next, we will increase the depth of the Pöschl–Teller potential by changing the value of l to observe where new bound states occur, as in section 2.3. To explore this, let us discretize the Pöschl–Teller potential using a seven-step potential where the widths of each step are the same, and the depths are changed to roughly follow the exact potential at each l (figure 10(a)). With the placement and potential value of each step set, T(E*) for different values of l can easily be computed and plotted, which is shown for two values of l in figure 10(b). If one plots T(E*) for the exact form of the Pöschl–Teller potential, it can be observed that new bound states form at the top of the well only at integer values of l. As seen in figure 10(b), this result is also observed in the discretized model. Resonance peaks approach and fall into the well only when the step depths correspond to approximately integer values of l.

Zoom In Zoom Out Reset image size

It should be pointed out that T(E* > 0) for the exact Pöschl–Teller potential does not contain oscillating resonances due to the shape of the potential. However, even with just seven steps, the transmission resonances only oscillate by ∼3%, and as more steps are added, the resonance oscillations diminish in amplitude. Lastly, as the l values increase above l = 3, the width of the Pöschl–Teller potential increases significantly along with the depth, and so it becomes necessary to vary the location and depth of the steps to approximate the potential accurately.