Quantum Tunneling Through a Potential Barrier

Quantum tunneling is one of the most striking manifestations of quantum mechanics with no classical analog. In classical physics, a particle cannot pass through a potential barrier if its energy is less than the barrier height. However, in quantum mechanics, there is a non-zero probability for particles to “tunnel” through barriers that would be insurmountable according to classical physics.

This phenomenon has numerous practical applications:

In this lecture, we’ll explore quantum tunneling by numerically solving the time-dependent Schrödinger equation for a wavepacket approaching a potential barrier. We’ll implement two different numerical methods:

  1. The Crank-Nicolson method: A finite difference approach that is unconditionally stable
  2. The Split-Step Fourier method: An efficient approach that leverages fast Fourier transforms

The Time-Dependent Schrödinger Equation

The time-dependent Schrödinger equation governs the evolution of quantum states:

\[\begin{equation} 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) \end{equation}\]

Where:

  • \(\Psi(x,t)\) is the wavefunction
  • \(\hbar\) is the reduced Planck constant
  • \(m\) is the particle mass
  • \(V(x,t)\) is the potential energy function

This equation can be represented in both position space (as above) and momentum space. The transformation between these representations is accomplished using the Fourier transform.

Position and Momentum Space Representations

The wavefunction in momentum space, \(\tilde{\Psi}(k,t)\), is related to the position-space wavefunction through the Fourier transform:

\[\begin{equation} \tilde{\Psi}(k,t) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \Psi(x,t)e^{-ikx} dx \end{equation}\]

And the inverse transformation:

\[\begin{equation} \Psi(x,t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}\tilde{\Psi}(k,t)e^{ikx}dk \end{equation}\]

The Schrödinger equation in momentum space becomes:

\[\begin{equation} i\hbar \frac{\partial \tilde{\Psi}}{\partial t} = \frac{\hbar^{2}k^{2}}{2m}\tilde{\Psi} + V\left(i\frac{\partial}{\partial k}\right) \tilde{\Psi} \end{equation}\]

This duality between position and momentum space will be particularly useful in our Split-Step method.

Method 1: Crank-Nicolson Scheme

The Crank-Nicolson method is a finite difference approach that offers second-order accuracy in time and is unconditionally stable. It effectively balances the forward and backward Euler methods.

Discretization of the Schrödinger Equation

Starting from the time-dependent Schrödinger equation, we can approximate the time derivative using centered differences:

\[ \frac{\Psi(x,t+\Delta t)-\Psi(x,t)}{\Delta t} \approx -\frac{i}{\hbar} \left[\frac{1}{2}H\Psi(x,t) + \frac{1}{2}H\Psi(x,t+\Delta t)\right] \]

Where \(H = \frac{-\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x)\) is the Hamiltonian operator.

Rearranging this equation:

\[ \left(1-\frac{\Delta t}{2i\hbar} H\right)\Psi(x,t+\Delta t) = \left(1+\frac{\Delta t}{2i\hbar} H\right)\Psi(x,t) \]

And solving for \(\Psi(x,t+\Delta t)\):

\[ \Psi(x,t+\Delta t) = \left(1-\frac{\Delta t}{2i\hbar} H\right)^{-1}\left(1+\frac{\Delta t}{2i\hbar} H\right)\Psi(x,t) \]

This gives us a formula for propagating the wavefunction forward in time.

Setting Up the Simulation Domain

First, we define our spatial domain and time step:

Initial Wavepacket

We create a Gaussian wavepacket with an initial momentum (represented by \(k_0\)) directed toward the barrier:

Potential Barrier Setup

We create a potential barrier that the wavepacket will encounter:

Constructing the Hamiltonian Matrix

We represent the Hamiltonian operator as a tridiagonal matrix:

Time Evolution Matrices

Following the Crank-Nicolson scheme, we create the matrices needed for time evolution:

Animation Setup

Now we’ll prepare for animating the wavepacket’s evolution:

Animation Execution

The animation will show the wavepacket approaching, partially reflecting from, and partially tunneling through the barrier:

# Setup the canvas for animation
canvas = Canvas(width=800, height=300, sync_image_data=False)
display(canvas)

# Run the animation
for i in range(200):
    # Evolve the wavefunction
    psi = evolution_matrix.dot(psi)
    prob = np.abs(psi)**2

    # Normalize
    norm = np.sum(prob)*dx
    prob /= norm
    psi /= np.sqrt(norm)

    # Update the plot
    fig.canvas.restore_region(background)
    ax.draw_artist(points)
    points.set_data(x, prob)

    fig.canvas.blit(ax.bbox)
    X = np.array(fig.canvas.renderer.buffer_rgba())

    with hold_canvas(canvas):
        canvas.clear()
        canvas.put_image_data(X)
        sleep(0.02)

Quantum Tunneling Animation (Crank-Nicolson)

Method 2: Split-Step Fourier Method

The Split-Step Fourier method leverages the efficiency of Fast Fourier Transforms (FFT) to solve the Schrödinger equation. It works by separately handling the kinetic and potential energy operators in their respective spaces where they are diagonal.

Mathematical Foundation

We can decompose the Hamiltonian into two parts:

\[\begin{equation} \hat{H} = \hat{T} + \hat{V} \end{equation}\]

Where:

  • \(\hat{T} = \frac{-\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\) is the kinetic energy operator
  • \(\hat{V} = V(x)\) is the potential energy operator

The formal solution to the Schrödinger equation is:

\[\begin{equation} \Psi(x, t+\Delta t) = e^{-i\hat{H}\Delta t/\hbar}\Psi(x, t) \end{equation}\]

Using the Trotter-Suzuki approximation, we can split this exponential:

\[\begin{equation} e^{-i(\hat{T}+\hat{V})\Delta t/\hbar} \approx e^{-i\hat{V}\Delta t/\hbar}e^{-i\hat{T}\Delta t/\hbar} + \mathcal{O}(\Delta t^2) \end{equation}\]

The Trotter-Suzuki approximation is a fundamental technique in quantum physics that allows us to handle non-commuting operators in exponential form. When we have a Hamiltonian that is the sum of two operators \(\hat{H} = \hat{A} + \hat{B}\) that don’t commute (\([\hat{A},\hat{B}] \neq 0\)), we cannot simply write \(e^{-i(\hat{A}+\hat{B})t} = e^{-i\hat{A}t}e^{-i\hat{B}t}\).

Instead, the first-order Trotter formula gives us:

\[e^{-i(\hat{A}+\hat{B})\Delta t} = e^{-i\hat{A}\Delta t}e^{-i\hat{B}\Delta t} + \mathcal{O}(\Delta t^2)\]

For better accuracy, we can use the second-order symmetric splitting:

\[e^{-i(\hat{A}+\hat{B})\Delta t} = e^{-i\hat{A}\Delta t/2}e^{-i\hat{B}\Delta t}e^{-i\hat{A}\Delta t/2} + \mathcal{O}(\Delta t^3)\]

In the context of the Schrödinger equation, this allows us to separately handle the kinetic energy operator (\(\hat{T}\)) and potential energy operator (\(\hat{V}\)) in their respective spaces where they are diagonal, making numerical solutions much more efficient.

This gives us:

\[\begin{equation} \Psi(x, t+\Delta t) \approx e^{-i\hat{V}\Delta t/\hbar}e^{-i\hat{T}\Delta t/\hbar}\Psi(x, t) \end{equation}\]

The key insight is that \(\hat{T}\) is diagonal in momentum space, while \(\hat{V}\) is diagonal in position space. This leads to the algorithm:

  1. Apply the potential evolution operator: \(\Psi_1(x) = e^{-i\hat{V}\Delta t/\hbar}\Psi(x, t)\)
  2. Transform to momentum space: \(\tilde{\Psi}_1(k) = \mathcal{F}[\Psi_1(x)]\)
  3. Apply the kinetic evolution operator: \(\tilde{\Psi}_2(k) = e^{-i\hat{T}\Delta t/\hbar}\tilde{\Psi}_1(k)\)
  4. Transform back to position space: \(\Psi(x, t+\Delta t) = \mathcal{F}^{-1}[\tilde{\Psi}_2(k)]\)

Where \(\mathcal{F}\) represents the Fourier transform.

In quantum mechanics, an operator is diagonal in a particular basis if it acts on states in that basis by simply multiplying them by a scalar value.

Potential Energy Operator in Position Space: The potential energy operator \(\hat{V}\) is defined by its action on a wavefunction:

\[\hat{V}\Psi(x) = V(x)\Psi(x)\]

This is a simple multiplication operation in position space. If we represent the wavefunction as a vector of values at different positions, the potential operator becomes a diagonal matrix with \(V(x_i)\) as its diagonal elements. This is why we can apply the potential evolution operator directly by multiplication:

\[e^{-i\hat{V}\Delta t/\hbar}\Psi(x) = e^{-iV(x)\Delta t/\hbar}\Psi(x)\]

Kinetic Energy Operator in Momentum Space: The kinetic energy operator in position space is a differential operator:

\[\hat{T} = \frac{-\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\]

This is not diagonal in position space (it’s a second derivative). However, in momentum space, it takes a simple form:

\[\hat{T}\tilde{\Psi}(k) = \frac{\hbar^2 k^2}{2m}\tilde{\Psi}(k)\]

In momentum space, the kinetic energy operator simply multiplies each momentum component by \(\frac{\hbar^2 k^2}{2m}\), making it diagonal. This is why we can efficiently apply the kinetic evolution operator in momentum space:

\[e^{-i\hat{T}\Delta t/\hbar}\tilde{\Psi}(k) = e^{-i\hbar k^2\Delta t/(2m)}\tilde{\Psi}(k)\]

The Split-Step method exploits this property by applying each operator in the space where it’s diagonal, using Fourier transforms to switch between spaces.

Setting Up the Simulation Domain

Potential Barrier Setup

Initial Wavepacket

Setting Up the Fourier Transform Space

Precomputing Phase Factors

Both the potential and kinetic evolution operators involve phase factors that remain constant throughout the simulation. We can precompute them:

Animation Setup

Animation Execution

The animation will show the wavepacket tunneling through the potential barrier:

# Setup the canvas for animation
canvas = Canvas(width=800, height=300, sync_image_data=False)
display(canvas)

# Run the animation
for i in range(100):
    fig.canvas.restore_region(background)
    ax.draw_artist(points)

    # Perform multiple time steps between animation frames for smoother results
    for j in range(10):
        # Apply potential evolution
        psi_modx = psi_modx * phase_x

        # Transform to momentum space
        psi_modk = fft(psi_modx)

        # Apply kinetic evolution
        psi_modk = psi_modk * phase_k

        # Transform back to position space
        psi_modx = ifft(psi_modk)

        # Convert back to standard wavefunction
        psi = psi_modx * np.exp(1j * k[0] * x) * np.sqrt(2 * np.pi) / dx

        # Normalize
        prob = np.abs(psi)**2
        norm = np.sum(prob) * dx
        psi = psi / np.sqrt(norm)
        prob = prob / norm

    # Update the plot
    points.set_data(x, prob)

    fig.canvas.blit(ax.bbox)
    X = np.array(fig.canvas.renderer.buffer_rgba())

    with hold_canvas(canvas):
        canvas.clear()
        canvas.put_image_data(X)
        sleep(0.01)

Quantum Tunneling Animation (Split-Step)

Physical Interpretation and Analysis

What’s Happening in the Simulation?

In both methods, we observe:

  1. Incident Wavepacket: The Gaussian wavepacket approaches the barrier from the left
  2. Reflection and Transmission: Upon reaching the barrier, the wavepacket splits:
    • A portion is reflected back (reflection)
    • A portion passes through the barrier (tunneling)
  3. Interference: After reflection, the incident and reflected components interfere, creating oscillations in the probability density

Factors Affecting Tunneling Probability

The tunneling probability depends on several factors:

  1. Barrier Height (\(V_0\)): Higher barriers reduce tunneling probability exponentially
  2. Barrier Width (\(a\)): Wider barriers reduce tunneling probability exponentially
  3. Particle Energy: Higher energy particles (larger \(k_0\)) have higher tunneling probabilities
  4. Particle Mass: Heavier particles have lower tunneling probabilities

Connection to Quantum Theory

For a plane wave with energy \(E = \frac{\hbar^2k^2}{2m}\) incident on a rectangular barrier of height \(V_0\) and width \(a\), the transmission coefficient in the case \(E < V_0\) is:

\[\begin{equation} T \approx \frac{16E(V_0-E)}{V_0^2} e^{-2\kappa a} \end{equation}\]

Where \(\kappa = \sqrt{\frac{2m(V_0-E)}{\hbar^2}}\) is the attenuation coefficient inside the barrier.

This exponential dependence on barrier width explains why quantum tunneling is primarily observed at microscopic scales.

Comparison of Numerical Methods

Both the Crank-Nicolson and Split-Step methods accurately simulate quantum tunneling, but they have different characteristics:

Crank-Nicolson Method

  • Pros: Unconditionally stable, second-order accurate in both space and time
  • Cons: Requires matrix inversion, which can be computationally expensive for large systems

Split-Step Fourier Method

  • Pros: Very efficient for periodic boundary conditions, leverages fast FFT algorithms
  • Cons: Requires care with boundary conditions, first-order accurate in time (though higher-order splitting schemes exist)

Further Explorations

Here are some suggestions for further exploration:

  1. Momentum Space Analysis: Examine how the wavepacket’s momentum distribution changes during tunneling
  2. Barrier Parameters: Investigate how tunneling probability changes with barrier height and width
  3. Wavepacket Parameters: Explore how the initial momentum and width of the wavepacket affect tunneling
  4. Different Potentials: Try other potential shapes, such as:
    • Double barriers (resonant tunneling)
    • Smooth potentials (e.g., Gaussian barriers)
    • Time-dependent barriers
  5. Quantum Reflection: Study reflection from potential wells instead of barriers

Where to go from here

  1. Modify the code to calculate and plot the transmission and reflection coefficients as a function of incident momentum
  2. Implement a double-barrier potential and observe resonant tunneling
  3. Visualize the wavefunction in momentum space alongside position space
  4. Implement a higher-order Split-Step method (e.g., second-order Strang splitting)
  5. Add measurement simulation to observe wavefunction collapse