Numerical solution of the 2d time-dependent Schrödinger equation
In this blog post, I focus on extending the time-dependent Schrödinger equation (TDSE) to a two-dimensional system, an essential step towards simulating more realistic quantum systems. My previous work on 1D quantum dynamics laid the foundation, and this extension enables me to tackle more complex potential landscapes that vary in both spatial dimensions.
The time-dependent Schrödinger equation in two dimensions is expressed as:
i\hbar \frac{\partial \psi(x,y,t)}{\partial t} = \left[ -\frac{\hbar^2}{2m} \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right) + V(x,y,t) \right] \psi(x,y,t)
Here, \psi(x,y,t) is the wave function, \hbar is the reduced Planck constant, m is the mass of the particle, and V(x,y,t) represents the potential energy. Analytical solutions to this equation are rare, so I rely on numerical methods for most cases.
Spatial discretization
To numerically solve the TDSE in two dimensions, I start by discretizing the Hamiltonian operator \mathbf{H}, which now includes second derivatives with respect to both x and y. The finite difference approximations for the second derivatives are:
\begin{aligned} \frac{\partial^2 \psi(x,y,t)}{\partial x^2} \approx \frac{\psi(x + \Delta x, y, t) - 2\psi(x, y, t) + \psi(x - \Delta x, y, t)}{\Delta x^2} \\ \frac{\partial^2 \psi(x,y,t)}{\partial y^2} \approx \frac{\psi(x, y + \Delta y, t) - 2\psi(x, y, t) + \psi(x, y - \Delta y, t)}{\Delta y^2} \end{aligned}
Substituting these into the Hamiltonian, I get:
\mathbf{H}_{i,j} = \frac{-\hbar^2}{2m} \left[ \frac{\psi_{i+1,j} - 2\psi_{i,j} + \psi_{i-1,j}}{\Delta x^2} + \frac{\psi_{i,j+1} - 2\psi_{i,j} + \psi_{i,j-1}}{\Delta y^2} \right] + V_{i,j} \psi_{i,j}
This matrix form allows me to handle both kinetic and potential energy contributions effectively.
Boundary conditions
I implement three types of boundary conditions in the simulation: periodic, hard wall, and permeable wall. Periodic boundary conditions connect opposite edges of the domain, simulating an infinite, repeating space. Hard wall boundaries reflect the wave function entirely, while permeable walls allow partial transmission, achieved by adding an imaginary component to the potential.
In my implementation, the boundary conditions are incorporated directly into the Laplacian matrix. For example, periodic boundaries are handled by modifying the matrix elements at the edges:
data[1][-1::self.Ny] = data[3][0::self.Ny] = 1 / self.dx**2
data[0][-self.Ny:] = data[4][:self.Ny] = 1 / self.dy**2
This code ensures that the wave function wraps around the boundaries, preserving continuity.
Temporal discretization
For time evolution, I use the Crank-Nicolson method, which is both stable and accurate. The time-dependent Schrödinger equation in 2D is discretized as:
\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{A} = \mathbf{I} - \frac{i\Delta t}{2\hbar} \mathbf{H} and \mathbf{B} = \mathbf{I} + \frac{i\Delta t}{2\hbar} \mathbf{H}.
The time evolution is then computed by solving the equation \mathbf{A} \psi^{n+1} = \mathbf{B} \psi^n. This approach ensures that my simulation remains stable even for relatively large time steps.
Initial wavefunction
I typically initialize the wavefunction \psi(x, y, 0) with a 2D Gaussian wavepacket, which provides control over the initial position and momentum. The wavepacket is given by:
\psi(x, y, 0) = \left( \frac{1}{2\pi\sigma_x \sigma_y} \right)^{1/2} \exp\left(-\frac{(x - x_0)^2}{4\sigma_x^2} - \frac{(y - y_0)^2}{4\sigma_y^2}\right) \exp\left(i\frac{k_x x + k_y y}{\hbar}\right)
This wavepacket evolves over time according to the Schrödinger equation, allowing me to simulate the quantum dynamics of a particle under various potentials.
Visualization
Throughout the simulation, I visualize different aspects of the wavefunction, such as the probability density |\psi|^2, the real and imaginary parts, and the phase. For example, the probability density is computed as:
|\psi(x, y)|^2 = \bar{\psi}(x, y) \cdot \psi(x, y)
This provides insight into the likelihood of finding the particle at different positions in the domain.
In summary, extending the Schrödinger equation to two dimensions allows me to simulate more complex quantum systems with spatially varying po
Conclusion
The techniques I discussed provide a robust framework for simulating the Schrödinger equation in scenarios where analytical solutions are unfeasible.
For more insights into this topic, you can find the details here.
You can access the full source on GitHub here.