Solving the Schrödinger equation numerically by expansion in eigenstates¶

Examples – Quantum Mechanics¶

By Jonas Tjemsland, Andreas Krogen og Jon Andreas Støvneng.

Last edited: January 20th 2019

---

Introduction¶

In this notebook we will be solving the one-dimensional Schrödinger equation,

$$i\hbar\frac{\partial\Psi(x, t)}{\partial t} = -\frac{\hbar^2}{2m}\frac{\partial^2 \Psi( x, t)}{\partial x^2}+V(x)\Psi( x, t) $$

numerically for an arbitrary initial condition $\Psi(x, 0)$. The eigenstates $\psi_n(x)$ and the eigenenergies $E_n$ of the system are found by solving the time-independent Schrödinger equation

$$-\frac{\hbar^2}{2m}\frac{\partial^2 \psi_n(x)}{\partial x^2}+V(x)\psi_n(x) = E_n\psi_n(x),$$

and normalizing the result. The inital condition $\Psi(x, 0)$ is expanded in terms of $\psi_n(x)$:

$$\Psi(x,0) = \sum_{i}\alpha_i\psi_i(x).$$

In turn, the solution at time $t$, $\Psi(x, t)$, is given by

$$\Psi(x, t) = \sum_n\alpha_n\psi_n(x)\exp\left(-i\frac{E_n}{\hbar}t\right).$$

As an example, we will be propagating an electron given by a gaussian wave packet towards a potential barrier. A similar example is studied in our notebook on One-Dimensional Wave Packet Propagation, but with a quite different approach.

The numerical scheme that is used is developed and explained in detail in the appendix at the end of this notebook. The reader is adviced to read through this before reviewing the notebook.

We start by importing packages, setting common figure parameters and defining physical parameters.

Potential¶

As mentioned in the introduction, we will be propagating an electron towards a potential barrier in one dimension. We will be considering a domain $x\in[0,L]$. Let us use $\Delta x = 1\text Å$, which is a typical diameter of an atom. In turn, the width of the barrier is being decided by the number of discretization points it consists of. We want each side of the potential barrier to be large, so that the electron is not influenced by the barrier or the edges at $t=0$. We choose $N=10$ discretization points for the barrier, and 50 times that for each of the sides. The barrier has a height $V_0=1.5\cdot 1.6\cdot 10^{-19}J = 1.5\text{eV}$.

Play around with other parameters and potential barriers. The code in this notebook works even for arbitrary potentials!

Wave packet¶

We will be representing the initial electron as a gaussian wavepacket, $$\Psi(x,0)=C\exp\left(-\frac{(x-x_0)^2}{4\sigma^2}+i\frac{p_0x}{\hbar}\right),$$ where $p_0=\sqrt{2mE_0}$ is the momentum of the wave packet, $E_0$ the energy of the electron, $x_0$ is the initial expectation value, and $\sigma$ is some parameter specifying the width of the wave packet.

It will not be unreasonable to choose $E_0\sim V_0$. As we will see, this will give a good visualization of transmission and reflection. We start by choosing an energy a bit higher than the potential height, $E_0=1.39V_0$. $x_0$ to be in the middle of the left part of the domain and $\sigma$ (one standard deviation) to 1/8 of the left part. Play around with different parameters!

We now visualize the initial wave packet and the potential (with a suitable scaling).

Solving the eigenvalue problem (Schrödinger equation)¶

Now that all the parameters are settled, we can finally solve the Schrödinger equation. This is done by solving an eigenvalue problem. We are using a real symmetric band matrix solver (You could of course also use numpy.linalg.eigh, but this requires the initialization of the whole matrix, mostly consisting of zeros). We thus need to initialize the diagonal and the sub- and superdiagonal. This is explained in detail in the appendices.

Note that this is the most computationally demanding part of these computations.

Let us visualize some of the eigenstates and eigenenergies!

Here is a quick question for the reader: why is $\psi_0(x)$ and $\psi_1{x}$ almost equal in the right part of the domain? Would we expect the same result for other pairs $(\psi_n, \psi_{n+1})$? Hint: nearly degenerate.

Finding the expansion coefficients¶

We now calculate the expansion coefficents as explained in the introduction and in the appendices.

Computing $\Psi(x,t)$¶

Now, everything is set to compute the wave function at some arbitrary time $t$ given the inital condition and potential. To do this, we create a function performing the calculation as explained in the introduction and in the appendices.

Finding a suitable time step - Ehrenfest's theorem¶

To find a suitable time step $\Delta t$, we will be using Ehrenfest's theorem. That is, the quantum mechanical expectation values obeys the classical equations of motion. For zero potential, (the expectation value of) the particle will thus have a velocity $$v = \frac{p(x)}{m} = \sqrt{\frac{2E_0}{m}}.$$ We will thus use $\Delta t \sim \sqrt{m/2E_0}\Delta x$.

Let us plot the result for some $t$'s!

Tunneling, reflection and transmission¶

There are many things one can learn from this simple exercise. For example, note that we have used an energy that is higher than the potential barrier, $E_0>V_0$. In classical mechanics we would expect total transmission, but from the plot above, we see that there is a probability for reflection! On the other hand, if $E_0<V_0$ we would classically expect total reflection, but there is some probability for transmission (test for yourself)! This is called tunneling. These concepts are explained in more detail in our notebook on One-Dimensional Wave Packet Propagation, and the different probabilities are explicitly calculated.

Note how the wave function has a high peak at the barrier. This is again due to reflection and transmission. In quantum mechanics we will have some reflection both when the potential is lowered and raised (check for yourself with a potential well!). The peak is thus due to constructive interference between different parts of the wave function being reflected repeatedly.

Exercises and further work¶

Investigate the problem further by yourself!

• What are the advantages and disadvantages using this method (opposed to the more direct method used in our notebook on One-Dimensional Wave Packet Propagation)?
• Compute numerically the transmission and reflection coefficient for different barrier widths and different barriers.
• Implement periodic boundary conditions. (Hint: Take a look at the matrix in the appendices and consider the boundary condition at the edges. We need to add two new non-zero matrix elements. These are located in the upper right and lower left corner. What are they? Note that we also need to use a sparse matrix or general eigenvalue solver, e.g. numpy.eigh)
• Explain why we have dispersion of the wave packet (it is spreading out).
• Calculate (you can make approximations if it is necessary) how long it takes for the electron to pass the barrier, reflect the right boundary, pass the barrier again and return to its initial position. Verify your calculations using the Python codes in this notebook.
• Generalize the method to two dimensions. (Hint: Use the same finite difference method as in the appendix on the two-dimensional Schrödinger equation. For simplicity, use $\Delta x = \Delta y = h$. To write the resulting approximation as a matrix, use the reindexing $i,j\to i + (j-1)N$. Treat carefully the boundaries! The easiest boundary condition is probably the Dirichlet boundary condition.)

Animation¶

Let us make an animation to visualize the propagating electron! It may also be instructive to calculate the probabilities for the particle to be in the different parts of the domain. When the particle has propagated through the barrier, the probability that the particle is on the right side of the barrier should be approximately equal to the transmission coefficient.