**Pseudo-spectral method** for solving the** time-dependent Schrödinger equation** are powerful numerical techniques that leverage the efficiency of spectral methods for spatial discretization combined with suitable time-stepping schemes. The aim is to approximate the solution of the Schrödinger equation, which describes the quantum mechanical behavior of a particle or system of particles. The time-dependent Schrödinger equation in one dimension can be written as:

\[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 \(i\) is the imaginary unit, \(\hbar\) is the reduced Planck’s constant, \(m\) is the mass of the particle, \(V(x)\) is the potential energy as a function of position, and \(\psi(x, t)\) is the wave function of the particle, which contains all the information about the system’s state.

## Steps for Pseudo-Spectral Method

**1. Spatial Discretization:** The first step involves discretizing the spatial domain into grid points. The wave function \(\psi(x, t)\) and potential \(V(x)\) are represented on these grid points. Spectral methods typically use Fourier transforms or Chebyshev polynomials for the spatial representation, allowing differential operators to be applied in the spectral (transformed) domain.

**2. Spectral Representation:** Apply a Fourier transform (or another suitable spectral transform) to \(\psi(x, t)\) to move from the spatial domain to the spectral domain. This transformation makes the application of the differential operator (related to the kinetic energy part of the Schrödinger equation) straightforward because differentiation in the spatial domain corresponds to multiplication in the spectral domain.

For a Fourier spectral method, the transformation is given by:

\[\hat{\psi}(k, t) = \mathcal{F}\{\psi(x, t)\} = \int_{-\infty}^{\infty} e^{-ikx}\psi(x, t) dx\]

and the inverse transform is:

\[\psi(x, t) = \mathcal{F}^{-1}\{\hat{\psi}(k, t)\} = \frac{1}{2\pi}\int_{-\infty}^{\infty} e^{ikx}\hat{\psi}(k, t) dk\]

**3. Time Evolution:** The Schrödinger equation is split into its kinetic (\(T\)) and potential (\(V\)) energy components. The time evolution can be treated using operator splitting methods, allowing the wave function to evolve under each part separately within a time step \(\Delta t\). A common approach is to use the Crank-Nicolson scheme or exponential time-differencing methods for this purpose.

For the kinetic part, in the spectral domain, the evolution is given by:

\[\hat{\psi}(k, t+\Delta t) = e^{-i\frac{\hbar k^2}{2m}\Delta t}\hat{\psi}(k, t)\]

For the potential part, in the spatial domain, the evolution is given by:

\[\psi(x, t+\Delta t) = e^{-i\frac{V(x)}{\hbar}\Delta t}\psi(x, t)\]

**4. Iterative Solution: **The solution is advanced in time by alternating between evolving the wave function under the kinetic and potential parts, typically using a leap-frog or other symmetric time-stepping method to ensure stability and accuracy. After each time step, if necessary, an inverse Fourier transform is applied to return to the spatial representation.

**5. Boundary Conditions and Normalization:** Care must be taken to apply appropriate boundary conditions, often periodic or infinite potential walls, and to normalize the wave function at each time step to ensure physical validity.

## Summary

This is a basic outline of using pseudo-spectral methods to solve the time-dependent Schrödinger equation. In practice, the choice of spectral basis, time-stepping scheme, and treatment of boundary conditions can significantly affect the accuracy and efficiency of the solution. The pseudo-spectral method’s advantage lies in its ability to accurately capture the wave function’s evolution with relatively few grid points compared to finite difference methods, especially in systems with smooth potential functions.

## An Example – Quantum Tunneling

To simulate quantum tunneling through a step barrier for an electron using the pseudo-spectral method in a 1D system, we’ll write a Python code that follows these steps:

**Initialize parameters**: Set up the physical constants, simulation parameters, and the step barrier potential.**Prepare the initial wave function**: Use a Gaussian wave packet.**Implement the pseudo-spectral method**: Use FFT for spatial derivatives and split-step method for time evolution.**Visualize the results**: Plot the probability density to observe the tunneling effect.

The initial Gaussian wave packet, which is a localized wave function representing a quantum particle with a certain probability of being found at different positions, (x), and is taken as:

\[\psi(x,0) = \frac{1}{\sqrt{\sigma\sqrt{\pi}}} e^{-\frac{(x-x_0)^2}{2\sigma^2}} e^{ik_0x}\]

where:

*x*is the initial position of the center of the packet,_{0}*σ*is the width of the packet,*k*is the initial wave number corresponding to the particle’s momentum._{0}

**Here’s a simplified Python code that demonstrates this process:**

**##### Import Libraries**
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, ifft
from scipy.constants import hbar, pi
**##### Parameters**
m = 9.10938356e-31 # electron mass
L = 5e-8 # size of the spatial domain (50 nm)
N = 1024 # number of spatial points
dx = L / N # spatial resolution
x = np.linspace(0, L, N) # spatial grid
V0 = 40*1.6e-19 # height of the step barrier (60 eV)
barrier_width = L / 20 # width of the step barrier
barrier_start = L / 2 - barrier_width / 2 # start of the barrier
dt = 1e-18 # time step (10 fs)
t_max = 5e-15 # total time of evolution (7 ps)
k0 = 4.5e10 # initial wave number
x0 = L/3 # initial position of the wave packet
sigma0 = L / 50 # width of the Gaussian wave packet (narrower)
**##### Potential**
V = np.zeros_like(x)
V[(x >= barrier_start) & (x <= barrier_start + barrier_width)] = V0
**##### Initial wave function: Gaussian wave packet**
psi_x = np.exp(-0.5 * ((x - x0) / sigma0) ** 2) * np.exp(1j * k0 * x)
psi_x /= np.sqrt(np.sum(np.abs(psi_x) ** 2) ) # normalize
**##### Precompute kinetic and potential operators for efficiency**
k = 2 * pi * np.fft.fftfreq(N, d=dx)
T_operator = np.exp(-1j * hbar * k**2 / (2 * m) * dt)
V_operator = np.exp(-1j * V / hbar * dt)
**##### Time evolution**
for t in range(int(t_max / dt)):
# Split-step Fourier method
psi_x = ifft(T_operator * fft(psi_x)) # Kinetic evolution
psi_x *= V_operator # Potential evolution
Rho = np.abs(psi_x) ** 2
if t % 100== 0:
# Plotting the probability density |psi|^2 at specific times
plt.figure(figsize=(10, 6))
plt.plot(x*1e9, Rho, label=f"t = {t*dt*1e15:.2f} ps")
plt.plot(x*1e9, 0.05*V/V0, label="Normalized Potential Barrier", color='red') # Plot V to show barrier
plt.xlabel('x (nm)')
plt.ylabel(r'$|\psi(x)|^2$ (Arb. unit)')
plt.title('Quantum Tunneling Simulation')
plt.legend()
plt.show()

## Simulation Results

We have explored the quantum tunneling phenomenon by an electron having kinetic energy *E _{0}* = 78 eV across different barrier heights: (

**1)**,

*V*= 30 eV_{0}**(ii)**, and

*V*= 60 eV_{0}**(iii)**. The animation below showcases the evolution of the wave packet.

*V*=80 eV_{0}### Case 1: Time evolution of the wave packet encountering a barrier with a height of *V*_{0} = 40 eV:

_{0}

### Case 2. Time evolution of the wave packet encountering a barrier with a height of *V*_{0} = 60 eV:

_{0}

### Case 3. Time evolution of the wave packet encountering a barrier with a height of *V*_{0} = 80 eV:

_{0}