A particle in a one-dimensional finite square well potential¶

This notebook evaluates the allowed energies and wavefunctions for bound states of a particle confined to a finite square well potential.

Introduction¶

Consider a particle of mass, $m$, that is moving in the $x$-direction and experiences a finite square potential of width, $L$, and height, $V_0$. We will assume that the energy of this particle, $E$, is less than $V_0$. Graphically, this potential has the form:

In [1]:
from IPython.display import Image
Image(url = "https://raw.githubusercontent.com/edeprince3/super_coding_fun_time/main/square_well_potential/square_well_potential.png")
Out[1]:

Mathematically, we define this potential in a piecewise way as

Region I:¶

$V(x) = V_0$, $x < 0$

Region II:¶

$V(x) = 0$, $0 \le x \le L$

Region III:¶

$V(x) = V_0$, $x > L$

Wavefunctions¶

Wavefunctions that satisfy the time-independent Schrodinger equation will also have a piecewise form, and we can determine the wavefunction by solving the Schrodinger equation in each region. We have

Region I:¶

$(\frac{-\hbar^2}{2m}\frac{d^2}{dx^2} + V_0) \psi_I(x) = E \psi_I(x)$

which has solutions of the form

$\psi_I(x) = C~{\rm exp}(\alpha x) + D~{\rm exp}(-\alpha x)$, $x < 0$

Region II:¶

$\frac{-\hbar^2}{2m}\frac{d^2}{dx^2} \psi_{II}(x) = E \psi_{II}(x)$

which has solutions of the form

$\psi_{II}(x) = A~{\rm cos}(\beta x) + B~{\rm sin}(\beta x)$, $0 \le x \le L$

and

Region III:¶

$(\frac{-\hbar^2}{2m}\frac{d^2}{dx^2} + V_0) \psi_{III}(x) = E \psi_{III}(x)$

which has solutions of the form

$\psi_{III}(x) = F~{\rm exp}(\alpha x) + G~{\rm exp}(-\alpha x)$, $x > l$

Above, $\alpha$ and $\beta$ are real numbers, defined by

$\alpha = ( 2 m [V_0 - E] / \hbar^2)^{1/2}$

$\beta = (2 m E / \hbar^2)^{1/2}$

Boundary conditions¶

If we consider that the wave function should be finite in the limit that $x$ tends to $\pm \infty$, then we immediately find that the coefficients $D$ and $F$ must be zero.

As for the other unknown coefficients, we can determine these through the application of various boundary conditions:

  1. The wavefunction should be continuous between regions I and II (at $x=0$). This condition leads us to

    $C = A$

  2. The derivative of the wavefunction should be continuous between regions I and II (at $x=0$). This condition leads us to

    $B = [(V_0 - E)^{1/2} / (E)^{1/2}] A$

  3. The wavefunction should be continuous between regions II and III (at $x=L$). This condition leads us to

    $G = A ~[ {\rm cos}(\beta L) + \alpha~ / \beta ~{\rm sin}(\beta L) ] ~{\rm exp}(\alpha L)$

  4. The derivative of the wavefunction should be continuous between regions II and III (at $x=L$). This condition leads us to a trancendental equation for the energy:

    ${\rm tan}[(2mE/\hbar^2)^{1/2} L] = 2 (V_0-E)(E)^{1/2}/(2 E-V_0)$

Lastly, the coefficent A can be determined by normalization.

Procedure for visualizing the allowable wavefunctions¶

In order to determine the allowable energies and corresponding wavefunctions for the particle in a finite square well potential, we will follow the following steps:

  1. We should specify parameters $m$, $L$, and $V_0$ that define our problem.

  2. We should plot the trancendental equation for the energy to (a) determine how many bound states we have and (b) to obtain reasonable guesses for these energies.

  3. We should numerically solve the trancendental equation using some functionality in the scipy package.

  4. Once we have the allowable energies, we can evaluate the corresponding wavefunction parameters defined above.

  5. Given the wave function parameters, we can visualize the wavefunction or its square modulus.

In [2]:
# step 1: specify the parameters, in atomic units (hbar = 1):

m = 1.0
V0 = 10.0
L = 1.0
In [3]:
# step 2: plot the trancendental equation for the energy

import numpy as np
import matplotlib.pyplot as plt

# energies should be less than V0 (we're looking at bound states)
dE = 0.01
E = np.arange(0, V0, dE)

# LHS = tan[sqrt(2mE/hbar^2)*L]
# RHS = 2*sqrt(V0-E)sqrt(E)/(2*E-V0)
lhs = []
rhs = []

for i in range(len(E)):
    lhs.append( np.tan(np.sqrt(2.0 * m * E[i]) * L) )
    rhs.append( 2.0 * np.sqrt(V0 - E[i]) * np.sqrt(E[i]) / ( 2.0 * E[i] - V0) )

plt.plot(E, lhs) 
plt.plot(E, rhs)
plt.ylim(-10,10)
ax = plt.gca()

plt.show()

/Users/deprince/miniconda3/envs/pq/lib/python3.7/site-packages/ipykernel_launcher.py:17: RuntimeWarning: divide by zero encountered in double_scalars
  app.launch_new_instance()

From above, we can see that there are two allowable energies (points at which the blue and orange curves cross, not counting the divergence in the blue curve). These energies look like they are roughly 2.5 and 8.5 (atomic units). We now need to numerically solve the trancendental equation. We're going to use the function:

optimize.fsolve

from the scipy package for this purpose. We will have to pass some objective function $f(E) = 0$ for which this function will find the optimal $E$ value(s).

In [4]:
def allowable_energies(E, *data):
    
    """
    
    function defining allowable energies:
    
    tan[sqrt(2mE/hbar^2)*L] - 2*sqrt(V0-E)sqrt(E)/(2*E-V0) = 0
    
    :param E: the energy
    :param data: m, V0, and L
    
    :return value: the function value
    
    """
    
    m, V0, L = data
    
    return np.tan(np.sqrt(2.0 * m * E) * L) - 2.0 * np.sqrt(V0 - E) * np.sqrt(E) / ( 2.0 * E - V0)

Now, we're ready to solve for the allowable energies.

In [5]:
# step 3: solve for allowable energies, with initial guess informed by plot above

from scipy import optimize
E = optimize.fsolve(func = allowable_energies, x0 = [2.5, 8.0], args = (m, V0, L) )

print(E)

[2.29499075 8.13714776]

Given the allowable energies, we can calculate the corresponding wavefunction parameters. Since all of the parameters can be expressed in terms of $A$, we start with the choice $A = 1$, but we'll rescale all of the coefficients such that the wavefunction is normalized.

In [6]:
# step 4. evaluate wavefunction parameters for each of the allowable energy levels

alpha = np.sqrt( 2.0 * m * (V0 - E) )
beta = np.sqrt( 2.0 * m * E )
    
A = 1 # we'll normalize, but start with A = 1
C = A
B = np.sqrt(V0 - E) / np.sqrt(E) * A
G = A * ( np.cos(beta * L) + alpha / beta * np.sin(beta * L) ) * np.exp(alpha * L)

Now, we can scale the coefficients ($A$, $B$, $C$, and $G$) such that the wavefunction is normalized, i.e., such that

$\int_{-\infty}^{\infty} |\psi(x)|^2 dx = 1$

To normalize the wave function, we should scale $\psi(x)$ such that

$\psi(x) \to \psi(x) / \left ( \int_{-\infty}^{\infty} |\psi(x)|^2 dx \right )^{1/2}$

and will need to evaluate the integral in a piecewise manner, as

$\int_{-\infty}^{\infty} |\psi(x)|^2 dx = \int_{-\infty}^{0} |\psi_I(x)|^2 dx + \int_{0}^{L} |\psi(x)_{II}|^2 dx + \int_{L}^{\infty} |\psi(x)_{III}|^2 dx$

Each of these integrals can be evaluated analytically as

$\int_{-\infty}^{0} |\psi_I(x)|^2 dx = |A|^2 / ( 2 \alpha)$

$\int_{0}^{L} |\psi(x)_{II}|^2 dx = (2 \beta L (A^2 + B^2) + (A - B) (A + B) {\rm sin}(2 \beta L) - 2 A B {\rm cos}(2 \beta L) + 2 A B)/(4 \beta)$

and

$\int_{L}^{\infty} |\psi(x)_{III}|^2 dx = {\rm exp}(-2 \alpha L) / (2 \alpha)$

In [7]:
# step 4 continued. normalize the wave function
    
nrm  = A*A / (2.0 * alpha) 
nrm += (2 * beta * L * (A**2 + B**2) + (A - B)*(A + B)*np.sin(2.0 * beta * L) - 2.0 * A * B * np.cos(2 * beta * L) + 2.0 * A * B)/(4.0 * beta)
nrm += G*G * np.exp(-2.0 * alpha * L) / (2.0 * alpha)
    
nrm = 1.0 / np.sqrt(nrm)
    
A *= nrm
B *= nrm
C *= nrm
G *= nrm
In [8]:
# step 5. plot the wavefunction (and the shape of the potential)

dx = 0.01
left = np.arange(-1, 0, dx)
inside = np.arange(0, L + dx, dx)
right = np.arange(L+dx, 2 + dx, dx)

V = []
psi = []
psi2 = []
x = []

# psi(x) = C exp(alpha x), x < 0
for i in range(len(left)):
    V.append(1.0)
    x.append(left[i])
    val = C * np.exp(alpha * left[i])
    psi.append(val)
    psi2.append(val*val)

# psi(x) = A cos(beta x) + B sin(beta x), 0 < x < L
for i in range(len(inside)):
    V.append(0.0)
    x.append(inside[i])
    val = A * np.cos(beta * inside[i]) + B * np.sin(beta * inside[i]) 
    psi.append(val)
    psi2.append(val*val)
    
# psi(x) = G exp(-alpha x), x > L
for i in range(len(right)):
    V.append(1.0)
    x.append(right[i])
    val = G * np.exp(-alpha * right[i])
    psi.append(val)
    psi2.append(val*val)

plt.plot(x, V)
plt.plot(x, psi) 
ax = plt.gca()

plt.show()

Now what?¶

Now, you should re-run this notebook with different values for the mass ($m$), the width of the well ($L$), and the height of the potential outside of the well ($V_0$). Ask yourself the following questions:

  1. Does increasing $V_0$ increase or decrease the number of bound states?

  2. Does increasing $L$ increase or decrease the number of bound states?

  3. Does increasing $m$ increase or decrease the number of bound states?

  4. For a given $m$, $L$, and $V_0$, which bound states have the highest / lowest tunneling probability (the probability of finding the particle outside of the well)? Can you rationalize this result?

  5. Can the number of bound states ever be zero, for certain combinations of $m$, $L$, and $V_0$?