Simulating Wavepacket Evolution with Crank-Nicolson

Hello! Welcome back to our module on numerical methods.

In our last lesson, you successfully implemented a finite difference solver for the time-independent Schrödinger equation. By turning the TISE into a matrix eigenvalue problem, you were able to find the stationary states and energy levels for any 1D potential, including the double-well. This gave us a static picture of the allowed quantum states.

Today, we move from the static to the dynamic. Our goal is to answer the question: if we start a quantum system in some initial state (that is not a stationary state), how does it evolve in time? To do this, we will tackle the full time-dependent Schrödinger equation (TDSE).

Your learning outcome for this lesson is to:

• Implement the Crank-Nicolson method to simulate wavepacket evolution for the 1D time-dependent Schrödinger equation.

This powerful technique is a workhorse in computational physics for simulating quantum dynamics, making it an essential addition to your numerical toolkit, especially given your interest in modeling new systems.

1. The Challenge of Time Evolution

The evolution of a quantum state $\psi (x,t)$ is governed by the 1D time-dependent Schrödinger equation:

$$i\hslash \frac{\partial \psi (x,t)}{\partial t}=\hat{H}\psi (x,t)$$

where $\hat{H}=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}+V(x)$ is the same Hamiltonian operator you discretized in the last lesson.

Our task is to step the wavefunction forward in time. A simple approach might be to use a forward difference for the time derivative (the "Forward Time, Centered Space" or FTCS method):

$$i\hslash \frac{\psi (x,t+\Delta t)-\psi (x,t)}{\Delta t}\approx \hat{H}\psi (x,t)$$

Rearranging this gives an explicit update rule:

$$\psi (x,t+\Delta t)\approx \left(I-\frac{i\Delta t}{\hslash }\hat{H}\right)\psi (x,t)$$

While simple, this method has a major flaw: it's not unitary. The operator $(I-i\Delta t\hat{H}/\hslash )$ does not preserve the norm of the wavefunction, meaning the total probability $\int |\psi {|}^{2}dx$ will not be conserved. Over many time steps, this numerical error accumulates, leading to an unstable solution that can either vanish or blow up.

We need a more stable method. The Crank-Nicolson method solves this by being smarter about how it approximates the derivatives.

2. The Crank-Nicolson Scheme

The core idea of the Crank-Nicolson method is to be more balanced. Instead of evaluating the right-hand side (the $\hat{H}\psi$ term) only at the current time $t$, it uses an average of the values at the current time $t$ and the next time step $t+\Delta t$.

The time derivative is still approximated with a central difference, but it's centered at $t+\Delta t/2$:

$$\frac{\partial \psi }{\partial t}{|}_{t+\Delta t/2}\approx \frac{\psi (x,t+\Delta t)-\psi (x,t)}{\Delta t}$$

The Hamiltonian part is also evaluated at this midpoint by averaging:

$$\hat{H}\psi {|}_{t+\Delta t/2}\approx \frac{1}{2}\left(\hat{H}\psi (x,t+\Delta t)+\hat{H}\psi (x,t)\right)$$

The following article provides a clear derivation of the method from these starting points.

As you saw in the reading, putting these approximations into the TDSE and rearranging terms to separate the unknown future state (${\psi }^{n+1}$) from the known current state (${\psi }^n$) gives the central equation:

$$\left(I+\frac{i\Delta t}{2\hslash }\hat{H}\right){\psi }^{n+1}=\left(I-\frac{i\Delta t}{2\hslash }\hat{H}\right){\psi }^n$$

This is an implicit method because the unknown ${\psi }^{n+1}$ appears on both sides of the original discretized equation. The final form is a linear system $A{\psi }^{n+1}=B{\psi }^n$, where we must solve for ${\psi }^{n+1}$ at each step.

The key advantage here is that the operator that evolves the state, $U_{CN}=A^{-1}B$, is unitary. This guarantees that the norm of the wavefunction is conserved at each step, making the algorithm unconditionally stable.

This idea of averaging an explicit and implicit step is very general. The following video explains it in the context of the diffusion equation, which is mathematically analogous to the Schrödinger equation in imaginary time.

The structure of the update rule can be visualized with a "stencil" that shows which grid points are connected by the equation.

3. Implementation: From Matrix Equation to Algorithm

We now have a matrix equation to solve at each time step: $A{\psi }^{n+1}=B{\psi }^n$.

From the last lesson, you know that the discretized Hamiltonian $\hat{H}$ is a tridiagonal matrix. Because the matrices $A$ and $B$ are defined as $I\pm c\hat{H}$, they are also tridiagonal.

Our task is to solve this tridiagonal system. A direct inversion of matrix $A$ would be computationally expensive ($\mathcal{O}(N^3)$). Fortunately, for tridiagonal systems, there is a much faster method called the Thomas Algorithm (or the Tridiagonal Matrix Algorithm, TDMA), which is a specialized form of Gaussian elimination with a cost of only $\mathcal{O}(N)$.

The article you read earlier continues with a perfect explanation of how to construct the matrices $A$ and $B$ and how the Thomas Algorithm is used to solve the system.

Now you have all the theoretical and algorithmic pieces. It's time to put them together in code.

4. A Python Implementation

Implementing the Crank-Nicolson method with the Thomas algorithm from scratch is a fantastic exercise. The following section of the matterwavex.com article provides a complete, well-commented Python script that does exactly this.

The simulation in the article demonstrates a classic quantum phenomenon: a wavepacket impinging on a potential barrier that is higher than the particle's kinetic energy. You will see part of the wavepacket reflect and part of it tunnel through, a purely quantum effect.

The video below uses a slightly different numerical method (a simple explicit one) but provides an excellent visualization of what the wavepacket dynamics should look like. It also offers a great alternative way to check your results: by expanding the initial state in the energy eigenbasis (which you learned to find in the last lesson!).

Test your understanding! Apply your knowledge.

Now, let's connect this lesson with the previous one. Your task is to simulate quantum tunneling in the double-well potential.

Problem: Simulate the time evolution of a Gaussian wavepacket in the double-well potential you studied in the last lesson's exercise:
$V(x)=V_0{\left({\left(\frac{x}{a}\right)}^2-1\right)}^2$

Instructions:

1. Adapt the Python code from the matterwavex.com article.
2. Replace the square barrier potential with the double-well potential. Use the same parameters as last time (e.g., $a=2$, $V_0=10$ in atomic units).
3. Modify the initial wavepacket so that it starts centered in one of the wells (e.g., at $x=-a$), with zero initial momentum ($k_0=0$). Make the packet narrow enough to fit mostly within that well.
4. Run the simulation. You should observe the wavepacket tunneling back and forth between the two wells. This is a beautiful visualization of the superposition of the symmetric ground state and the anti-symmetric first excited state.
5. Bonus: Try to measure the time it takes for the probability to move from one well to the other and back again (the tunneling period, $T_{\text{tunnel}}$). In the previous lesson, you found the energies $E_0$ and $E_1$ for the two lowest states. How does your measured $T_{\text{tunnel}}$ compare to the theoretical prediction $T=h/(E_1-E_0)$?

Show hints on what to expect

You should see the probability density "slosh" from the left well to the right well and then back again. The rate of this tunneling is determined by the energy splitting $\Delta E=E_1-E_0$. A smaller energy splitting (which you found occurs for higher barriers) corresponds to a longer tunneling period, as the particle is more "trapped" in each well. This exercise beautifully connects the static eigenstate picture (energy splitting) with the dynamic picture (tunneling time).

Conclusion

In this lesson, you've added a crucial method for simulating quantum dynamics to your computational arsenal. By moving beyond simple explicit methods, you can now reliably model how quantum systems change over time.

Key Takeaways:

• The Time-Dependent Schrödinger Equation can be solved numerically to simulate quantum dynamics.
• Simple explicit methods like FTCS are often unstable because they don't conserve probability.
• The Crank-Nicolson method is an implicit, unconditionally stable method that reformulates the TDSE into a sequence of linear systems, $A{\psi }^{n+1}=B{\psi }^n$.
• Because the discretized Hamiltonian is tridiagonal, the resulting system can be solved very efficiently in $\mathcal{O}(N)$ time using the Thomas Algorithm.

Preview of the Next Lesson:

So far, we have focused on single-particle quantum mechanics. However, your primary interest lies in many-body systems. In the next lesson, we will make the leap into this domain. You will apply exact diagonalization to calculate the ground state and low-lying excitations of a small quantum many-body system, such as a Heisenberg spin chain. This technique is a cornerstone of computational condensed matter physics.