Time evolution in quantum mechanics. The particle in a box model
! WORK IN PROGRESS !
Roberto Di Remigio, Luca Frediani
In this notebook, we will consider how the time evolution of physical systems is treated in quantum mechanics. We will again make use of the simple particle in a box model to highlight the basic ideas.
The equation of motion for quantum mechanical systems is the time-dependent Schrödinger equation:
\begin{equation} H\Psi(x, t) = \mathrm{i}\hbar\frac{\partial\Psi(x, t)}{\partial t} \end{equation}where $H$ is the Hamiltonian operator and $\Psi(x, t)$ is the wavefunction describing the system. For a general, time-dependent Hamiltonian the solution of the time-dependent Schrödinger equation is far from trivial and is usually achieved using perturbation theory.
In case the Hamiltonian is not time-dependent, we can solve the time-dependent Schrödinger equation by separation of the space ($x$) and time variables ($t$):
\begin{equation} \Psi(x, t) = \psi(x)\theta(t) \end{equation}Inserting this Ansatz into the equation above:
\begin{equation} [H\psi(x)]\theta(t) = \psi(x)\mathrm{i}\hbar\frac{\partial\theta(t)}{\partial t} \end{equation}followed by division by $\psi(x)$ and $\theta(t)$ yields:
\begin{equation} \frac{H\psi(x)}{\psi(x)} = \frac{1}{\theta(t)}\mathrm{i}\hbar\frac{\partial\theta(t)}{\partial t} \end{equation}This means that the left- and right-hand sides have to be equal to a constant value, which is the energy of the system:
\begin{equation} \frac{H\psi(x)}{\psi(x)} = E = \frac{1}{\theta(t)}\mathrm{i}\hbar\frac{\partial\theta(t)}{\partial t} \end{equation}What we have obtained is thus two equations. The first one is time-independent Schrödinger equation:
\begin{equation} H\psi_n(x) = E_n\psi_n(x) \end{equation}whose eigenvalues indeed are the allowed energies for the system. The second equation determines the time evolution of the eigenfunctions:
\begin{equation} -\frac{\mathrm{i}E_n}{\hbar}\theta(t) = \frac{\partial\theta(t)}{\partial t} \end{equation}this is not an eigenvalue equation. The solution to this differential equation is given by the exponential function:
\begin{equation} \theta(t) = \exp(-\frac{\mathrm{i}E_nt}{\hbar}) = \cos(\frac{E_nt}{\hbar}) -\mathrm{i}\sin(\frac{E_nt}{\hbar}) \end{equation}Eventually, the solution to the time-dependent problem is:
\begin{equation} \Psi_n(x, t) = \psi_n(x)\theta(t) = \exp(-\frac{\mathrm{i}E_nt}{\hbar}) \psi_n(x) \end{equation}Notice that this form of the solution is general. It is the same for all quantum systems were the Hamiltonian doesn't depend explicitly on time. You will also notice that the time evolving part of the wavefunction is given as a simpe phase factor depending on the energy of the eigenfunction to the time-independent problem. This means that the probability density for the eigenstates is unchanged with the passing of time:
\begin{equation} |\Psi_n(x, t)|^2 = \Psi_n^*(x, t)\Psi_n(x, t) = \exp(+\frac{\mathrm{i}E_nt}{\hbar}) \psi^*_n(x) \exp(-\frac{\mathrm{i}E_nt}{\hbar}) \psi_n(x) = \psi^*_n(x)\psi_n(x) =|\psi(x)|^2 \end{equation}which also means that expectation values of time-independent operators are unchanged with the passing of time:
\begin{equation} \langle O \rangle = \int\mathrm{d} x \Psi_n^*(x, t)O(x)\Psi_n(x, t) = \int\mathrm{d} x \exp(+\frac{\mathrm{i}E_nt}{\hbar}) \psi^*_n(x)O(x) \exp(-\frac{\mathrm{i}E_nt}{\hbar}) \psi_n(x) = \int\mathrm{d} x \psi^*_n(x)O(x)\psi_n(x) \end{equation}States for which all the observables are unchanged by time evolution are called stationary. All eigenvectors of a time-independent Hamiltonian are stationary states.
In the following, our toy system will be the particle in a box (both one-dimensional and two-dimensional) For the one-dimensional model the (normalized) eigenfunctions are: \begin{equation} \psi_n(x) = \sqrt{\frac{2}{L}}\sin(\frac{n\pi x}{L}) \quad\quad \forall n \neq 0 \end{equation} with energies: \begin{equation} E_n = \frac{h^2n^2}{8ML^2} \quad\quad \forall n \neq 0 \end{equation}
For the two-dimensional, square model the (normalized) eigenfunctions are: \begin{equation} \psi_{nm}(x,y) = \frac{2}{L}\sin(\frac{n\pi x}{L})\sin(\frac{m\pi y}{L}) \quad\quad \forall n,m \neq 0 \end{equation} with energies: \begin{equation} E_{nm} = \frac{h^2}{8ML^2}(n^2 + m^2) \quad\quad \forall n, m \neq 0 \end{equation}
Eigenfunctions and probability densities
The first thing we can try is to plot wavefunctions and the associated probability densities at different times. Start with the one-dimensional model and look at the eigenfunctions for the ground and first excited state at different times. Use the functions defined below to achieve this.
At what time points should we make the plots? Given that these are stationary states and that the time-evolving factor in the wavefunction is an oscillatory function, there will be a characteristic oscillation period for each eigenfunction. At what time $T$ (the period) does a given eigenfunction become equal to its starting value? In other words, for which $T$ is the following equation satisfied: \begin{equation} \Psi(x, t = 0) = \Psi(x, t = T) \end{equation} or more explicitly: \begin{equation} \psi_n(x) = \exp(-\frac{\mathrm{i}E_nT}{\hbar}) \psi_n(x) \end{equation} Once you have find the value of the period, you can plot at time points that are multiples of the period.
Warning The eigenfunctions are now complex functions! Plot both the real and imaginary part to get something meaningful! Use the NumPy real and imag functions to get them.
In [14]:
import numpy as np
from scipy.constants import *
def theta(En, t):
""" The time evolving factor for a stationary state.
This is a complex function!
En -- the energy of the stationary state
t -- time point
"""
return np.exp(-1j * (En / hbar) * t)
def eigenfunction1D(n, L, x):
""" Normalized eigenfunction for the 1D particle in a box.
n -- the NumPy array with the x values
L -- the quantum number
x -- the size of the box
"""
if n <= 0:
raise TypeError('Quantum number has to be greater than 0!')
normalization = np.sqrt(2.0 / L)
return normalization * np.sin(n * pi * x / L)
def eigenvalue1D(n, M, L):
""" Eigenvalue for the 1D particle in a box.
Returned values is in Joule.
n -- the quantum number
M -- the mass of the particle
L -- the size of the box
"""
if n <= 0:
raise TypeError('Quantum number has to be greater than 0!')
gs_energy = h**2 / (8 * M * L**2)
return n**2 * gs_energy
The example below shows how things would look for the ground state eigenfunction (once you fill in how to calculate the period ;))
In [37]:
import matplotlib.pyplot as plt
# make sure we see it on this notebook
%matplotlib inline
# The mass will be that of an electron
M = m_e
# The length of the box is pi
L = pi
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].set_xlabel(r'$x$')
axes[0].set_ylabel(r'$\psi$')
axes[0].set_title(r'1D particle in a box wavefunctions')
axes[0].set_xlim([0, L])
axes[0].spines['right'].set_color('none')
axes[0].spines['top'].set_color('none')
axes[0].xaxis.set_ticks_position('bottom')
axes[0].spines['bottom'].set_position(('data',0)) # set position of x spine to x=0
axes[0].yaxis.set_ticks_position('left')
axes[0].spines['left'].set_position(('data',0)) # set position of y spine to y=0
E1 = eigenvalue1D(1, M, L)
print('Eigenvalue 1 is {} Joule'.format(E1))
# What's the period?
T1 =
print('The oscillation period is {} seconds'. format(T1))
x = np.linspace(0, L, 1000)
axes[0].plot(x, np.real(theta(E1, 0)*eigenfunction1D(1, L, x)), label=r'$\Re(\psi_1(t=0))$')
axes[0].plot(x, np.imag(theta(E1, 0)*eigenfunction1D(1, L, x)), label=r'$\Im(\psi_1(t=0))$')
axes[0].plot(x, np.real(theta(E1, 0.1*T1)*eigenfunction1D(1, L, x)), label=r'$\Re(\psi_1(t=0.1T_1))$')
axes[0].plot(x, np.imag(theta(E1, 0.1*T1)*eigenfunction1D(1, L, x)), label=r'$\Im(\psi_1(t=0.1T_1))$')
axes[0].legend()
axes[1].set_xlabel(r'$x$')
axes[1].set_ylabel(r'$|\psi|^2$')
axes[1].set_title(r'1D particle in a box probability densities')
axes[1].set_xlim([0, pi])
axes[1].spines['right'].set_color('none')
axes[1].spines['top'].set_color('none')
axes[1].xaxis.set_ticks_position('bottom')
axes[1].spines['bottom'].set_position(('data',0)) # set position of x spine to x=0
axes[1].yaxis.set_ticks_position('left')
axes[1].spines['left'].set_position(('data',0)) # set position of y spine to y=0
# We need to use the absolute function from NumPy to get the correct probability
# distribution since the wavefunction is now complex
axes[1].plot(x, (np.absolute(theta(E1, 0)*eigenfunction1D(1, L, x)))**2, label=r'$|\psi_1(t=0)|^2$')
axes[1].plot(x, (np.absolute(theta(E1, 0.1*T1)*eigenfunction1D(1, L, x)))**2, label=r'$|\psi_1(t=0.1T_1)|^2$')
axes[1].legend()
Linear combinations
The next step is to look at the time evolution of a linear combination of eigenfunctions. Let us consider the following, normalized, linear combination: \begin{equation} \Psi(x, t) = \frac{1}{\sqrt{2}}(\exp(-\frac{\mathrm{i}E_1t}{\hbar})\psi_1(x) + \exp(-\frac{\mathrm{i}E_2t}{\hbar})\psi_2(x)) \end{equation} From the expression of the eigenvalues: \begin{equation} E_n = \frac{h^2n^2}{8ML^2} \end{equation} we can see that the energy of the excited states is a multiple of that for the ground state: \begin{equation} E_\mathrm{gs} = E_1 = \frac{h^2}{8ML^2} \quad\quad E_n = n^2E_\mathrm{gs} \end{equation} We can thus rewrite the linear combination above as: \begin{equation} \Psi(x, t) = \frac{1}{\sqrt{2}}(\exp(-\frac{\mathrm{i}E_1t}{\hbar})\psi_1(x) + \exp(-\frac{\mathrm{i}4E_1t}{\hbar})\psi_2(x)) =\frac{1}{\sqrt{2}}\exp(-\frac{\mathrm{i}E_1t}{\hbar})\left[\psi_1(x) + \exp(-\frac{\mathrm{i}3E_1t}{\hbar})\psi_2(x)\right] \end{equation}
In [ ]: