Numerical Solution to Schrodinger Equation in a 1D Periodic Potential

Summary:

The following article serves as a discussion of a program that I designed to serve as an educational tool for students of introductory solid state physics.  This program calculates the dispersion relation for a particle in periodic 1D potential and presents related information, such as band gap energies and movies of the time-evolved stationary states directly above and below the first band gap.  I begin by introducing the reader to the underlying physical concepts, such as energy gaps, dispersion curves, Bloch’s theorem, and the nearly-free electron model.  I then proceed with a discussion of the program and how it can be used to demonstrate a wide array of phenomena related to to these theoretical topics.  I conclude with a description of how the program operates, its limitations, and suggestions for how this work can be extended in the future.

Introduction:

In many solids, such as crystals and semiconductors, atoms are arranged in an organized and periodic manner. This periodicity in the organization of the atoms leads to an electric potential that is also periodic in space. Describing the behavior of particles, such as electrons, in such potentials is a primary topic in any introductory solid state physics course, and remains an important area of research in modern electronics.

An interesting phenomenon that occurs in periodic potentials, and a topic which is at the heart semiconductor device operation, is the presence of energy bands. In many solids, there are intervals of allowed energies which particles experiencing that potential can take. That is, if one were to measure the energy of an electron in a periodic potential, they might measure values anywhere between two energies, say, a and b, but nowhere between two other energies, say c and d. The allowed energy intervals are known as energy bands and the energy difference between two bands is known as a band gap [1].

In introductory solid state texts, such as “Solid State Physics,” by Kittel, the origin of the band gap is often presented in the context of the nearly-free electron model [1]. In the nearly-free electron model, electrons are described as free particles with an effective mass that depends on the electron’s wavenumber. The energy gap is then explained conceptually as being due to the existence of standing waves that form when the wavelength of the electron is some harmonic of the periodicity of the potential. Two distinct standing waves form, one in which the electron’s wavefunction is concentrated around regions of high potential energy, and one in which the electon’s wavefunction is concentrated around regions of low potential energy. The energy difference between these two standing wave states is taken as the bang gap energy.

While presenting band gaps in this way avoids grappling with the full complexity of quantum mechanics, band gaps are an essentially quantum mechanical phenomenon. By the time many students begin studying solid-state physics, they have already studied at least one of the archetypal quantum mechanics examples, such as the infinite square well or the hydrogen atom. Furthermore, they have learned that only certain energies have associated normalizable solutions to the Schrodinger equation, and that these are the measurable energies of the system. The only difference in infinite periodic potentials, is that instead of discrete energy levels, there are intervals of allowed energies. The goal of this project is to help students bridge the gap between the conceptual nearly-free electron model and the results of quantum mechanics.

Bloch’s Theorem:

One of the most important theorems involving solutions to the Schrodinger equation in a periodic potential is Bloch’s theorem. This states that the normalizable solutions to the time independent Schrodinger equation in a periodic potential have the form \psi_k=u_k(\vec {r})e^{i \vec{k} \cdot \vec{r}}, where \vec{r} is the position vector, \vec{k} is the wave vector, u_k(\vec {r})=u_k(\vec {r}+\vec{T}), and \vec{T} is a lattice vector (the position vector between two lattice sites) [1]. In one dimension, this reduces to the criterion that \psi_k=u_k(x)e^{ikx}, where u_k(x)=u_k(x+na), a is the periodicity of the potential, and n is an integer. It follows from Bloch’s theorem that if \psi_k is a solution to the Schrodinger equation, then \psi_{k+G} is also a solution, provided that G=\frac{m \pi}{a} for some integer m [1]. It is common to plot the allowed energies against k for \frac{- \pi}{a} < k < \frac{\pi}{a} (known as the first Brillouin zone), as this gives a complete account of the allowed energies of the system [1]. The energy of the allowed wavefunctions as a function of k is known as the dispersion relation or dispersion curve [1].

Project and Results:

I created a program that calculates the dispersion relation for three 1-dimensional periodic potentials, a square, triangle, and sinusoidal wave (figure 1).

Fig 1: Periodic potentials which are available which are available to be solved by simulation. This is understood to represent a single cell of the potential functions. The true potentials infinitely many cells sewn together at the edges.

Fig. 1: Periodic potentials available available in the program. This is understood to represent a single cell of the potential functions. The true potentials repeat infinitely many times in both the positive and negative spacial directions.

This program is intended as an educational tool for students of solid-state physics who are beginning to study the concept of bang gaps. It calculates the dispersion curves by finding the first several allowed energy eigenstates (of the potential of interest) for a discrete selection of k values in the first Brillouin zone. It then plots the energies of these states against the wave number, k, which gives the dispersion curve (figure 2). The first several band gaps are evident as the empty regions between dispersion curves.

tr

tr1

Fig. 2: Dispersion relations calculated by the program. The top image corresponds to the triangle-wave potential, the middle image corresponds to the sinusoidal potential, and the bottom image corresponds to the square-wave potential. Regions between the colored horizontal lines represent band gaps.

The program also produces a movie of the time-evolved eigenstates on either side of the first band gap (figures 3 and 4). In these movies, it is worth noting that the state just above the first energy band is concentrated around high-energy regions of the potential and the state just below the first band is grouped around low-energy regions of the potential. One can also see that the periodicity of the allowed wavefunctions on either side of the gap is indeed a harmonic of the periodicity of the potential (actually a sub-harmonic). These results demonstrate that the qualitative explanation of band gaps provided by the nearly-free electron model is in close agreement with the results of quantum mechanics.

www-gifcreator-me_jyqlj

Fig. 3: Animation of the real and imaginary parts of the wavefunction just below the first band gap (click image to see animation)

Fig. 4: Animation of the real and imaginary parts of the wavefunction just above the first band gap (click image to see animation)

Along with having several types of potentials to choose from, the height and periodicity of the potential are left as variable inputs for students to change and explore. By changing these values, students can begin to build and intuition for how the structure of the dispersion curves depend on the size and shape of the potential. For instance, by decreasing the height of a potential of fixed periodicity, students are able to watch the structure transition from being an insulator (large band gap), to a semiconductor (small band gap), to a metal (little or no band gap) [1].

Fig. 4 Dispersion relation for decreasing magnitude of the potential barriers. The decrease in energy gaps between dispersion curves represent the transition from insulator to conductor.

Fig. 5: Dispersion relations for decreasing magnitude of the potential barriers. The decrease in energy gaps between dispersion curves demonstrate the transition from insulator to conductor.

The program can also be used as a tool for students to build a better understanding of Bloch’s theorem. For instance, by instructing the program to calculate the allowed energies for values of k beyond the first Brillouin zone, we see that the dispersion relation is periodic in k (figure 6). This corresponds to the result from Bloch’s theorem that if \psi_k is a solution to the Schrodinger equation in a periodic potential, then \psi_{k+G} is also a solution [1].

Fig. 5 Dispersion relation, which is periodic in k

Fig. 6: Representative dispersion relation that extends beyond the first Brillouin zone.  Note the periodicity in k, as predicted by Bloch’s theorem.

Verification of Results –  The Kronnig-Penney Model:

In order to demonstrate that this program produces solutions that are in agreement with theory, I considered the square-wave potential. This potential corresponds to the Kronnig-Penney model, which has been studied extensively. There is an analytic solution for the allowed energies of this potential in the form of roots to a transcendental equation [1]. Much like the methods employed by this program, this expression is found using Bloch’s theorem to determine appropriate wavefunction boundary conditions.

This equation for the allowed energies has been used to find the analytic dispersion curve for the Kronnig-Penney model. This is plotted along with the dispersion curve determined by the numerical methods employed in this program (figure 7). Indeed, the two dispersion relations are in excellent agreement, which demonstrates that this program does reliably calculate the allowed energies for a 1D periodic potential.

Figure 1: The top image shows the dispersion curve for the square-wave potential calculated by numerical methods. The bottom image shows the dispersion curve calculated using analytic results obtained from Bloch's theorem.

Fig. 7: The top image shows the dispersion curve for a square-wave potential calculated by numerical methods. The bottom image shows the dispersion curve calculated using analytic results obtained from Bloch’s theorem.

A more in depth discussion of the Kronnig-Penney model can be found in “Solid State Physics,” by Kittel [1].

Program Operation:

In order to determine the allowed energies of a particle in a periodic potential, I employed the numerical methods discussed in “Numerical Solution of the 1D Schrodinger Equation: Bloch Wavefunctions,” by Dr. Constantino Diaz [2]. The theory and operation of the procedure outlined in this paper is structured in the following way. The time independent Schrodinger equation is a second order differential equation in x, so there exists two linearly independent solutions for each value of E, which we will denote as \psi_C(k,E,x) and \psi_S(k,E,x) [2,3]. In other words, any wavefunction with wavenumber k and energy E, is expressible as a linear combination of \psi_C(k,E,x) and \psi_S(k,E,x). That is, \psi_k(E,x)=A\psi_C(k,E,x)+B\psi_S(k,E,x) for some complex numbers A and B [2,3].

\psi_C(k,E,x) and  \psi_S(k,E,x) are determined for a chosen value of k by using the Schrodinger equation to calculate the wavefunction in discrete steps [2]. This is done by expressing the Schrodinger equation in finite difference form:

    \[\frac{\hbar^2}{2m} \frac{d^2}{dx^2} \psi+V\psi=E\psi \approx \frac{\hbar^2}{2m} \frac{\psi_{n+1}+\psi_{n-1}-2\psi_{n}}{\Delta x^2}+V_n \psi_n=E\si_n,\]

where V(x) is the potential function, m is the particle mass, \Delta x is the spacial step size, and n denotes the n^{th} spacial spacial step [3]. Rearranging terms gives:

    \[ \psi_{n+1}=\frac{2m(\Delta x^2)}{\hbar^2}(E-V_n)\psi_n-\psi_{n-1}+2\psi_n,\]

which can be used to calculate the wavefunction for some guess value of E [3].  Since the potential functions are symmetric, initial conditions of \psi_C(\frac{a}{2})=1, \frac{d\psi_C}{dx}(\frac{a}{2})=0 and \psi_S(\frac{a}{2})=0, \frac{d\psi_S}{dx}(\frac{a}{2})=1 can be chosen without loss of generality [2,3].

According to Bloch’s theorem, \psi(a)= u_k(a)e^{ika}=u_k(0)e^{ika}=\psi(0)e^{ika}, and \psi^{'}(a)= u_k^{'}(a)e^{ika}+iku_k(a)e^{ika}=u_k^{'}(0)e^{ika}+iku_k(0)e^{ika}=\psi^{'}(0)e^{ika} [1,2]. These boundary conditions can be expressed in terms of the coefficients A and B as A\Delta C+B\Delta S=0 and A\Delta C^{'}+B\Delta S^{'}=0, where \Delta C=\psi_C(a)e^{\frac{ika}{2}} -\psi_C(0)e^{\frac{-ika}{2}, \Delta S=\psi_S(a)e^{\frac{ika}{2}} -\psi_S(0)e^{\frac{-ika}{2}, \Delta C^{'}=\psi_C^{'}(a)e^{\frac{ika}{2}} -\psi_C^{'}(0)e^{\frac{-ika}{2}, and \Delta S^{'}=\psi_S^{'}(a)e^{\frac{ika}{2}} -\psi_S^{'}(0)e^{\frac{-ika}{2} [2]. For this system of equations to be satisfied, the determinant of the coefficients must be equal to zero. That is, the allowed energies for a given value of k satisfy the equation, G(E)=\Delta C \Delta S^{'} - \Delta S \Delta C^{'} =0. If G(E)=0, then E is taken as an allowed energy [2]. If not, then a new value E+\Delta E is treated as the energy and new values of \Delta C, \Delta S, \Delta C^{'}, and \Delta S^{'} are computed. This process is repeated many times so that several allowed energies are found for the chosen value of k. The whole procedure is then repeated for evenly spaced values of k in the first Brillouin zone.

Program Limitations:

The primary limitation of this program is its sensitivity to choices of step size.  If \Delta E is chosen to be too small, an insufficient number of roots are found for G(E) to create a viable dispersion curve.  If the energy step size is chosen too large, then the program is unable to locate the zeros of G(E), particularly at the ends of the first Brillouin zone.  In both of these instances, the program stops running and produces plot of G(E) for the last calculated value of k. If this graph has less than four roots, then the step size was chosen to be too small.  If the graph has many roots, but the program still failed, then the step size was chosen to be too big. The program seems to be most effective if \Delta E is chosen to be roughly one percent of the maximum of the potential functions.

The program also fails to produce reliable results if the spacial step size is chosen to be too large.  As long as the spacial step size is no more than five percent of the width of the potential barrier and 0.5 percent of the periodicity of the potential, the program seems to operate normally.  The program also produces poor results if the step size in k is chosen to be too large, but this only decreases the resolution of the dispersion curve.  An overly large choice in the step size of k will not result in a failure of the program unless it is chosen to be comparable to the size of the first Brillouin zone. As long as operators are aware of these problems, the program performs well.  Still, I recommend that only small changes in the parameters of the potential be made each time the program is run so that the user is able to adjust the step sizes gradually.

Future Work:

A good continuation of this project would be to alter the program to accommodate for potentials that are not symmetric about the point x=\frac{a}{2}. Suggestions on how to implement this modification can be found in “Numerical Solution of the 1D Schrodinger Equation: Bloch Wavefunctions,” by Dr. Constantino Diaz [2]. An even more interesting progression of this work would be to write a program that can solve the Schrodinger equation in a two or three-dimensional potential lattice. The particular method employed in this program is not generalizable to higher dimensions, but other methods that also apply Bloch’s theorem to provide appropriate boundary conditions can be applied to periodic potentials of higher dimension [2]. In the case of higher dimensional potentials, the dispersion relation is anisotropic with respect to spatial orientation. Qualitatively, this is due to the fact that the spacing between potential barriers is different along distinct spatial directions [1]. A program that could allow students to explore this directional dependence of the dispersion relation and plot the associated wavefunctions could be a helpful tool to solidify the effect of spatial anisotropy on wavefunctions in crystals.

Conclusions:

I created a program that calculates the dispersion relation for a particle in a 1D periodic potential, which can be used as a educational tool for students of solid state physics.  It uses the discrete version of the Schrodinger equation along with boundary conditions provided by Bloch’s theorem in order to find the first several allowed energies for values of k in the first Brillouin zone.  Comparison of the results of this program to the analytic solution for the dispersion relation of Kronnig-Penney demonstrates that this program provides valid results. The program demonstrates several key properties of wavefunctions in a periodic potential, including energy bands, the distribution of wavefunctions on either side of the band gap around high and low regions of the potential function, the dependency of energy gaps on the height of potential function, and the periodicity of the dispersion curve predicted by Bloch’s theorem. A good continuation of this work would be to incorporate non-symmetric potential functions or to construct a procedure for calculating the dispersion relation in higher dimensional potentials.

References:

[1] Kittel, Charles. Solid State Physics. 8th ed. New Caledonia: John Wiley & Sons, 2005. Print.

[2] Uteras Diaz, Constantino A. “Numerical Solution of the 1D Schrodinger Equation: Bloch Wavefunctions.” Numerical Solution of the 1D Schrodinger Equation: Bloch Wavefunctions. ARXIV, 23 Mar. 2005. Web. 10 Dec. 2016. <http://cds.cern.ch/record/828907/files/0503171.pdf>.

[3] Giordano, Nicholas J. Computational Physics. 2nd ed. Upper Saddle River, NJ: Prentice Hall, 2006. Print.

 

Share