Free Particle (Python Notebook)

Describing the Free Particle

We start by describing a free particle: a particle that is not under the influence of a potential. As any other particle, the state of a Free Particle is described with a wavefunction \(\psi(x)\). In order to learn about this particle (i.e., measure its properties) we must construct Hermitian Operators. For example, what is the momentum operator \(\hat P\)?

The momentum operator most obey the following property (eigenfunction/eigenvalue equation):

\[\hat P \psi_k(x) = \lambda  \psi_k(x) \label{1}\]

where \(\lambda\) is an eigenvalue of a Hermitian operator and therefore it is a real number.

In the \(x\) representation, using the momentum operator as \(\hat P =-i\hbar \frac{\partial }{\partial x}\), we can solve Equation \ref{1} by proposing a function to represent \( \psi_k(x) \) as \(\psi_k(x) = c\ e^{ikx}\), where \(k\) is a real number.

Let's see if this functional form for the wavefunction works (i.e., solved Equation \ref{1}):

\[\hat P \psi_k(x) =p \psi_k(x)\]

\[-i\hbar \frac{\partial {c\ e^{ikx}}}{\partial x} =-i\hbar\ c\ ik\ e^{ikx} \]

\[\hbar k\ c\ e^{ikx} = \hbar k\ \psi_k(x) \label{2}\]

with \(p=\hbar k\). Although \(\psi_k(x)\) can not be normalized, since it is not a square-integral function (i.e., the integral of the square of the absolute value is not finite). However, the constant \(c\) is chosen in such a way (see note below) that the function becomes:

\[\large{\psi_k(x) = \frac{1}{\sqrt{2\pi}} e^{ikx}}\]

The Normalization Factor and Dirac's Delta Function

The normalization of a bra-ket is given by

\[\left<\psi_{k'}(x)|\psi_k(x)\right> = \left<k'|k \right> = \int^{+\infty}_{-\infty} \psi_{k'}^*(x)\psi_k(x) dx = \int^{+\infty}_{-\infty} c^* e^{-ik'x} ce^{ikx} dx = |c|^2 \int^{+\infty}_{-\infty} e^{i(k-k')x} dx = |c|^2 2\pi\ \delta(k-k')\]

Where we used the definition of the Dirac's delta function:

\[\delta(x) = \frac{1}{2\pi} \int^{+\infty}_{-\infty} e^{ivx} dv \nonumber\]

which obeys the following properties:

\[\delta(x-a)=0\ if x \neq a \nonumber\]

\[\int^\infty_{-\infty} f(x)\ \delta(x-a) \ dx =f(a) \nonumber\]

The normalization constant \(c\) is chosen as \(c = \frac{1}{\sqrt{2\pi}}\) in such a way that:

\[\left<\psi_{k'}(x)|\psi_k(x)\right> = \delta(k-k') \nonumber\]

Once we know the normalized eigenfunctions, we can ask some questions:

Questions

• Q1: What is the momentum of a free particle?
• Q2: What are the units of momentum?
• Q3: Is the momentum of a Free Particle quantized?

In SI units, \(k\) has units of \(\frac{1}{m}\) and \(\hbar\) has units of \(J.s =\frac{kg.m^{2}}{s}\), therefore from Equation \ref{2} we can see that \(p=\hbar k\) and the units are \(\frac{J.s}{m} = \frac{Kg.m}{s}\).

To visualize the wavefunction \(\psi_k(x)\) we use Euler relationship

\[e^{i\theta}=\cos(\theta)+i\sin(\theta)\]

and graph separately the Real and Imaginary contributions to the wavefunction (even though the wavefunction is complex, the \(\hat P\) is a Hermitian operator, thus the eigenvalues are measurable and real numbers).

On the other hand, the Probability Density (\(\psi^*_k(x) \psi_k(x)\)) is a constant, since:

\[ \psi^*_k(x) \psi_k(x) = \frac{1}{\sqrt{2 \pi}}e^{-ikx}\frac{1}{\sqrt{2 \pi}}e^{ikx} =\frac{1}{2\pi}\]

You can use the cell below to enter different values for the momentum wavevector \(k\) and the length of space \(x\)

Interactive Element

import matplotlib.pyplot as plt
import numpy as np

# Defining Psi and PsiC (complex conjugate) functions
def psi(x,k): return (1.0/np.sqrt(2.0*np.pi))*(np.cos(k*x)+np.sin(k*x)*1j)
def psiC(x,k): return (1.0/np.sqrt(2.0*np.pi))*(np.cos(k*x)-np.sin(k*x)*1j)

# Reading the input variables from the user
k = float(input("Enter the value for k (in Angstroms-1) = "))
xmax = float(input("Enter the maximum value for x (in Angstroms) = "))

# Generating the wavefunction graph
plt.rcParams.update({'font.size': 18, 'font.family': 'STIXGeneral', 'mathtext.fontset': 'stix'})
x = np.linspace(-xmax, xmax, 900)
fig, ax = plt.subplots()
lim1=1/np.sqrt(2*np.pi) # Maximum value of the wavefunction
ax.axis([-xmax,xmax,-lim1*1.1,lim1*1.1]) # Defining the limits to be plot in the graph
ax.plot(x, psi(x,k).real, label="Real", color="blue") # Ploting the real part of the wavefunction
ax.plot(x, psi(x,k).imag, label="Imag.", color="red") # Ploting the real part of the wavefunction
ax.hlines(0.0, -xmax, xmax, linewidth=1.8, linestyle='--', color="black") # Adding a horizontal line at 0
# Now we define labels, legend, etc
ax.legend(loc=2);
ax.set_xlabel(r'$x$ (Angstroms)')
ax.set_ylabel(r'$\psi_k(x)$')
str1=r"$k$ = "+str(k)+r" A$^{-1}$"
plt.title('Wavefunction for '+str1)
plt.legend(bbox_to_anchor=(1.1, 1), loc=2, borderaxespad=0.0)

# Generating the probability density graph
fig, ax = plt.subplots()
ax.axis([-xmax,xmax,0.0,lim1*lim1*1.4]) # Defining the limits to be plot in the graph
ax.plot(x, psi(x,k)*psiC(x,k),label="Probability Density", color="green") # Plotting the probability density
# Now we define labels, legend, etc
ax.legend(loc=2);
ax.set_xlabel(r'$x$ (Angstroms)')
ax.set_ylabel(r'$\left|\psi_k(x)\right|^2$')
str1=r"$k$ = "+str(k)+r" A$^{-1}$"
plt.title('Probability Density for '+str1)
plt.legend(bbox_to_anchor=(1.1, 1), loc=2, borderaxespad=0.0)

# Show the plots on the screen once the code reaches this point
plt.show()

What do we learn from this graphical representaton?: The first plot shows that that the wavefunction is completely delocalized in the x coordinate. However, wherever we look, the value of the probability density (\(|\psi^*_k(x) \psi_k(x) |\)) is the same.

An important question to ask is: what is the probability to find a free particle within a \(\Delta x\) region in space?

This is crucial for measurements, since we would never measure ALL space at the same time, but instead, our instrument will measure a section of \(x\) space

\[\Psi^*_k(x) \Psi_k(x) \Delta x = \left|\Psi_k(x)\right|^2 \Delta x = \left(\frac{1}{\sqrt{2\pi}} e^{ikx}\right)^*\frac{1}{\sqrt{2\pi}} e^{ikx} \Delta x = \frac{1}{2\pi} \Delta x \]

The probability of finding it within a \(\Delta x\), at any place in space, is the same.

Thus when the particle has a well-defined momentum (\(\hbar k\)), the probability of finding it anywhere in space is the same.

Localization of the particle can be obtained by considering not just a single value of momentum (\(\hbar k\)), but a range of values.

Consider the superposition of waves with a range of momenta given by

\[p=p_0 \pm \Delta p \nonumber \]

or

\[\hbar k = \hbar(k_0 \pm \Delta k) \nonumber\]

Since the momentum varies continuously, we use an integration for the superposition of waves, all equally weighted:

\[\Psi_{\Delta k}(x) =N\int^{k_o+\Delta k}_{k_o-\Delta k}\ e^{ikx}dk \nonumber\]

\[\Psi_{\Delta k}(x) = \frac{2N\sin(\Delta kx)}{x} e^{ik_ox} \nonumber\]

where \(N\) is a normalization constant and is equal to

\[N=\frac{1}{2\sqrt{\pi\Delta k}} \nonumber\]

To understand this new wavefunction, we plot its Real and Imaginary components separately.

Interactive Element

import matplotlib.pyplot as plt
import numpy as np

# Defining functions
def psi_contour(x,dk): return (np.sin(dk*x)/(np.sqrt(np.pi*dk)*x))
def psi(x,k,dk): return psi_contour(x,dk)*(np.cos(k*x)+np.sin(k*x)*1j)

# Reading the input variables from the user
#print 'To understand this new wavefunction, we plot the Real and Imaginary components separately'
k = float(input("Enter the value for k_o (in Angstroms-1) = "))
dk = float(input("Enter the value for Delta k (in Angstroms-1) = "))
xmax = float(input("Enter the maximum value for x (in Angstroms) = "))

# Generating the wavefunction graph
lim1 = np.sqrt(dk/np.pi)
plt.rcParams.update({'font.size': 18, 'font.family': 'STIXGeneral', 'mathtext.fontset': 'stix'})
x = np.linspace(-xmax, xmax, 900)
fig, axes = plt.subplots(1, 2, figsize=(13,4))
str1=r"$k_o  \pm \Delta k$ = "+str(k)+r" $\pm$ "+str(dk)+r" A$^{-1}$"
# axes[0] is the graph at the left
axes[0].axis([-xmax,xmax,-1.1*lim1,1.1*lim1])
axes[0].plot(x, psi(x,k,dk).real, label="Real", color="blue", linewidth=1.8)
axes[0].plot(x, psi_contour(x,dk), label="", linestyle ="--",color="blue", linewidth=1.2)
axes[0].plot(x, -psi_contour(x,dk), label="", linestyle ="--",color="blue", linewidth=1.2)
axes[0].hlines(0.0, -xmax, xmax, linewidth=1.8, linestyle='--', color="black")
axes[0].set_xlabel(r'$x$ (Angstroms)')
axes[0].set_ylabel(r'$Real(\Psi_{\Delta k}(x))$')
axes[0].set_title('Real contribution to $\Psi_{\Delta k}(x)$ \n for '+str1)
# axes[1] is the graph at the right
axes[1].axis([-xmax,xmax,-1.1*lim1,1.1*lim1])
axes[1].plot(x, psi(x,k,dk).imag, label="Imag.", color="red", linewidth=1.8)
axes[1].plot(x, psi_contour(x,dk), label="", linestyle ="--",color="red", linewidth=1.2)
axes[1].plot(x, -psi_contour(x,dk), label="", linestyle ="--",color="red", linewidth=1.2)
axes[1].hlines(0.0, -xmax, xmax, linewidth=1.8, linestyle='--', color="black")
axes[1].yaxis.set_label_position("right")
axes[1].set_xlabel(r'$x$ (Angstroms)')
axes[1].set_ylabel(r'$Imag(\Psi_{\Delta k}(x))$')
axes[1].set_title('Imaginary contribution to $\Psi_{\Delta k}(x)$ \n for '+str1)

# Generating the probability density graph
fig, ax = plt.subplots()
ax.axis([-xmax,xmax,0.0,lim1*lim1*1.1]) # Defining the limits to be plot in the graph
ax.plot(x, (psi(x,k,dk).real)**2+(psi(x,k,dk).imag)**2,label="Probability Density", color="green") # Plotting the probability density
# Now we define labels, legend, etc
ax.legend(loc=2);
ax.set_xlabel(r'$x$ (Angstroms)')
ax.set_ylabel(r'$\left|\Psi_{\Delta k}(x)\right|^2$')
plt.title('Probability Density \n for '+str1)
plt.legend(bbox_to_anchor=(1.1, 1), loc=2, borderaxespad=0.0)

# Show the plots on the screen once the code reaches this point
plt.show()

Assuming \(k_o > \Delta k\), the \(\cos(k_o x)\) and \(\sin(k_o x)\) components oscillate with a larger frequency (\(k_o\)) and they are modulated by a sinusoidal component with smaller frequency (\(\Delta k\)).

If \(\Delta k\) is very small, the particle is still delocalized. As \(\Delta k\) increases, the particle becomes more localized in a specific region around \(x = 0\). This leads us to conclude that as the uncertainty in the momentum increases (larger \(\Delta k\)), the uncertainty in the position decreases, and the particle becomes more localized.

The superposition of waves can be accomplished by giving different contributions to different values of momenta.

For example, we can construct the superposition with a Gaussian weighting function. In this case, we modulate the contribution of each k-wave with a value given by a Gaussian distribution (maximum contribution from \(k_o\), with smaller contributions from waves with other \(k\) values).

\[\Psi_{\Delta k}(x) =N\int^{+\infty}_{-\infty}{e^{-\frac{(k-k_o)^2}{2\Delta k^2}} \ e^{ikx}}dk, \nonumber\]

where \(\Delta k\) is related to the width of the Gaussian distribution.

After integration, we obtain:

\[\Psi_{\Delta k}(x) = \sqrt{2\pi}N\Delta k e^{-\frac{1}{2}x^2\Delta k^2} e^{ik_ox} \nonumber\]

where \(N\) is a normalization constant and is equal to

\[N=\frac{1}{\sqrt{2\Delta k\sqrt{\pi^3}}} \nonumber\]

Interactive Element

import matplotlib.pyplot as plt
import numpy as np

# Defining functions
def psi_contour(x,dk): return dk*np.exp(-0.5*x*x*dk*dk)/np.sqrt(dk*np.sqrt(np.pi))
def psi(x,k,dk): return psi_contour(x,dk)*(np.cos(k*x)+np.sin(k*x)*1j)

# Reading the input variables from the user
print ("We plot the Real and Imaginary components separately")
k = float(input("Enter the value for k_o (in Angstroms-1) = "))
dk = float(input("Enter the value for Delta k (in Angstroms-1) = "))
xmax = float(input("Enter the maximum value for x (in Angstroms) = "))

# Generating the wavefunction graph
lim1 = dk/np.sqrt(dk*np.sqrt(np.pi))
plt.rcParams.update({'font.size': 18, 'font.family': 'STIXGeneral', 'mathtext.fontset': 'stix'})
x = np.linspace(-xmax, xmax, 900)
fig, axes = plt.subplots(1, 2, figsize=(13,4))
str1=r"$k_o  \pm \Delta k$ = "+str(k)+r" $\pm$ "+str(dk)+r" A$^{-1}$"
# axes[0] is the graph at the left
axes[0].axis([-xmax,xmax,-1.1*lim1,1.1*lim1])
axes[0].plot(x, psi(x,k,dk).real, label="Real", color="blue")
axes[0].plot(x, psi_contour(x,dk), label="", linestyle ="--",color="blue", linewidth=1.8)
axes[0].plot(x, -psi_contour(x,dk), label="", linestyle ="--",color="blue", linewidth=1.8)
axes[0].hlines(0.0, -xmax, xmax, linewidth=1.8, linestyle='--', color="black")
axes[0].set_xlabel(r'$x$ (Angstroms)')
axes[0].set_ylabel(r'$Real(\Psi_{\Delta k}(x))$')
axes[0].set_title('Real contribution to $\Psi_{\Delta k}(x)$ \n for '+str1)
# axes[1] is the graph at the right
axes[1].axis([-xmax,xmax,-1.1*lim1,1.1*lim1])
axes[1].plot(x, psi(x,k,dk).imag, label="Imag.", color="red")
axes[1].plot(x, psi_contour(x,dk), label="", linestyle ="--",color="red", linewidth=1.8)
axes[1].plot(x, -psi_contour(x,dk), label="", linestyle ="--",color="red", linewidth=1.8)
axes[1].hlines(0.0, -xmax, xmax, linewidth=1.8, linestyle='--', color="black")
axes[1].set_xlabel(r'$x$ (Angstroms)')
axes[1].set_ylabel(r'$Imag(\Psi_{\Delta k}(x))$')
axes[1].set_title('Imaginary contribution to $\Psi_{\Delta k}(x)$ \n for '+str1)

# Generating the probability density graph
fig, ax = plt.subplots()
ax.axis([-xmax,xmax,0.0,lim1*lim1*1.1]) # Defining the limits to be plot in the graph
ax.plot(x, (psi(x,k,dk).real)**2+(psi(x,k,dk).imag)**2,label="Probability Density", color="magenta") # Plotting the probability density
# Now we define labels, legend, etc
ax.legend(loc=2);
ax.set_xlabel(r'$x$ (Angstroms)')
ax.set_ylabel(r'$\left|\Psi_{\Delta k}(x)\right|^2$')
plt.title('Probability Density \n for '+str1)
plt.legend(bbox_to_anchor=(1.1, 1), loc=2, borderaxespad=0.0)

# Show the plots on the screen once the code reaches this point
plt.show()

We can compare the effect of the Gaussian distribution with the equally-weigthed \(k\) values:

Here are a couple questions to consider graphically:

Questions

• Q1: What do you expect to see for a particle with a momentum (\(k_o\)) with contributions from waves with very different momenta? (large value of \(\Delta_k\))
• Q2: What do you expect to see for a a particel with a momentum (\(k_o\)) with contributions from waves with similar momenta? (small value of \(\Delta_k\))

Interactive Element

import matplotlib.pyplot as plt
import numpy as np

# Defining functions
def psi_contour(x,dk): return np.sin(dk*x)*np.sin(dk*x)/(np.pi*dk*x*x)
def psi_contourG(x,dk): return dk*dk*np.exp(-x*x*dk*dk)/(dk*np.sqrt(np.pi))

# Reading the input variables from the user
print ("What is the probability density of finding the particle as a function of position?")
k = float(input("Enter the value for k_o (in Angstroms-1) = "))
dk = float(input("Enter the value for Delta k (in Angstroms-1) = "))
xmax = float(input("Enter the maximum value for x (in Angstroms) = "))

# Generating the probability density graphs
lim0 = dk/np.pi
lim1 = dk*dk/(dk*np.sqrt(np.pi))
plt.rcParams.update({'font.size': 18, 'font.family': 'STIXGeneral', 'mathtext.fontset': 'stix'})
x = np.linspace(-xmax, xmax, 900)
fig, axes = plt.subplots(1, 2, figsize=(13,4))
str1=r"$k_o \pm \Delta k$ = "+str(k)+r" $\pm$ "+str(dk)+r" A$^{-1}$"
axes[0].plot(x, psi_contour(x,dk), label="", linestyle ="-",color="green", linewidth=1.8)
axes[0].hlines(0.0, -xmax, xmax, linewidth=1.8, linestyle='--', color="black")
axes[0].axis([-xmax,xmax,-0.1*lim0,1.1*lim0])
axes[0].set_xlabel(r'$x$ (Angstroms)')
axes[0].set_ylabel(r'$\left|\Psi_k(x)\right|^2$')
axes[0].set_title('Probability Density for equally weigthed k \n '+str1)
axes[1].plot(x, psi_contourG(x,dk), label="", linestyle ="-",color="magenta", linewidth=1.8)
axes[1].hlines(0.0, -xmax, xmax, linewidth=1.8, linestyle='--', color="black")
axes[1].axis([-xmax,xmax,-0.1*lim1,1.1*lim1])
axes[1].set_xlabel(r'$x$ (Angstroms)')
axes[1].set_ylabel(r'$\left|\Psi_k(x)\right|^2$')
axes[1].set_title('Probability Density for Gaussian-weigthed k \n '+str1)

# Show the plots on the screen once the code reaches this point
plt.show()

What can we conclude: From the graphical representation we can learn that for \(\Delta k\) very small compared to \(k_o\) (try plotting for \(k_o = 4\) and \(\Delta k = 0.04\)), the particle is still delocalized. As \(\Delta k\) increases the particle becomes more localized in a specific region around \(x=0\) (i.e \(k_o = 4\) and \(\Delta k = 0.4\)). When \(\Delta k\) is large in comparison to \(k_o\) (i.e. \(k_o = 4\) and \(\Delta k = 16\)), we observe the probability density concentrated around \(x=0\).

This leads us to conclude that as the uncertainty in the momentum increases (larger \(\Delta k\)), the uncertainty in the position decreases, and the particle becomes more localized. This is the exactly the behavior predicted by the Heisenberg's uncertainty principle.

What other properties can we measure?

A Free particle has no potential (no acting force), thus the Hamiltonian only has a Kinetic term. The Schrödinger equation is:

\[-\frac{\hbar^2}{2m} \frac{\partial^2\psi(x)}{\partial x^2} =E\psi(x)\]

Since we already know an eigenfunction of the momentum operator, let's work a bit with commutator properties, and how they can be applied to the Free Particle wavefunction.

One commutator property says that if two Hermitian operators (\(\hat A\) and \(\hat B\)) commute, then they share a common set of eigenfunctions.

Let's now check if the Hamiltonian and the momentum operators commute: is

\(\left[\hat P, \hat H \right] =0\)?.

with \(\hat H =\frac{\hat P^2}{2m}\). This required evaluating

\[\begin{align*} \left[\hat P, \hat H \right] &= \hat P \cdot \hat H - \hat H \cdot \hat P \\[4pt] &=\hat P \cdot \frac{\hat P^2}{2m} - \frac{\hat P^2}{2m}\cdot \hat P \\[4pt] &= \frac {1}{2m}\left(\hat P \cdot \hat P^2-\hat P \cdot \hat P^2 \right) = 0 \end{align*}\]

The two operators commute, leding to the existance of a complete set of eigenfunctions common to both operators. Since we already know the eigenfunctions of $\hat P$, we can try to see if they are also eigenfunctions of \(\hat H\). For \(\left|k\right> = \frac{1}{\sqrt{2\pi}}e^{ikx}\), let's evaluate what is the result of \(\hat H \left|k\right>\):
(note that we switch to a notation where \(\Psi_k(x)= \left|k\right>\)) :

\[\begin{align*} \hat H\left|k\right> &= -\frac{\hbar^2}{2m} \frac{\partial^2\left|k\right>}{\partial x^2} = -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2}\left(\frac{1}{\sqrt{2\pi}}e^{ikx} \right) \\[4pt] &= -\frac{\hbar^2}{2m\sqrt{2\pi}} \frac{\partial^2e^{ikx}}{\partial x^2} = -\frac{\hbar^2}{2m\sqrt{2\pi}} \left(ik\right)^2 e^{ikx} \\[4pt] &= \frac{\hbar^2k^2}{2m} \left(\frac{1}{\sqrt{2\pi}}e^{ikx} \right) = \frac{p^2}{2m}\left|k\right> \end{align*}\]

The first and last terms in this derivation show that the momentum eigenfunctions are also eigenfunctions of the Hamiltonian, with eigenvalues equal to the classical form of the kinetic energy. However, one has to be beware that although the value of the energy might look classical, the behavior of the particles are very different. For the quantum mechanical systems, the particle can never be at rest (if \(\Delta k = 0\) the particle is deloclaized; if the particle is very localized, then its momentum has many different values).

We are now ready to tackle "A Particle in a box with infinite-potential walls"