In this Python project, we will look at quantum tunneling. Even if you are completely new to quantum mechanics, you must have heard of this mysterious, quantum property of matter. Particles can, apparently, move through solid barriers! What is this quantum woo all about? This, we will find out by solving a standard textbook problem — that of a one-dimensional potential barrier. Our goal will be to derive the tunneling probability. We will then implement a simulation of a particle near such a potential barrier and look at the numerical solution of the problem.

## The Schrödinger equation

What do we do when we want to study the evolution of a macroscopic physical system — a bullet, a ping-pong ball, the Solar system — in time? We reach for an equation of motion, such as the second law of Newton. To describe a system in classical mechanics we only need the equation of motion and two pieces of data: the initial position, and the initial velocity of all parts of the system. Armed with these, we can predict the state of the system at any arbitrary point in time.

Quantum mechanics differs from classical mechanics in the equation of motion and the required initial conditions. A quantum particle is described not by its position and velocity, but rather by a complex function — called the wave function — spread out over all of space. Its motion (or change), in turn, is dictated by Schrödinger’s equation of motion:

The capital H with the fancy hat is the Hamiltonian operator. It is related to the Hamiltonian we discussed some time ago in relation to the classical Hamiltonian description of the double pendulum. The difference here is that instead of just a scalar number (which can be interpreted as the total energy of the system in certain cases), it is an operator — a mathematical entity that acts on a wave function and changes it, producing a different wave function.

Furthermore, the above equation is clearly a first-order partial differential equation in time. Hence, only the initial state of the system is required as the initial condition.

Finally, it would be nice to know the shape of the Hamiltonian operator. Generally, one can go from classical mechanics to quantum mechanics just by putting hats on the quantities (yes, it’s often that easy). For example, a classical particle in a spatially changing potential (think electric potential energy or gravitational potential) has the following Hamiltonian:

The corresponding quantum-mechanical Hamiltonian operator is:

Now we just plug in the expressions of the position and momentum operators in position space (read regular space): and . We find:

Our Schrödinger equation therefore looks like this:

What are its solutions ? Since the left hand side of the equation is independent of time, we could try to look for a separated solution. Moreover, the right hand side would then represent a simple first-order differential equation in time. It makes sense, then, to look for a solution in the form:

Plugging this Ansatz into our equation we have that:

Carrying out the time derivation and dividing by the complex exponential we get:

This is called the time independent Schrödinger equation. Its solutions are called eigenstates — they represent wave functions that, upon action of the Hamiltonian operator, remain identical up to a scalar multiplicative coefficient: the eigenvalue. This scalar value can in turn be identified with the total energy corresponding to the state. In simpler terms, think of the eigenstates as of the basic building blocks of the solution, while the eigenvalues are the corresponding energies. Let us now solve the Schrödinger equation for a model system illustrating the essence of quantum tunneling!

## Rectangular Potential Barrier

Consider a one-dimensional potential barrier described by a potential :

This naturally splits the physical system into three parts. Let us designate the solution for all by , for by , and for all even greater by .

The time independent Schrödinger equation for the region is:

the solution of which is trivial. In the absence of a potential, the solution is a linear combination of sines and cosines. However, we might instead prefer the complex exponential form:

Plugging this into the left hand side of our equation of motion we find that:

and thus:

The solution in region is, naturally, very similar:

Finally, we can look at the region of space inside the potential barrier. While the potential here is still constant, it is non-zero. This, however, only results in a change of the wavenumber; the solution itself is mathematically the same as in the other regions. We have that:

where:

How can we interpret the above solution? In the regions to both sides of the potential barrier, the solution represents a superposition of harmonic waves travelling in either direction. Inside the potential barrier, the solution depends on the energy of the particle compared to the barrier height. If , the particle “flies above” the barrier. This has a direct analogy in classical physics, like a ball flying above a physical wall. However, if , the wavenumber is imaginary, and the solution becomes a superposition of real exponentials. While exponentials are quickly decreasing functions, they still lead to a non-zero wave function (and thus probability density) inside the potential barrier! Particles flying *through* walls!

The solution presented just now includes several free parameters. In total, there are seven of them: , and the energy of the particle . On the other hand side, the wave function has to satisfy certain boundary conditions. For one thing, it must be a smoothly varying function of the spatial coordinate. Hence, we get two continuity conditions on either side of the barrier. For another, the derivative of a wave function must also be continuous. Both the continuity of the wave function, and that of its derivative, stem from the mathematical structure of the Schrödinger equation — it is a second order differential equation in space. As long as the potential is finite, both continuity conditions must hold. Long story short, this gives us another two boundary conditions. Finally, we have the normalization condition. Unlike in the case of spatially confined systems, where the integral of the wave function is a finite number we can normalize, here we must use a different normalization criterion. However, the exact choice does not matter. The important point is that we get a fifth constraint.

This leaves us with 7 parameters and 5 constraints. Hence, the energy spectrum of the solution is continuous and doubly degenerate. In other words: for any choice of energy , there exist two solutions, representing a particle moving in either direction.

Now, let’s see where we stand. The solution to the problem is:

The continuity of the wave function and of its derivative yields four boundary conditions:

Plugging in the expressions for the wave functions we get:

This is the system of equations we need to solve. However, this represents a very general case.

Let us consider a particle moving in from the left with a unit amplitude . We will neglect the case of a particle moving towards the barrier from the right by setting . We now rename the parameter , as it represents the reflected amplitude, and , since this corresponds to the transmitted part. We now have:

Four equations, four parameters. We can get rid of and , and find the expressions for the reflected and transmitted amplitudes. Their squares, and , give us the reflection and transmission coefficients, respectively.

As this post is getting quite long, let me stop here. In the next instalment we shall solve the above system of equations, find the two parameters and , and calculate the probability with which the particle tunnels through our potential barrier.

Until next time!