In previous discussions, I explored the theory of the Schrödinger equation and its analytical solutions for specific cases such as the particle in a box. The objective of this page is to describe a generic method to approximate the Schrödinger equation when the potential is not simple and requires a numerical solution.
The time-dependent Schrödinger equation (TDSE) in one dimension is given by:
i\hbar \frac{\partial \psi(x,t)}{\partial t} = \left[ -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x,t) \right] \psi(x,t) = \mathbf H \psi(x,t)
where:
Analytical solutions to the TDSE are often possible only for simple potentials. However, for more complex potentials, we must resort to numerical methods. This article focuses on a numerical approach to solve the TDSE, providing a robust framework for approximating the wave function \psi(x,t) under arbitrary potential landscapes.
The first step is to discretize the right-hand side (RHS) of the Schrödinger equation, which involves the Hamiltonian operator \mathbf{H}.
Using the finite-difference scheme, the second derivative of the wave function \psi with respect to x can be approximated as:
\frac{d^2 \psi(x,t)}{dx^2} \approx \frac{\psi(x + \Delta x, t) - 2\psi(x, t) + \psi(x - \Delta x, t)}{\Delta x^2}
This allows us to discretize the Hamiltonian operator in terms of a band matrix.
The corresponding Hamiltonian matrix for a 1D system, considering the kinetic and potential energy terms, is given by:
\mathbf{H} = \begin{bmatrix} V_1 - \frac{\hbar^2}{m \Delta x^2} & \frac{\hbar^2}{2m \Delta x^2} & 0 & \cdots & 0 \\ \frac{\hbar^2}{2m \Delta x^2} & V_2 - \frac{\hbar^2}{m \Delta x^2} & \frac{\hbar^2}{2m \Delta x^2} & \cdots & 0 \\ 0 & \frac{\hbar^2}{2m \Delta x^2} & V_3 - \frac{\hbar^2}{m \Delta x^2} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & V_N - \frac{\hbar^2}{m \Delta x^2} \end{bmatrix}
Where V_i represents the potential energy at position i, and \Delta x is the spatial step size.
The corresponding Python code to construct this Hamiltonian matrix using sparse matrices is:
# Discretize the space
x = np.arange(x_min, x_max, dx)
N = len(x)
# Potential energy array (can vary)
V = create_potential(x)
# Create potential energy matrix
V = sparse.diags(V_array, 0)
# Create kinetic energy matrix
T = sparse.diags([1, -2, 1], [-1, 0, 1], shape=(N, N)) * (-hbar**2 / (2 * m * dx**2))
# Hamiltonian matrix (kinetic + potential)
H = T + V
This code constructs the Hamiltonian matrix \mathbf{H} as a sparse matrix, which includes both the kinetic and potential energy terms.
For any numerical problem, boundary conditions need to be specified. Since we are dealing with a second derivative, we need to specify both \psi(0) and \frac{d\psi}{dx} at the borders. In this formulation, we assume both are zero. This assumption is reasonable if the boundaries are relatively far from the region of interest, as \psi must go to zero as x approaches infinity to ensure the wave function is normalizable (the probability must sum to 1).
Therefore, our boundary conditions are:
\begin{aligned} & \psi(0, t) = 0 \\ & \frac{d\psi}{dx} \bigg|_{x=0} = 0 \end{aligned}
These conditions ensure that the wave function and its derivative are zero at the boundaries, which is essential for maintaining the normalizability of the wave function.
However, this is not always the case. For example, in crystal lattices, periodic boundary conditions are often used:
\begin{aligned} & \psi(0, t) = \psi(L, t) \\ & \frac{d\psi}{dx} \bigg|_{x=0} = \frac{d\psi}{dx} \bigg|_{x=L} \end{aligned}
Where L is the length of the crystal. These conditions ensure that the wave function and its derivative are periodic, reflecting the repeating structure of the crystal lattice.
Having dealt with the spatial part, it is now possible to discretize the time. The time-dependent Schrödinger equation can be expressed as:
i\hbar \frac{\partial \psi(x,t)}{\partial t} = \mathbf{H} \psi(x,t)
To numerically solve this, we discretize the time using a time step \Delta t. One common method for time discretization is the Crank-Nicolson method, which is an implicit method known for its stability and accuracy. The Crank-Nicolson method averages the Hamiltonian at the current and next time steps.
The time evolution operator can be defined using the following equations:
\psi(t + \Delta t) = \psi(t) + \Delta t \frac{\partial \psi}{\partial t}
Using the finite-difference approximation, we can write:
\psi^{n+1} = \psi^n + \Delta t \frac{\partial \psi}{\partial t}
To solve the above equation numerically, we need to formulate it in a matrix form. The Crank-Nicolson method gives us:
\left( \mathbf{I} - \frac{i\Delta t}{2\hbar} \mathbf{H} \right) \psi^{n+1} = \left( \mathbf{I} + \frac{i\Delta t}{2\hbar} \mathbf{H} \right) \psi^n
Where \mathbf{I} is the identity matrix. This can be written as:
\mathbf{A} \psi^{n+1} = \mathbf{B} \psi^n
Where:
\begin{aligned} & \mathbf{A} = \mathbf{I} - \frac{i\Delta t}{2\hbar} \mathbf{H} \\ &\mathbf{B} = \mathbf{I} + \frac{i\Delta t}{2\hbar} \mathbf{H} \end{aligned}
The corresponding Python code for this approach is:
# Convert to dense matrix for operations
H_dense = H.toarray()
# Identity matrix
In = np.eye(N)
# Time evolution operators (Crank-Nicolson method)
A = In - 1j * dt / (2 * hbar) * H_dense
B = In + 1j * dt / (2 * hbar) * H_dense
# Initial wave function (example: Gaussian)
x = np.linspace(0, N*dx, N)
psi = np.exp(-0.5 * ((x - N*dx/2) / (0.1*N*dx))**2)
psi = psi / np.linalg.norm(psi) # Normalize
# Time evolution
for _ in range(steps):
psi = np.linalg.solve(A, B @ psi)
This code constructs the matrices \mathbf{A} and \mathbf{B} for the Crank-Nicolson method and evolves the wave function \psi in time by solving the matrix equation at each time step. By using np.linalg.solve(A, B @ psi), we efficiently and stably compute the next time step’s wave function without explicitly inverting the matrix \mathbf{A}.
For increased precision in time evolution, higher-order methods such as Runge-Kutta can be employed. The scipy.integrate.solve_ivp function allows for such methods, including RK23, RK45, and DOP853 (an 8th-order method). These methods provide more accurate solutions at the cost of increased computational effort.
The following code snippet demonstrates how to implement time evolution using these higher-order methods:
from scipy import integrate
# Define the time evolution function for solve_ivp
def schrodinger_rhs(t, psi, H):
dpsi_dt = 1j / hbar * (H @ psi)
return dpsi_dt
# Define the time span and evaluation points for performing many steps with a single calculation
t_span = [0, steps * dt]
t_eval = np.linspace(0, steps * dt, steps + 1)
# Solve the Schrödinger equation using a higher-order Runge-Kutta method
sol = integrate.solve_ivp(
schrodinger_rhs, t_span, psi, method='DOP853', # 'RK23', 'RK45', or 'DOP853' are available
t_eval=t_eval, args=(H,)
)
# Update psi to the last computed value
psi = sol.y[:, -1]
The solve_ivp function in scipy requires a function that returns the right-hand side of the differential equation. In this case, it is the schrodinger_rhs function. The time span and evaluation points are defined to span the desired time range, and the method parameter allows for the selection of different Runge-Kutta methods.
By using solve_ivp, we can achieve higher accuracy in the time evolution of the wave function. The optional methods (RK23, RK45, DOP853) offer various levels of precision and computational cost, allowing to balance accuracy and performance according to the needs.
For the time-dependent Schrödinger equation, boundary conditions must also be specified for the time evolution of the wavefunction. Typically, the initial wavefunction \psi(x, 0) is provided, and the wavefunction’s evolution is determined by the Hamiltonian operator. Unlike the spatial boundary conditions, which often assume \psi(0) = 0 and \frac{d\psi}{dx}(0) = 0 for simplicity, the time evolution boundary conditions ensure the continuity and smoothness of \psi(x, t) over time.
In contrast, the wave equation, which describes the propagation of waves such as sound or light, often requires initial conditions for both the displacement \psi(x, 0) and the velocity \frac{\partial \psi}{\partial t}(x, 0). These conditions specify the state of the system at t = 0 and are necessary to solve the second-order time derivative present in the wave equation. While the Schrödinger equation involves a first-order time derivative, leading to a complex-valued wavefunction evolution driven by the Hamiltonian, the wave equation’s second-order time derivative results in real-valued solutions representing physical displacements or field intensities over time. This fundamental difference in the nature of the time derivatives underscores the distinct physical phenomena governed by these equations.
To effectively simulate the dynamics of a quantum particle, a wavepacket can be used as the initial wavefunction. A Gaussian wavepacket is particularly useful as it allows for precise control over the initial momentum and position of the particle.
The wavepacket is defined as:
\psi(x, 0) = \left( \frac{1}{2\pi\sigma^2} \right)^{1/4} \exp\left(-\frac{(x - x_0)^2}{4\sigma^2}\right) \exp\left(i\frac{p x}{\hbar}\right)
Where: - x_0 is the initial position - \sigma is the width of the wavepacket - p is the initial momentum - \hbar is the reduced Planck constant
The following Python function creates a Gaussian wavepacket with a given initial position, width, and momentum:
import numpy as np
def create_wavepacket(x):
# Gaussian shape
gaussian = np.exp(-(x - p.x0) ** 2 / (2 * p.sigma ** 2))
# Adding the wave component (momentum)
psi0 = gaussian * np.exp(-1j * p.p * x / p.hbar)
# Normalize
psi0 /= np.sqrt(np.sum(np.abs(psi0) ** 2) * p.dx)
return psi0
Additionally, the kinetic energy of the wavepacket can be computed using the following function, which calculates the second derivative of the wavepacket and integrates the kinetic energy density over space:
def T_wavepacket(psi):
d2psi_dx2 = np.gradient(np.gradient(psi, p.dx), p.dx)
T_density = -0.5 * p.hbar ** 2 / p.m * np.conj(psi) * d2psi_dx2
T = np.sum(T_density.real) * p.dx
return T
Using this wavepacket allows for precise initialization of momentum and provides a clear visualization of the wave evolution. By computing the average initial momentum of the wavepacket, insights can be gained into expected outcomes when interacting with potential barriers.
I prepared various potential energy functions, which can be used in simulations to study the behavior of a quantum particle under different conditions. The create_potential function generates these potentials based on the specified type.
Free Space:
When the potential type is set to 0, the function returns a zero potential across the entire spatial domain, simulating a free particle with no external forces acting on it.
def create_potential(x):
if p.potential == 0:
return np.zeros(len(x))
Harmonic Oscillator: For the harmonic oscillator (potential type 1), the potential is given by V(x) = \frac{1}{2} m \omega^2 x^2. This potential describes a particle in a quadratic well, which is a common model in quantum mechanics.
def create_potential(x):
if p.potential == 1:
return 0.5 * p.m * p.omega**2 * x ** 2
Infinite High Barrier:
When the potential type is 2, the function creates an infinite barrier centered around x = 0 with symmetric width defined. Outside this region, the potential is set to a very high value, effectively creating an impenetrable barrier and therefore leaving the wavepacket confined within the well.
def create_potential(x):
if p.potential == 2:
V = np.zeros(len(x))
V[x <= -p.Vx_bar] = p.Vx_bar * 1000
V[x >= p.Vx_bar] = p.Vx_bar * 1000
return V
Infinite Right Barrier:
For the infinite right barrier (potential type 3), the function creates a potential that is zero until the start of the barrier and a fix value after until the right end of the boundary, simulating a step on the right side of the domain.
def create_potential(x):
if p.potential == 3:
V = np.zeros(len(x))
V[x >= p.Vx_bar] = p.V_barrier
return V
Finite Right Barrier:
The finite right barrier (potential type 4) is created by setting a fix potential value between two spatial points, simulating a barrier of finite height and width.
def create_potential(x):
if p.potential == 4:
V = np.zeros(len(x))
V[(x >= p.Vx_bar) & (x <= p.Vx_finite_bar)] = p.V_barrier
return V
By using this function, various potential landscapes can be simulated, enabling comprehensive analysis and comparison with theoretical predictions.
The program allows for the visualization of different aspects of the wavefunction during the simulation. Depending on the selected parameters, the probability density, real and imaginary parts of the wavefunction, or the modulus and phase of the wavefunction can be visualized.
Probability Density: The probability density |\psi|^2 can be visualized. This represents the likelihood of finding a particle in a particular region of space.
|\psi(x)|^2 = \bar \psi(x) \cdot \psi(x)
ax.plot(x, abs(self._psi.conjugate() * self._psi))
Wavefunction (Real \Re{(\psi)} and Imaginary \Im{(\psi)} Parts):
# Plot Modulus
ax.plot(x, abs(psi))
# Plot Real Part
ax.plot(x, self._psi.real)
# Plot Imaginary Part
ax.plot(x, self._psi.imag)
Wavefunction (Modulus and Phase):
Alternatively, the modulus and phase of the wavefunction can be visualized. The phase is mapped to the HSL color space.
\begin{aligned} & \text{Modulus: } |\psi(x)| \\ & \text{Phase: } \arg(\psi(x)) \end{aligned}
The phase is normalized to the range [0, 1] for HSL mapping.
\text{Normalized Phase: } \frac{\arg(\psi) + \pi}{2\pi}
# Plot Modulus
ax.plot(x, abs(psi))
# Extract the phase of psi
phase = np.angle(psi)
# Normalize phase to [0, 1] for HSL
normalized_phase = (phase + np.pi) / (2 * \pi)
# Create vertices for PolyCollection
verts = [np.column_stack([x, np.zeros_like(x)])]
verts.append(np.column_stack([x, abs(psi)]))
# Update the PolyCollection with phase colors
polys = [np.column_stack(
[[x[j], x[j + 1], x[j + 1], x[j]],
[0, 0, abs(psi[j + 1]), abs(psi[j])]])
for j in range(len(x) - 1)]
self._poly.set_verts(polys)
self._poly.set_array(normalized_phase[:-1])
By setting the appropriate visualization parameters, it is possible to display:
This flexibility allows for a comprehensive analysis of the wavefunction’s behavior over time, providing valuable insights into the dynamics of the quantum system.
The complete code and results for the simulation of the time-dependent Schrödinger equation in 1D, including various potential types and visualization options, are available on GitHub. The repository contains detailed instructions on how to set up the environment, run the simulations, and visualize the results.
The repository is available at the following link here and includes: