Wavepacket Scattering: Transmission & Reflection

Hello! Welcome to our next lesson on numerical methods in quantum mechanics.

In our last session, we explored Variational Monte Carlo, a powerful stochastic method for finding the ground state properties of a time-independent system. We saw how to use the Metropolis algorithm to sample a trial wavefunction and estimate the ground state energy.

Today, we shift our focus from static properties to dynamics. We will tackle the problem of how a quantum particle, described by a localized wavepacket, interacts with a potential barrier over time. This is a fundamental process in many areas of physics, from alpha decay to electron transport in semiconductor devices.

Your learning outcome for this lesson is to numerically simulate wavepacket scattering from a potential barrier to extract transmission and reflection coefficients. To achieve this, you will learn how to solve the time-dependent Schrödinger equation (TDSE) on a computer and then analyze the results to obtain physical observables.

1. From Continuous Evolution to Discrete Steps

The starting point for any quantum dynamics simulation is the time-dependent Schrödinger equation. In one dimension, it is:
$i\hslash \frac{\partial \psi (x,t)}{\partial t}=\hat{H}\psi (x,t)=\left[-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}+V(x)\right]\psi (x,t)$

To solve this on a computer, we must discretize both space and time. We replace the continuous variables $x$ and $t$ with a grid of points $x_j=j\Delta x$ and $t_n=n\Delta t$. The wavefunction $\psi (x,t)$ becomes a set of values ${\psi }_j^n=\psi (x_j,t_n)$.

Our goal is to find a reliable rule to calculate the state of the system at the next time step, ${\psi }^{n+1}$, given the current state, ${\psi }^n$. A particularly robust and widely used method for this is the Crank-Nicolson algorithm. It's an implicit method, which means it solves for the future state ${\psi }^{n+1}$ using a system of equations involving both ${\psi }^n$ and ${\psi }^{n+1}$. This gives it excellent stability properties and, crucially for quantum mechanics, it is unitary, meaning it conserves the total probability $\int |\psi {|}^{2}dx$ over time.

The core result of the Crank-Nicolson method is that the time-evolution is reduced to solving a linear system of equations, ${\psi }^{n+1}=U_1^{-1}{U}_2{\psi }^n$, at each time step. Since the matrices $U_1$ and $U_2$ are tridiagonal, this system can be solved very efficiently, making the method practical.

2. Setting Up a Realistic Simulation

Before we can run the simulation, we need to address a few practical but crucial details.

The Initial State: A Gaussian Wavepacket

Instead of an unphysical, infinitely spread-out plane wave, we'll start with a localized particle. A Gaussian wavepacket is a common and convenient choice:
$\psi (x,0)={\left(\frac{1}{2\pi {\sigma }_0^2}\right)}^{1/4}{e}^{-(x-x_0{)}^2/(4{\sigma }_0^2)}{e}^{i{k}_{0}x}$
This describes a particle initially centered at $x_0$, with a spatial width characterized by ${\sigma }_0$, and moving with an average momentum $\hslash k_0$. Your simulation will start by setting up an array representing $\psi (x,0)$ on your spatial grid.

The Boundaries: Avoiding Unphysical Reflections

A major challenge in simulating scattering on a finite grid is preventing the wavepacket from reflecting off the artificial boundaries of your simulation domain. If the wave reflects from the edge and comes back, it will interfere with the part of the wavepacket that was reflected from the potential barrier, contaminating your results.

A sophisticated solution is to use Absorbing Boundary Conditions (ABCs). These are special mathematical conditions applied at the grid's edges that are designed to "soak up" any incoming wave, effectively mimicking an infinite space.

Numerical Stability and Accuracy

Even with a stable algorithm like Crank-Nicolson, the accuracy of your simulation depends on your choice of grid spacing $\Delta x$ and time step $\Delta t$. They must be small enough to resolve the most rapid oscillations in the wavefunction, both in space and time.

• Spatial Grid ($\Delta x$): The grid must be fine enough to resolve the shortest de Broglie wavelength in the simulation. The wavelength is given by $\lambda =h/p=h/\sqrt{2mK}$, where $K$ is the kinetic energy. Higher kinetic energy means a shorter wavelength, requiring a smaller $\Delta x$.
• Time Step ($\Delta t$): The time step must be short enough to resolve the fastest time evolution, which is governed by the total energy $E$ through the relation $E/\hslash$.

3. Running the Simulation and Extracting Physics

After setting up the initial state, potential, and boundary conditions, you will run the time-evolution loop, repeatedly solving ${\psi }^{n+1}=U_1^{-1}{U}_2{\psi }^n$ for many time steps. The output will be a large array, PSI[j, n], containing the value of the wavefunction at every grid point j for every time step n. A good way to visualize this is with a heatmap of the probability density, $|\psi (x,t){|}^2$.

Once the simulation has run long enough for the interaction to complete, how do we calculate the transmission ($T$) and reflection ($R$) coefficients?

There are two common methods:

1. Probability Current: A more rigorous method is to calculate the probability current, $J(x,t)=\frac{\hslash }{m}\text{Im}({\psi }^{\ast }\frac{\partial \psi }{\partial x})$, at a point $x_{trans}$ far beyond the barrier. The total transmitted probability is the time integral of this current, $P_{trans}={\int }_0^{\infty }J(x_{trans},t)dt$. The transmission coefficient is then $T=P_{trans}/P_{initial}$, where $P_{initial}$ is the total probability of the initial packet (which should be 1 if normalized).

2. Integrated Probability: A simpler and more intuitive method is to simply wait until the wavepacket has fully separated into reflected and transmitted components. At a final time $t_f$, you integrate $|\psi (x,t_f){|}^2$ over the region to the right of the barrier to get $P_{trans}$, and over the region to the left to get $P_{refl}$. Since the initial probability was 1, we have $T=P_{trans}$ and $R=P_{refl}$.

For your implementation, the second method is more straightforward.

4. Your Task: Simulating Scattering from a Rectangular Barrier

Now it's time to put all the pieces together. Your task is to write a Python script that simulates a Gaussian wavepacket scattering from a simple rectangular potential barrier.

Here is a step-by-step guide for your implementation. Given your strong background in Python and Numpy, you are well-equipped to translate these steps into code.

1. Setup:

   • Define your physical constants ($\hslash ,m$) and simulation parameters ($L$, n_points, $\Delta x$, $\Delta t$, n_steps). Choose them carefully based on the stability considerations from section 2.
   • Create your spatial grid x and an empty array PSI to store the results.

2. Initial and Potential Conditions:

   • Define your initial Gaussian wavepacket $\psi (x,0)$ using the formula from section 2. Choose parameters $x_0,{\sigma }_0,k_0$ that place the packet far to the left of the barrier and give it enough momentum to move towards it.
   • Define your rectangular potential barrier $V(x)$. It should be zero everywhere except for a finite region where it has height $V_0$.

3. Crank-Nicolson Implementation:

   • Construct the tridiagonal matrices $U_1$ and $U_2$ as described in resource LINK. You can use scipy.sparse.spdiags to do this efficiently.
   • (Optional, but highly recommended for better results) Modify the corners of your matrices to include Absorbing Boundary Conditions.
   • Pre-compute the LU decomposition of $U_1$ using scipy.sparse.linalg.splu. This will make the time-stepping loop much faster.

4. Time-Evolution Loop:

   • Loop for n from 0 to n_steps:
     • Calculate the right-hand side of the equation: b = U2.dot(psi_current).
     • Solve for the next state: psi_next = LU.solve(b).
     • Update psi_current = psi_next.
     • Store psi_next in your PSI array for later analysis/visualization.

5. Analysis:

   • After the loop finishes, take the final wavefunction $\psi (x,t_f)$.
   • Find the indices corresponding to the regions to the left and right of the barrier.
   • Calculate $R$ by summing $|\psi {|}^2\Delta x$ over the reflection region.
   • Calculate $T$ by summing $|\psi {|}^2\Delta x$ over the transmission region.
   • As a check, verify that $T+R\approx 1$.

You can find a good set of starting parameters that are known to work in the resource Computational Solution to Discrete Time-Dependent ... (ID LINK) in the section "Simulate Wave hitting a Potential Barrier".

Conclusion

In this lesson, we have moved from the static, time-independent world of VMC to the dynamic evolution of quantum systems. You have learned the theoretical and practical steps required to build a powerful simulation tool from scratch.

Key Takeaways:

• The time-dependent Schrödinger equation can be solved numerically by discretizing space and time and using a time-stepping algorithm.
• The Crank-Nicolson method is a stable and unitary algorithm that transforms the TDSE into a sequence of matrix equations.
• Practical considerations are critical for a successful simulation: defining a localized initial state, handling grid boundaries with Absorbing Boundary Conditions, and choosing stable grid parameters ($\Delta x,\Delta t$).
• Physical observables like transmission and reflection coefficients are extracted by post-processing the final state of the simulated wavefunction.

Preview of the Next Lesson:

Today, we simulated the quantum mechanics of a single particle in continuous space. This is a powerful tool, but many problems in your field of interest, condensed matter physics, are naturally described on a lattice of atoms.

In our next lesson, we will explore a different but related framework: "Implement the Landauer-Büttiker formalism for a 1D tight-binding model to calculate quantum conductance." We will shift from a time-domain simulation to a steady-state, energy-domain picture of transport, connecting the scattering concepts we've learned today to the macroscopic property of electrical conductance. This will be a key step towards understanding quantum transport in nano-scale devices.