# Uncomment ONE line to choose matplotlib backend
# if using Jupyter Notebook, use interactive "notebook" backend for best results
# if using Jupyter Lab, use interactive "widget" backend for best results
# if both fail, use static "inline" backend
%matplotlib notebook
#%matplotlib widget
#%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
Units in numerical computations¶
An alternative to SI units is atomic units. Here, $\hbar$, the electron mass $m_e$, the elementary charge $e$ and $4\pi\varepsilon_0$ are all defined to have the value $1$. Then, energies will be in the unit hartree, $E_h \equiv \hbar^2 / m_e a_0^2 \approx 4.36 \cdot 10^{-18} \text{ J}$, and lengths will be in the unit $a_0$, the Bohr radius. Yet another option would be to use $\text{nm}$ and $\text{eV}$ for lengths and energies. 1 hartree corresponds to ca 27.2 eV.
Use the units you prefer, but be aware of the limited precision of numbers on a computer.
# Assign values to parameters that will not change, e.g., hbar = 1 etc (if you use atomic units)
Numerical solution of the time independent Schrödinger equation (TUSL)¶
To solve the Schrödinger equation $$ \hat{H} \psi = -\frac{\hbar^2}{2 m} \psi'' + V \psi = E \psi, $$ amounts to finding the energy eigenvalues $E$ and the corresponding eigenfunctions $\psi(x)$ for a given potential $V(x)$. Here, we will use an elegant and general method to solve this equation for an arbitrary potential (in one dimension). The method is also discussed in the lectures.
We assume that a particle is constrained to move between the two positions $x_0$ and $x_{N+1}$, and we discretize this interval in the points $x_0, x_1, \ldots, x_N, x_{N+1}$, with a constant distance $\Delta x$ between neighbouring points. Outisde this interval, the potential is $V(x \leq x_0) = V(x \geq x_{N+1}) = \infty$, and hence $\psi(x \leq x_0) = \psi(x \geq x_{N+1}) = 0$. In other words, the wavefunction needs to be determined on the $N$ positions $x_1, \ldots, x_N$. To each position we assign the values $\psi_i = \psi(x_i)$ and $V_i = V(x_i)$ of the wavefunctions and the potential, respectively, and we refer to the values of these functions in all positions with the vectors $\boldsymbol{V} = [V_1, \ldots, V_N]^T$ and $\boldsymbol{\psi} = [\psi_1, \ldots, \psi_N]^T$. (Here, $T$ denotes the transpose.)
An intuitive and simple approximation of the derivative of a function is the central difference $$ \psi'(x) = \frac{\psi(x + \Delta x / 2) - \psi(x - \Delta x / 2)}{\Delta x}. $$ If you use this approximation twice, we can approximate the second derivative as $$ \psi''(x) = \frac{\psi'(x + \Delta x / 2) - \psi'(x - \Delta x / 2)}{\Delta x} = \frac{\psi(x + \Delta x) - 2 \psi(x) + \psi(x - \Delta x)}{\Delta x^2}.$$ We insert this into the TUSL and obtain $$-\frac{\hbar^2}{2 m} \frac{\psi_{i+1} - 2 \psi_i + \psi_{i-1}}{{\Delta x}^2} + V_i \psi_i = E \psi_i \quad ; \quad i = 1, 2, \ldots, N.$$
We may write these $N$ equations in a compact form by introducing an $N \times N$ Hamiltonian matrix $H$ with elements $$ H_{i j} = \begin{cases} \hbar^2 / (m \Delta x^2) + V_i & \text{for} \,\, i = j & \text{(on the diagonal)} \\ -\hbar^2 / (2 m \Delta x^2) & \text{for} \,\, i = j \pm 1 & \text{(on the semidiagonals)} \\ 0 & \text{ellers} \\ \end{cases} $$ and using the vector (array) $\boldsymbol{\psi} = [\psi_1, \ldots, \psi_N]^T$. We may then write $$ H \boldsymbol{\psi} = E \boldsymbol{\psi}. $$ The energies $E$ and the eigenfunctions $\boldsymbol{\psi}$ are now eigenvalues and eigenvectors of the matrix $H$.
Numerical libraries have well tested functionality for computing eigenvalues and eigenvectors of matrices. There are also specific effective functions for diagonalizing special matrices, for example the tridiagonal (as well as real and symmetric) matrix $H$.
Construct a function that calculates all eigenvalues $E$ and the corresponding eigenvectors $\boldsymbol{\psi}$ for a particle with mass $m$ in a given potential $\boldsymbol{V}$. If you use Python, the function eigh_tridiagonal from scipy.linalg is tailored for our purposes. Check that the eigenvectors form a complete and orthonormal set.
Next, we will use this numerical method with various potentials. If possible, we will compare numerical values with analytical expressions. At this stage, we recommend that you write a rather generic function for plotting, to be used in all questions below.
Write a function that displays a potential $\boldsymbol{V}$ as well as one or more energy eigenvalues $E$ and eigenfunctions $\boldsymbol{\psi}$ (or $|\boldsymbol{\psi}|^2$) as functions of $x$. The function will also be used to compare numerical and analytical energies and wavefunctions. It is probably a good idea to make this function in parallel with the remaining questions, and adapt it to what is needed.
Particle in a box¶
This is typically the first system introduced in the lectures. The potential, the normalized eigenfunctions, and the energy eigenvalues are:
$$ V(x) = \begin{cases}0 & \text{for}\,\, 0 \leq x \leq L \\ \infty & \text{ellers} \end{cases}, \quad \psi_n(x) = \sqrt{\frac{2}{L}} \sin{\frac{n \pi x}{L}}, \quad E_n = \frac{n^2 \pi^2 \hbar^2}{2 m L^2}, \quad \quad n = 1, 2, \ldots . $$
Compare numerical and anlytical values of a few energies (e.g. in a table) and energy eigenfunctions (in a figure) for an electron in a potential box. What is the energy difference between the ground state and the first excited state if the electron is in a narrow box 1 nm wide? What is the corresponding energy difference for en elecron in a wide box 1 $\mu$m wide? Use the unit meV (millielectronvolt). What is the maximum system temperature in each case resulting in a thermal energy $k_B T$ pr particle not larger than the energy separation between the ground state and the 1. excited state?
Harmonic oscillator¶
Another standard quantum mechanical problem is the harmonic oscillator with $$ V(x) = \frac{1}{2}m \omega^2 x^2, \quad \psi_n(x) = \frac{1}{\sqrt{2^n n!}} \cdot \left(\frac{m \omega}{\pi \hbar}\right)^{1/4} \cdot \exp{\left(-\frac{m \omega x^2}{2 \hbar}\right)} \cdot H_n\left(\sqrt{\frac{m \omega}{\hbar}}x\right), \quad E_n = \left(n+\frac{1}{2}\right)\hbar \omega, \quad \quad n = 0, 1, 2, \ldots $$ The functions $H_n(y)$ (with dimensionless $y$) are called (the physicist's) Hermite polynomials. They are available in numerical libraries, but they can also be computed with the recursion relation $$ H_n(x) = 2 x H_{n-1}(x) - 2 (n-1) H_{n-2}(x), \quad H_0(x) = 1, \quad H_1(x) = 2 x.$$
The harmonic oscillator is particularly interesting in two-particle systems, where a problem involving e.g. two atoms with mass $m_1$ and $m_2$ in a diatomic molecule can be reduced to an equivalent one-particle problem with a reduced mass $m = m_1 m_2 / (m_1 + m_2)$. The reduced mass $m$ and the angular frequency $\omega$ determine the value of the spring constant that governs the vibrational movement of the two atoms.
Compare numerical and analytical values of a few energies and energy eigenfunctions in a harmonic oscillator. The effective spring constant in a nitrogen molecule (N$_2$) is 2.3 kN/m. What is the corresponding angular frequency $\omega$? Compare the energy difference $\hbar\omega$ between the quantized vibrational states of the molecule with the available thermal energy $k_B T$ pr molecule at room temperature. Argue why the molar heat capacity $C_V$ for a nitrogen gas at room temperature is $5R/2$, and not $7R/2$ as expected from the classical equipartition theorem. (Here, $R$ is the molar gas constant.)
Particle in finite potential wells¶
A third well known example is a single potential well $$V(x) = \begin{cases} -V_0 & \text{for}\,\, 0 < x < w \\ 0 & \text{otherwise} \end{cases}$$ with width $w$ and depth $V_0 > 0$. In this potential, there are no simple analytical solutions for the bound state energy eigenvalues.
The single well is easily generalized to a periodic potential consisting of $N_w$ such wells placed next to each other, with barriers of width $b$ (and potential $V=0$) between the wells. This is a simple model of the periodic potential experienced by an electron in a solid material with a regular crystalline structure, relevant for metals, semiconductors and insulators.
Display the bound states for an electron, both in a single well and in a periodic potential consisting of several wells. Choose values of the well depth and width so that you obtain 3-4 bound states in the single well. How are the energy levels distributed in the periodic potential when compared with the single well? Based on your results, can you explain what we mean by the band structure of a solid material, using concepts like band width and bandgap?
Voluntary extra exercise:
Use (for example) a simple box potential with a large width and zero potential and construct a gaussian and normalized initial state $\Psi(x,0)$ centered in a suitable position $x_0$ and with a suitable group velocity and corresponding momentum. Express $\Psi(x,0)$ as a linear combination of (all) stationary states and compute the coefficients $c_n$ ($n=1,2,\dots,N$). Next, compute $\Psi(x,t)$ for selected times $t$ and display the wave packet at each of these times. Use Re $\Psi$, Im $\Psi$, $|\Psi|$ or $|\Psi|^2$. Does the wave packet move back and forth between the two hard walls more or less like a classical particle? For the ambitious: Make an animation of the movement of the wave packet.