Split-Step Fourier Method for TDSE: Step-by-Step Guide

Last Updated on August 18, 2024 by Max

The Time-Dependent Schrödinger Equation (TDSE) describes how the quantum state of a physical system changes over time.

Solving this equation analytically is often impossible for complex systems, so numerical methods like the Split-Step Fourier Method (SSFM) are used.

The SSFM splits the evolution of the wave function into manageable steps, involving separate evolutions for kinetic and potential energy operators.

In this guide, we’ll walk through each step of the Split-Step Fourier Method for TDSE, with clear explanations for beginners.

The Time-Dependent Schrödinger Equation (TDSE)

The 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)\right] \psi(x,t),\]

where, \( \psi(x,t) \) represents the wave function of the particle, \( m \) denotes the mass of the particle, \( \hbar \) is the reduced Planck’s constant, and \( V(x) \) is the potential energy function.

The Core Idea of the Split-Step Fourier Method

The SSFM is based on the idea that the time evolution can be split into two main parts: kinetic and potential.

The evolution is approximated by:

Half-step evolution by the potential energy operator.
Full-step evolution by the kinetic energy operator.
Another half-step evolution by the potential energy operator.

Mathematically, this is expressed as [1, 2],

\[\psi(x,t + \Delta t) \approx e^{-i \hat{V} \Delta t / 2\hbar} e^{-i \hat{T} \Delta t / \hbar} e^{-i \hat{V} \Delta t / 2\hbar} \psi(x,t),\]

where \( \hat{T} \) and \( \hat{V} \) are the kinetic and potential energy operators.

This approach allows separate and simpler handling of kinetic and potential terms. The method is second-order accurate in time step \( \Delta t \), providing a good balance between computational efficiency and accuracy.

Steps in Split-Step Fourier Method in One Dimension

Step 1: Discretization of Space and Time

Begin by discretizing space into \( N \) points with spacing \( \Delta x \) and time into intervals \( \Delta t \).

Let:

\( x_j = j\Delta x \) for \( j = 0, 1, 2, \dots, N-1 \)
\( t_n = n\Delta t \) for \( n = 0, 1, 2, \dots \)

Step 2: Initial Wave Function

Set the initial wave function \( \psi(x,0) = \psi_0(x) \). Ensure it satisfies boundary conditions and is well-behaved across the domain.

Typically, the initial wave function is chosen as a Gaussian wave packet.

Step 3: First Half-Step Evolution by the Potential Energy Operator

The first half-step evolution applies the potential energy operator:

\[\psi(x,t + \Delta t/2) = \psi(x,t) \exp\left(-i \frac{V(x)\Delta t}{2\hbar}\right)\]

This step accounts for the influence of the potential energy on the wave function.

Step 4: Full-Step Evolution by the Kinetic Energy Operator

Next, transform to momentum space using the Fourier transform:

\[\tilde{\psi}(k,t + \Delta t/2) = \mathcal{F}\{\psi(x,t + \Delta t/2)\}\]

Apply the kinetic energy operator in momentum space:

\[\tilde{\psi}(k,t + \Delta t/2) = \tilde{\psi}(k,t) \exp\left(-i \frac{\hbar k^2 \Delta t}{2m}\right)\]

Then, transform back to position space:

\[\psi(x,t + \Delta t/2) = \mathcal{F}^{-1}\{\tilde{\psi}(k,t + \Delta t/2)\}\]

Step 5: Second Half-Step Evolution by the Potential Energy Operator

Finally, apply the second half-step evolution by the potential energy operator:

\[\psi(x,t + \Delta t) = \psi(x,t + \Delta t/2) \exp\left(-i \frac{V(x)\Delta t}{2\hbar}\right)\]

This step completes the time evolution for a single interval \( \Delta t \).

Step 6: Iterate for the Desired Time Evolution

Repeat Steps 3 to 5 for successive time steps until the desired final time is reached.

Step 7: Visualization and Analysis

Visualize the probability density \( |\psi(x,t)|^2 \) and analyze the results, such as the evolution of wave packets, interference patterns, or tunneling behavior.

Example: Simulation of Quantum Tunneling

Let’s now simulate and visualize quantum tunneling of an electron through a one-dimensional square potential barrier using the Split-Step Fourier Method, as demonstrated above.

Initial Wave Function: Gaussian Wave Packet

In the simulation, an electron is represented by a Gaussian wave packet propagating towards the potential barrier, described by the following equation:

\[\psi(x, t=0) = \exp\left(-\frac{(x – x_0)^2}{2\sigma_x^2}\right) \exp(i k_0 x),\]

where \( x \) is the position, \( x_0 \) is the initial center of the wave packet, \( \sigma_x \) is the width of the Gaussian envelope, \( k_0 \) is the initial wave vector (\( k_0 = \frac{2\pi}{\lambda} \)), the Gaussian term \( \exp\left(-\frac{(x – x_0)^2}{2\sigma_x^2}\right) \) localizes the wave packet, and \( \exp(i k_0 x) \) represents its momentum.

Figure 1: Schematic diagram of the quantum tunneling simulation setup.

Potential Barrier

The potential energy \( V(x) \) as a function of position \( x \) is given by:

\[V(x) =\begin{cases}V_0 & \text{if } |x| < \frac{a}{2}, \\0 & \text{otherwise},\end{cases}\]

where \( V_0 \) is the height of the potential barrier, and \( a \) is the width of the potential barrier. In the simulation code, the height \( V_0 \) of the barrier is defined as:

\[V_0 = 0.8 \times \frac{\hbar^2 k_0^2}{2 m_e},\]

where \( \hbar \) is the reduced Planck constant, \( k_0 \) is the wave number associated with the initial wave packet, and \( m_e \) is the mass of the electron.

Values of all the simulation parameters are listed in Table I.

Parameter	Value
Reduced Planck constant \(( \hbar )\)	\( 1.0545718 \times 10^{-34} \) J·s
Electron mass \(( m_e )\)	\( 9.10938356 \times 10^{-31} \) kg
de Broglie wavelength \(( \lambda )\)	\( 5 \times 10^{-9} \) m
Wave number \(( k_0 )\)	\( 2\pi / \lambda \) m\(^{-1}\)
Simulation domain size \(( L_x )\)	\( 1 \times 10^{-6} \) m
Number of grid points \(( N_x )\)	1024
Grid spacing \(( \Delta x )\)	\( L_x / N_x \)
Time step \(( \Delta t )\)	\( 3 \times 10^{-15} \) s
Number of time steps \(( N_t )\)	1001
Wave packet width \(( \sigma_x )\)	\( 60 \times 10^{-9} \) m
Initial position \(( x_0 )\)	\( -2 \times 10^{-7} \) m
Barrier width	\( 4 \times 10^{-8} \) m
Barrier height	\( 0.8 \times \frac{\hbar^2 k_0^2}{2m_e} \) J

Table I: Values of the parameters used in the simulation.

Tracking Quantum Tunneling

The simulation tracks the wave packet’s evolution in time, showing how it partially tunnels through the barrier. The code calculates and displays the transmission probability at different time steps, along with the corresponding probability density and potential barrier. This allows the study of quantum tunneling phenomena in a visually intuitive and computationally efficient manner.

Figure 2: Distribution of the probability density \(|\psi(x,t)|^2\) after interaction with the potential barrier.

Simulation Code in Python

Wave Packet Simulation


import numpy as np
import matplotlib.pyplot as plt

# Constants
hbar = 1.0545718e-34  # Reduced Planck constant (Joule second)
m_e = 9.10938356e-31  # Mass of electron (kg)
wavelength = 5e-9  # de Broglie wavelength (meters)
k0 = 2 * np.pi / wavelength  # Wave number

# Grid parameters
Lx = 1e-6  # Simulation domain size (meters)
Nx = 1024  # Number of grid points
dx = Lx / Nx  # Grid spacing
x = np.linspace(-Lx / 2, Lx / 2, Nx)

# Time parameters
dt = 3e-15  # Time step (seconds)
Nt = 1001  # Number of time steps

# Wave packet parameters
sigma_x = 60e-9  # Width of the Gaussian packet (meters)
x0 = -2e-7  # Initial position of the packet (meters)
kx0 = k0  # Initial wave vector

# Barrier parameters
barrier_width = 4e-8  # Width of the square barrier (meters)
barrier_height = 0.8 * hbar**2 * k0**2 / (2 * m_e)  # Height of the barrier (Joules)

# Initial Gaussian wave packet
psi0 = np.exp(-(x - x0) ** 2 / (2 * sigma_x ** 2)) * np.exp(1j * kx0 * x)
psi0 = psi0.astype(np.complex128)

# Normalize the initial wave packet
psi0 /= np.sqrt(np.sum(np.abs(psi0) ** 2) * dx)

# Potential energy (Square barrier in the center)
V = np.zeros_like(x)
V[np.abs(x) < barrier_width / 2] = barrier_height

# Fourier space components
kx = np.fft.fftfreq(Nx, dx) * 2 * np.pi
K2 = kx ** 2

# Split-step Fourier method
psi = psi0.copy()
transmission_prob = 0

for t in range(Nt):
    # (a) 1/2 Evolution for the potential energy in real space
    psi *= np.exp(-1j * V * dt / (2 * hbar))
    
    # (b) Forward transform
    psi_k = np.fft.fft(psi)
    
    # (c) Full evolution for the kinetic energy in Fourier space
    psi_k *= np.exp(-1j * hbar * K2 * dt / (2 * m_e))
    
    # (d) Inverse Fourier transform
    psi = np.fft.ifft(psi_k)
    
    # (e) 1/2 Evolution for the potential energy in real space
    psi *= np.exp(-1j * V * dt / (2 * hbar))
    
    # Calculate transmission probability
    if t >= 0:
        Rho = np.abs(psi) ** 2
        transmission_prob = np.sum(Rho[x > barrier_width / 2]) * dx
    
    # Visualization of wave packet evolution
    if t % 50 == 0:
        plt.figure(figsize=(8, 6))
        Rho = np.abs(psi) ** 2
        plt.plot(x * 1e9, Rho / np.max(Rho), color='hotpink', lw=2)
        plt.plot(x * 1e9, V / np.max(V), color='black', lw=2)
        plt.title(f'Time step: {t}\nTransmission Probability: {transmission_prob:.3f}')
        plt.xlabel(r'$x$ (nm)', fontsize=14)
        plt.ylabel(r'Probability Density', fontsize=14)
        plt.pause(1)  # Pause for 1 second
        plt.close()

plt.show()

Animation of Quantum Tunneling: Get the Python Code

Split-Step Fourier Method in Two Dimension

To implement the Split-Step Fourier Method in two spatial dimensions, the method involves extending the 1D approach to both \( x \) and \( y \) directions, applying the Fourier transform in two dimensions.

This process is detailed in one of our earlier articles, “Time Evolution of a Free Particle Wave Packet: Numerical Simulation of Wave Packet Dispersion,” where the application of this method to 2D wave packet simulations is discussed [3].

Conclusion

The Split-Step Fourier Method (SSFM) is an effective numerical approach for solving the Time-Dependent Schrödinger Equation (TDSE).

By splitting the evolution into half-steps for potential energy and full-steps for kinetic energy, the SSFM efficiently handles complex quantum systems.

This guide detailed its application to one-dimensional quantum tunneling, with potential extensions to higher dimensions.

The SSFM offers intuitive visualization of quantum phenomena like tunneling and interference, making it accessible for both educational and research purposes. Its stability and efficiency make it a valuable tool for studying time-dependent quantum systems.

Additional Resources

Below is the Python code to generate the quantum tunneling animation: