Visualizing The Potential Barrier Problem Using Finite Difference Time Domain Method

1. Introduction

For this project I wanted to visualize the time development of the potential barrier problem. I found that one way to do this, that did not need too much code or insight, was using the finite difference time domain method (FDTD). By using this method both quantum tunneling and quantum reflection can be visualized.

2.1 Potential Barrier

The potential barrier is a problem in quantum mechanics that demonstrates the phenomena: quantum tunneling and reflection. A free particle travels through one-dimensional space and approaches a potential. Classically there will be two scenarios:

     - If E < V, the particle will bounce back from the potential.

     - If E > V, the particle will pass the potential.

However, in quantum mechanics there is a probability that the opposite will happen, so:

     - if E < V, there is a probability that the particle will pass the potential.

     - if E > V, there is a probability that the particle bounce back.

The energy of a free particle is given by:

                                    bilde                                  (1)

2.2 Finite Difference Time Domain Method

FDTD is usually a method for solving problems in electrodynamics, where equations are solved iteratively using discrete derivative operators. However, it can also be used for solving the time-dependent Schrodinger equation (TDSE).

                         bilde                          (2)

To avoid imaginary numbers, the wave function is separated into one real part and one imaginary part:

                              bilde                              (3)

Then TDSE is divided into two real functions corresponding to each part:

                            bilde                             (4)

                             bilde                            (5)

TDSE is dependent of x and t and their derivatives, so to implement the FDTD both x and t must be discretized. This is done by separating each spatial point by a distance Δx, and each temporal point by a difference Δt. The wave function is then set in a grid with points (xi, tn).

                                    bilde                                    (6)

                                    bilde                                    (7)

And can be written:

                                   bilde                                  (8)

The derivatives in (4) and (5) are then found by central differencing scheme that gives:

                              bilde                              (9)

                             bilde                              (10)

                               bilde                                (11)

                            bilde                             (12)

Where the imaginary part is half an temporal interval from the real part. This gives the two TDSEs:

                  bilde                 (13)

                   bilde                  (14)

Which can be found iteratively. These are also called the update functions of the simulation since it moves the wave functions in the spatial direction for each time step. Here Ψ(i+1) is called the future state, Ψ(i) is called the present state and Ψ(i-1) is called the past state. ΨI is first found since it is used to calculate ΨR. Equation (13) and (14) takes on a simpler form by introducing the two constants c1 and c2.

                       bilde                     (15)

                       bilde                     (16)


                                bilde                                  (17)

                                   bilde                                   (18)

It turns out that a bigger Δt, gives a more unstable simulation. Therefore a sufficient value has to be found so that the simulation is guaranteed to be stable. This is done by adding constraints to equation (14) and (15).

Solutions to TDSE of a free particle takes the form:

                                bilde                                 (19)

With D=1 and using Euler,s formula, this gives a real and imaginary part with iterative values of x and t:

                             bilde                                (20)

                             bilde                                (21)

The notation is made easier by introducing the constants A=kiΔx - ωnΔt, B=kΔx and C=ωΔt which gives:

                                bilde                                 (22)

                             bilde                               (23)

                              bilde                               (24)

                            bilde                              (25)

By putting these into equation (15) a constraint can be found to the discrete values Δx and Δt.

               bilde                (26)

To maintain a stable simulation constraints have to be added on c1 and c2 such that equation (26) is only satisfied by real values of A, B and C. The left side of equation (26) has the constraints -1 ≤ cos(A-C) ≤ 1, which leads to the constraint[1]:

                                bilde                                  (27)

that gives:

                                 bilde                                 (28)

By using a value for Δt that satisfy this, the simulation is guaranteed to be stable. It is also worth mentioning that the computational cost grows as Δt gets smaller, so it is not necessarily better to choose small values.

Writing this as code gives a one-dimensional simulation of a free particle traveling against a potential barrier through a domain where 0 ≤ n ≤ N and 0 ≤ t ≤ T. The wave function is distributed as a gaussian wave packet given by:

                               bilde                                (29)

The wave function is plotted as the probability given by equation (30), it is also normalized so that the total probability is set to be equal to 1.

                                 bilde                                  (30)

For convenience ℏ and the mass of the particle is set to be equal to one since the simulations will only focus on whether the particle energy is greater or less than the magnitude of the potential barrier, not the exact values of the energies.

3. Results

The project files are available here:

By changing the values marked as simulation variables in main.cpp, the simulation gets different properties. Only T, Tcheck and E was changed for this project. Where T is the total amount of time steps, Tcheck is amount of time steps before plotting the wave functions and E is the particle energy.

One simulation was done for each of the five situations:

     - E << V

     - E < V

     - E = V

     - E > V

     - E >> V

A wave packet with higher energy will travel faster than one with lower energy, so different values of T and Tcheck in main.cpp was used to keep the animations around the same duration for each simulation. Therefore, the velocity of each wave packet in the animations is not comparable to their actual velocities, or to each others.

E << V


E < V


E = V


E > V


E >> V


4. Conclusion

The outcome of each simulation turned out as presumed. For big energy differences the outcome was more or less as expected from a classical particle, where it either got reflected by the potential or passed through it with no probability of the opposite happening. For small energy differences the outcome was as expected from a quantum particle, with a probability present of both being reflected and transmitted.

The wave functions were normalized, but while checking the total probability after every iteration, the value was found to fluctuate between 0.98 and 1.02, and thus the simulation is not entirely correct. The reason for this is probably that the precision of double values are not accurate enough, giving the deviation from 1.

Looking at the simulations, it is noticeable that the peak of the wave packet shrinks with time, while the width increases. This is expected[2] as the width should increase quadratically as a function of time.

To conclude, the finite difference time domain method proved a good way to visualize the potential barrier problem.


[1] Nagel, J. R. (2009) A Review and Application of the Finite-Difference Time-Domain Algorithm Applied to the Schrodinger Equation, Aces Journal, Vol. 24, No. 1, Available at: Link

[2] Unknown Evolution of a free Gaussian wavepacket, Ohio State University, Available at: Link

Soriano, A. (2004) Analysis of the finite difference time domain technique to solve the Schrodinger equation for quantum devices, Journal of Applied Physics, Available at: Link

Hemmer, P. C. (2015) Kvantemekanikk, Fagbokforlaget.