Numerical Solution Of The 1d Time-Dependent Schrödinger Equation

Numerical solution of the 1d time-dependent Schrödinger equation

In my exploration of the time-dependent Schrödinger equation in one dimension, I focus on numerical methods to handle cases where the potential V(x, t) complicates analytical solutions. This approach is essential for understanding quantum mechanical phenomena in more realistic scenarios, such as variable potential landscapes in quantum wells or during particle interaction scenarios.

Mathematical formulation

The time-dependent Schrödinger equation in one dimension is represented as:

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)

where \psi(x,t) is the wave function, \hbar is the reduced Planck constant, m is the mass of the particle, and V(x,t) represents the potential energy, which can be complex and time-varying.

Spatial discretization

To approximate the solution, I discretize the Hamiltonian operator \mathbf{H} using a finite difference method for the spatial derivatives. The second derivative is approximated by:

\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 construction of the Hamiltonian matrix in terms of kinetic and potential energy components, crucial for evolving the wave function in time.

Numerical implementation

Here is how I set up the Hamiltonian matrix in Python using sparse matrices for efficiency:

x = np.linspace(x_min, x_max, N)  # Discretized space
V = np.vectorize(potential_function)(x)  # Potential function evaluated at each point

# Kinetic energy matrix using finite differences

T = sparse.diags([1, -2, 1], [-1, 0, 1], shape=(N, N)) * (-hbar**2 / (2 * m * dx**2))

# Hamiltonian matrix

H = T + sparse.diags(V, 0)

Boundary conditions

For my simulations, I use zero boundary conditions, assuming the wave function vanishes at the boundaries of the domain. This is typical for isolated systems where the particle is not expected to be at the edges.

Temporal discretization

The temporal aspect of the Schrödinger equation is handled using the Crank-Nicolson method, which is semi-implicit and offers stability and accuracy:

from scipy.linalg import solve

I = np.eye(N)  # Identity matrix
A = I - 1j * dt / (2 * hbar) * H.toarray()
B = I + 1j * dt / (2 * hbar) * H.toarray()

# Initial wave function
psi = np.exp(-0.5 * ((x - x0) ** 2) / (sigma ** 2))
psi = psi / np.linalg.norm(psi)  # Normalization

for t in range(num_steps):
    psi = solve(A, B @ psi)

Optionally an higher order explicit Runge-Kutta can be selected instead:

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]

Visualization

Visualizing the probability density |\psi|^2 and the real and imaginary parts of \psi helps in understanding the dynamics and evolution of the wavefunction. This is key for interpreting physical phenomena like tunneling and interference patterns.

import matplotlib.pyplot as plt

plt.figure()
plt.plot(x, np.abs(psi)**2, label='Probability Density')
plt.plot(x, psi.real, label='Real Part')
plt.plot(x, psi.imag, label='Imaginary Part')
plt.legend()
plt.show()

Conclusion

The techniques I discussed provide a robust framework for simulating the Schrödinger equation in scenarios where analytical solutions are unfeasible. This is crucial for advancing our understanding of quantum systems in complex environments.

For more insights into this topic, you can find the details here.

If you’re interested in experimenting with the code and running your own simulations, you can access the full source on GitHub here.