T. Scott Palmer II

T. Scott Palmer II

Phys 472

Computational Project

Integrating the One-Dimensional Time-Independent Schrӧdinger Equation for the Infinite “V-well” Potential

In this project, the one-dimensional Time-Independent Schrӧdinger Equation (TISE) is used to calculate the wavefunction as a function of distance, ѱ(x), for the potential V(x)=α|x|, where α is the slope with units eV/Å. This potential is referred to as the “V-well” potential due to its symmetric, linear form. The TISE is written as

, (1)

where m is the mass of the particle, ħ is the reduced Planck constant, and E is the initial energy of the particle in eV. Before integration can begin, we must first dimensionalize the TISE so that the program does not integrate values with physical units. This is a means of avoiding excessive round-off error and catastrophic cancellation in the integration. Also, since the fourth-order Runge-Kutta method (RK4) can be used to integrate second order ODEs after the second order ODE is written as a coupled system of first order ODEs, we must define a variable as the rate of change of ѱ with respect to x. These new constants are

, , (2)

where β is set equal to 1.0 (eVÅ)-1 (Finding Energy Eigenvalues for Symmetric Potentials from Numerical Integration of the Schrӧdinger’s Equation 4). The new variable η allows us to write the TISE as

, (3)

so that we can define a function g(x,η(x)), which the RK4 will use to integrate ѱ, as

. (4)

This is needed because in order to integrate, the RK4 needs to estimate the slope locally so that it can approximate the next value of ѱ using a combination of linear approximation (Euler’s Method) and midpoint rule. For n integrations, then,

(5)

, (6)

where Δx is the size of the step taken between points in the integration. The step size is calculated as the reciprocal of the number of steps taken for a full integration, and k1, k2, k3, k4 are defined as

, (7)

, (8)

, (9)

. (10)

The first and fourth k values integrate using Euler’s method of linear integration, whereas the second and third k values integrate using a midpoint approximation. The weighted average of these k values are then used to integrate ѱ to the next step, x+Δx.

Having written the RK4 using the parameters defined above, the program is then calibrated for a potential that can be solved analytically. The chosen case was the symmetric infinite square well, with β=1.0, E=5.0, η0=0.0, ѱ0=1.0, because given these values the resulting ѱ should be cosine in form and represent the ground state solution to the infinite square well potential (Finding Energy Eigenvalues for Symmetric Potentials from Numerical Integration of the Schrӧdinger’s Equation 5). Using the RK4 with these initial conditions produces a wavefunction which is plotted below in Figure 1.

Figure 1

As can be seen in Fig. 1, ѱ is a cosine function which is zero-valued at the classical turning point of x=1.405 Å. To avoid accumulating excessive round-off error, wavefunctions are integrated from the middle of the infinite square well, x=0, to beyond the classical turning point instead of integrating for the full well width. Since it is impossible to give the potential a truly infinite value beyond the classical turning point, it is only necessary to determine that the wavefunction is effectively zero-valued at the classical turning point, and a cosine before the turning point, to check for ground state symmetry. Determining the degree to which this cosine function is precisely the analytical cosine solution is a means of estimating the relative error compounded by the program during integration. So, to determine the optimal timestep to be used, I defined the machine error, ε, as

, (9)

where the amplitude and wavenumber of the cosine function are determined from the incident energy and the value of the full well width. The resulting error plot is show below in Fig. 2.

Figure 2

As can be seen in Fig. 2, the error peaks periodically, which is a consequence of the periodicity of the cosine function. The height of the peaks increases as the number of steps increases as well, but this is caused from the longer integration time, which compounds more round-off error as the step size decreases. The optimal number of time steps chosen was 1724, since it rested between large spikes of error. The error at this number of time steps was 0.004977, the lowest order of magnitude for error to be compounded.

Having calibrated the program, the absolute value potential V(x)=α|x| was then substituted into Eqn. 1, where α was arbitrarily set at 5.0 eV/Å, as shown in Fig. 3 below. To find the first four wavefunctions for this potential that share the same odd/even parity as the excited state wavefunctions for the infinite square well, the input energies corresponding to these states must first be determined.

Figure 3

Since there exists a full continuum of states between the energies that produce the desired parity, randomly guessing the correct input energy rarely produces the desired wavefunction. To guide the intuitive guesswork, the Variational Principle was applied to find the maximum energy possible in the ground state configuration by assuming the resulting wavefunction is Gaussian in form (Griffiths Problem 7.1(a)). This gives an estimate of

(10)

for the maximum ground state energy possible. This estimate depends only on the slope parameters for V(x) and our normalization constant, which is already of order unity. Using E1 as the input energy produces the wavefunction shown below in Fig. 4.

Figure 4

This ѱ follows the expected parity of the ground state wavefunction with the inflection point at the classical turning point x=1.405 Å, and decays exponentially beyond the classical turning point. Although the wavefuction was assumed to be Gaussian when applying the Variational Principle, the estimation is still valid because the input energy produced the wavefunction with the expected parity. What causes the decay to continue and then diverge is the round-off error introduced by the precision with which E1 can be calculated (Finding Energy Eigenvalues for Symmetric Potentials from Numerical Integration of the Schrӧdinger’s Equation 2). This same effect occurs for solutions that follow the parity for the first three excited states, shown below in Figs. 5-7.

Figure 5

Figure 6

Figure 7

To find the even-parity functions, the initial conditions used were η0=0.0, due to the location of maxima/minima of even-parity functions at x=0Å, and ѱ0=1.0, due to the linearity of the Schrӧdinger equation which allows for any arbitrary finite value to be used as the initial value. For the odd-parity functions, ѱ0 must equal zero since the odd-parity wavefunctions change sign at x=0Å, but because the linearity of the TISE still holds, η0 can still be set to 1.0 (Finding Energy Eigenvalues for Symmetric Potentials from Numerical Integration of the Schrӧdinger’s Equation 2). The trial energies E2=6.840709eV, E3=9.494709eV, E4=11.949709eV were determined by iterating guesses of energy values higher than E1until the expected parity was produced. This method was able to find precise input energies whose solutions followed the expected parity between the classical turning points, but it did take many trials to guess more appropriate input energies. If a trial energy was too high, the divergence occurred at smaller x values—even before the classical turning point in some cases. If a trial energy was too low, the wavefunction would turn away from the limiting value of zero, and it would diverge in the opposite direction. By iterating guesses and using higher and higher orders of decimal precision for the energies, solutions were generated that were more stable for longer integrations in x and demonstrated the exponential decay beyond the classical turning points more clearly than before. Although it is not mathematically rigorous, this method can be considered useful because the continuum of energy states that exist between states corresponding to a given parity means that there are many similar solutions in a small range of energies. As expected, the classical turning points move out to farther x as the input energy increases. This is an artifact of the symmetry and linearly increasing behavior of the “V-well” potential. To improve the numerical calculation in the future, a method of exactly determining the excited state energies would need to be found to avoid the divergences evident in Figs. 5-7.

By integrating the Time-Independent Schrӧdinger Equation using the fourth-order Runge-Kutta technique, the ground state and first three excited state wavefunctions were found for the “V-well” potential and trial energies of E1=2.99471eV, E2=6.840709eV, E3=9.494709eV, E4=11.949709eV.