Are there simple ways to numerically solve the time-dependent Schrödinger equation?

Question

I would like to run some simple simulations of scattering of wavepackets off of simple potentials in one dimension.

Are there simple ways to numerically solve the one-dimensional TDSE for a single particle? I know that, in general, trying to use naïve approaches to integrate partial differential equations can quickly end in disaster. I am therefore looking for algorithms which

• are numerically stable,
• are simple to implement, or have easily accessible code-library implementations,
• run reasonably fast, and hopefully
• are relatively simple to understand.

I would also like to steer relatively clear of spectral methods, and particularly of methods which are little more than solving the time-independent Schrödinger equation as usual. However, I would be interested in pseudo-spectral methods which use B-splines or whatnot. If the method can take a time-dependent potential then that's definitely a bonus.

Of course, any such method will always have a number of disadvantages, so I would like to hear about those. When does it not work? What are common pitfalls? Which ways can it be pushed, and which ways can it not?

Answer 1

In the early 90s we were looking for a method to solve the TDSE fast enough to do animations in real time on a PC and came across a surprisingly simple, stable, explicit method described by PB Visscher in Computers in Physics: "A fast explicit algorithm for the time-dependent Schrödinger equation". Visscher notes that if you split the wavefunction into real and imaginary parts, $\psi=R+iI$, the SE becomes the system:

\begin{eqnarray}\frac{dR}{dt}&=&HI \\ \frac{dI}{dt}&=&-HR \\ H&=&-\frac{1}{2m}\nabla^2+V\end{eqnarray}

If you then compute $R$ and $I$ at staggered times ($R$ at $0,\Delta t,2\Delta t,...$ and $I$ at $0.5\Delta t, 1.5\Delta t,...)$, you get the discretization:

$$R(t+\frac{1}{2} \Delta t)=R(t-\frac{1}{2} \Delta t)+\Delta t HI(t)$$

$$I(t+\frac{1}{2} \Delta t)=I(t-\frac{1}{2} \Delta t)-\Delta t HR(t)$$

with $$\nabla^2\psi(r,t)=\frac{\psi(r+\Delta r,t)-2\psi(r,t)+\psi(r-\Delta r,t)}{\Delta r^2}$$ (standard three-point Laplacian).

This is explicit, very fast to compute, and second-order accurate in $\Delta t$.

Defining the probability density as

$$P(x,t)=R^2(x,t)+I(x,t+\frac{1}{2} \Delta t)I(x,t-\frac{1}{2} \Delta t)$$ at integer time steps and,

$$P(x,t)=R(x,t+\frac{1}{2} \Delta t)R(x,t-\frac{1}{2} \Delta t)+I^2(x,t)$$ at half-integer time steps

makes the algorithm unitary, thus conserving probability.

With enough code optimization, we were able to get very nice animations computed in real-time on 80486 machines. Students could "draw" any potential, choose a total energy, and watch the time-evolution of a gaussian packet.

Answer 2

The Schroedinger equation is effectively a reaction-diffusion equation $$ i\frac{\partial\psi}{\partial t}=-\nabla^2\psi+V\psi\tag{1} $$ (all constants are 1). When it comes to any partial differential equation, there's two ways to solve it:

1. Implicit method (adv: large time steps & unconditionally stable, disadv: requires matrix solver that can give bad data)
2. Explicit method (adv: easy to implement, disadv: requires small timesteps for stability)

For parabolic equations (linear in $t$ and 2nd order in $x$), the implicit method is often the better choice. The reason is the condition for stability for the explicit method requires $dt\propto dx^2$, which will be very small. You can avoid this issue by using the implicit method, which has no such limitation on the time-step (though in practice you don't normally make it insanely large because you can lose some of the physics). What I describe next is the Crank-Nicolson method, a common second order accurate (space & time) implicit scheme.

Starters

In order to computationally solve a PDE, you need to discretize it (make the variables fit onto a grid). The most straight-forward is a rectangular, Cartesian grid. Here on, $n$ represents the time index (and is always a super-script) and $j$ the position index (always a subscript). By employing a Taylor expansion for the position-dependent variable, Equation (1) becomes $$ i\frac{\psi^{n+1}_j-\psi_j}{dt} = -\frac12\left(\frac{\psi_{j+1}^{n+1}-2\psi_j^{n+1}+\psi_{j-1}^{n+1}}{dx^2}+\frac{\psi_{j+1}^{n}-2\psi_j^{n}+\psi_{j-1}^{n}}{dx^2}\right) \\ +\frac12\left(V_j\psi_j^{n+1} +V_j\psi_j^n\right) $$ Where we have assumed that $V=V(x)$. What happens next is a grouping of like spatial and temporal indices (you may want to double check the math): $$ \frac12\frac{dt}{dx^2}\psi_{j+1}^{n+1}+\left(i-\frac{dt}{dx^2}-\frac12V_j\right)\psi_j^{n+1}+\frac12\frac{dt}{dx^2}\psi_{j-1}^{n+1}=\\ i\psi_j^n-\frac12\frac{dt}{dx^2}\left(\psi_{j+1}^n-2\psi_j^n+\psi_{j-1}^n\right)+\frac12V_j\psi_j^n\tag{2} $$ This equation has the form $$ \left(\begin{array}{ccccc}A_0 & A_- & 0 & 0 &\cdots\\ A_+ & A_0 & A_- & 0 &\cdots\\ 0 & A_+ & A_0 & A_- &\cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots\end{array}\right)\left(\begin{array}{c}\psi_0^{n+1} \\ \psi_1^{n+1}\\ \vdots \\ \psi_{J-1}^{n+1}\end{array}\right)=\left(\begin{array}{c}\psi_0^{n} \\ \psi_1^{n}\\ \vdots \\ \psi_{J-1}^{n}\end{array}\right) $$ Which is called a tri-diagonal matrix and has a known solution (plus working examples, including one written by me!). The explicit method scratches out the entire left side (or should I say top line?) of Equation (2) except for the $i\psi_j^{n+1}$ term.

Issues

The biggest issue that I have found with implicit methods is that they are strongly dependent on the boundary conditions. If you have poorly defined/implemented boundary conditions, you can get spurious oscillations in your cells that can lead to bad results (see my SciComp post on a similar topic). This leads to actually having 1st order accuracy in space, rather than 2nd that your scheme ought to give.

Implicit methods are also supposedly difficult to parallelize, but I have only used them for 1D heat equations and not needed parallel support, so I can neither verify nor deny the claim.

I am also unsure how the complex nature of the wave function will affect the calculations. The work that I have done uses the Euler fluid dynamic equations, and are thus entirely real with non-negative magnitudes.

Time-dependent potential

If you have an analytic time-dependent potential (e.g. $V\propto \cos(\omega t)$), then you'd simply use the current time, $t$, for the $V_j$ on the RHS of (2) and the future time, $t+dt$, on the LHS. I don't believe that this would create any problems, but I have not tested this so I cannot verify or deny this aspect too.

Alternatives

There are some interesting alternatives to the Crank-Nicolson method as well. The first up is the so-called "super-time-stepping" method. In this explicit method, you take the time step ($dt\propto dx^2$) and use the roots of Chebyshev polynomials to get an optimized set of time-steps that quickly sum to $dt$ faster than doing $dt/N$ steps $N$ times (effectively you get $\Delta T=N^2dt$ so that each step $N$ advances you $Ndt$ in time). (I employ this method in my research because you have a well-defined "flux" from one cell to another that is used for merging data from one processor to another, using the Crank-Nicolson scheme I was unable to do this).

EDIT One thing to note is that this method is first order accurate in time, but if you use a Runge-Kutta 2 method in conjunction, it will give you a 2nd order accurate scheme in time.

The other is called an alternating-direction explicit. This method requires you to have known and well-defined boundary conditions. It then proceeds to solve the equation by using the boundary directly in the computation (no need to apply it after each step). What happens in this method is you solve the PDE twice, once in an upward sweep and once in a downward sweep. The upward sweep uses $$ \frac{\partial^2\psi}{\partial x^2}\approx\frac{\psi_{j-1}^{n+1}-\psi_j^{n+1}-\psi_j^n+\psi_{j+1}^n}{dx^2} $$ while the downward sweep uses $$ \frac{\partial^2\psi}{\partial x^2}\approx\frac{\psi_{j+1}^{n+1}-\psi_j^{n+1}-\psi_j^n+\psi_{j-1}^n}{dx^2} $$ for the diffusion equation while the other terms would remain the same. The time-step $n+1$ is then solved by averaging the two directional sweeps.

Answer 3

Kyle Kanos's answer looks to be very full, but I thought I'd add my own experience. The split-step Fourier method (SSFM) is extremely easy to get running and fiddle with; you can prototype it in a few lines of Mathematica and it is, extremely stable numerically. It involves imparting only unitary operators on your dataset, so it automatically conserves probability / power (the latter if you're solving Maxwell's equations with it, which is where my experience lies). For a one-dimensional Schrödinger equation (i.e. $x$ and $t$ variation only), it is extremely fast even as Mathematica code. And if you need to speed it up, you really only need a good FFT code in your target language (my experience lies with C++).

What you'd be doing is a disguised version of the Beam Propagation Method for optical propagation through a waveguide of varying cross section (analogous to time varying potentials), so it would be helpful to look this up too.

The way I look at the SSFM/BPM is as follows. Its grounding is the Trotter product formula of Lie theory:

$$\tag{1}\lim\limits_{m\to\infty}\left(\exp\left(\mathcal{D}\,\frac{t}{m}\right)\,\exp\left(\mathcal{V}\,\frac{t}{m}\right)\right)^m = \exp((\mathcal{D+V}) t)$$

which is sometimes called the operator splitting equation in this context. Your dataset is an $x-y$ or $x-y-z$ discretised grid of complex values representing $\psi(x,y,z)$ at a given time $t$. So you imagine this (you don't have to do this; I'm still talking conceptually) whopping grid written as an $N$-element column vector $\Psi$ (for a $1024\times1024$ grid we have $N=1024^2=1\,048\,576$) and then your Schrödinger equation is of the form:

$$\tag{2}\mathrm{d}_t \Psi = K\Psi = (\mathcal{D+V}(t)) \Psi$$

where $K = \mathcal{D+V}$ is an $N\times N$ skew-Hermitian matrix, an element of $\mathfrak{u}(N)$, and $\Psi$ is going to be mapped with increasing time by an element of the one parameter group $\exp(K\,t)$. (I've sucked the $i\hbar$ factor into the $K = \mathcal{D+V}$ on the RHS so I can more readily talk in Lie theoretic terms). Given the size of $N$, the operators' natural habitat $\mathfrak{U}(N)$ is a thoroughly colossal Lie group so PHEW! yes I am still talking in wholly theoretical terms!. Now, what does $\mathcal{D+V}$ look like? Still imagining for now, it could be thought of as a finite difference version of $i\,\hbar\,\nabla^2/(2\,m) - i\hbar^{-1}V_0 + i\hbar^{-1}(V_0-V(x,y,z,t_0))$, where $V_0$ is some convenient "mean" potential for the problem at hand.

We let:

$$\tag{3}\begin{array}{lcl}\mathcal{D} &=& i\frac{\hbar}{2\,m} \nabla^2 - i\hbar^{-1}V_0\\ \mathcal{V}&=&i\hbar^{-1}(V_0-V(x,y,z,t))\end{array}$$

Why I have split them up like this will become clear below.

The point about $\mathcal{D}$ is that it can be worked out analytically for a plane wave: it is a simple multiplication operator in momentum co-ordinates. So, to work out $\Psi\mapsto\exp(\Delta t\,\mathcal{D}) \Psi$, here are the first three steps of a SSFM/BPM cycle:

1. Impart FFT to dataset $\Psi$ to transform it into a set $\tilde{\Psi}$ of superposition weights of plane waves: now the grid co-ordinates have been changed from $x,\,y,\,z$ to $k_x,\,k_y,\,k_z$;
2. Impart $\tilde{\Psi}\mapsto\exp(\Delta t\,\mathcal{D}) \tilde{\Psi}$ by simply multiplying each point on the grid by $\exp(i\,\Delta t (V_0-k_x^2+k_y^2+k_z^2)/\hbar)$;

3. Impart inverse FFT to map our grid back to $\exp(\Delta t\,\mathcal{D}) \Psi$

   .Now we're back in position domain. This is the better domain to impart the operator $\mathcal{V}$ of course: here $\mathcal{V}$ is a simple multiplication operator. So here is your last step of your algorithmic cycle:

4. Impart the operator $\Psi\mapsto\exp(\Delta t\,\mathcal{V}) \Psi$ by simply multiplying each point on the grid by the phase factor $\exp(i\,\Delta t\,(V_0-V(x,y,z,t))/\hbar)$

....and then you begin your next $\Delta t$ step and cycle over and over. Clearly it is very easy to put time-varying potentials $V(x,y,z,t)$ into the code.

So you see you simply choose $\Delta t$ small enough that the Trotter formula (1) kicks in: you're simply approximating the action of the operator $\exp(\mathcal{D+V}\,\Delta t)\approx\exp(\mathcal{D}\,\Delta t)\,\exp(\mathcal{V}\,\Delta t)$ and you flit back and forth with your FFT between position and momentum co-ordinates, i.e. the domains where $\mathcal{V}$ and $\mathcal{D}$ are simple multiplication operators.

Notice that you are only ever imparting, even in the discretised world, unitary operators: FFTs and pure phase factors.

One point you do need to be careful of is that as your $\Delta t$ becomes small, you must make sure that the spatial grid spacing shrinks as well. Otherwise, suppose the spatial grid spacing is $\Delta x$. Then the physical meaning of the one discrete step is that the diffraction effects are travelling at a velocity $\Delta x/\Delta t$; when simulating Maxwell's equations and waveguides, you need to make sure that this velocity is much smaller than $c$. I daresay like limits apply to the Schrödinger equation: I don't have direct experience here but it does sound fun and maybe you could post your results sometime!

A second "experience" point with this kind of thing - I'd be almost willing to bet this is how you'll wind up following your ideas. We often have ideas that we want to do simple and quick and dirty simulations but it never quite works out that way! I'd begin with the SSFM as I've described above as it is very easy to get running and you'll quickly see whether or not its results are physical. Later on you can use your, say Mathematica SSFM code check the results of more sophisticated code you might end up building, say, a Crank Nicolson code along the lines of Kyle Kanos's answer.

---

Error Bounds

The Dynkin formula realisation of the Baker-Campbell-Hausdorff Theorem:

$$\exp(\mathcal{D}\Delta t)\,\exp(\mathcal{V})\Delta t) = \exp\left((\mathcal{D}+\mathcal{V})\Delta t + \frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,\Delta t^2 + \cdots\right)$$ converging for some $\Delta t>0$ shows that the method is accurate to second order and can show that:

$$\exp(\mathcal{D}\Delta t)\,\exp(\mathcal{V})\Delta t)\,\exp\left(-\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,\Delta t^2\right) = \exp\left((\mathcal{D}+\mathcal{V})\Delta t + \mathcal{O}(\Delta t^3)\right)$$

You can, in theory, therefore use the term $\exp(\mathcal{V})\Delta t)\,\exp\left(-\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,\Delta t^2\right)$ to estimate the error and set your $\Delta t$ accordingly. This is not as easy as it looks and in practice bounds end up being instead rough estimates of the error. The problem is that:

$$\frac{\Delta t^2}{2}[\mathcal{D},\,\mathcal{V}] = -\frac{i\,\Delta t^2}{2\,m}\,\left(\partial_x^2 V(x,\,t) + 2 \partial_x V(x,\,t)\,\partial_x\right)$$

and there are no readily transformed to co-ordinates wherein $[\mathcal{D},\,\mathcal{V}]$ is a simple multiplication operator. So you have to be content with $\exp\left(-\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,\Delta t^2\right) \approx e^{-i\,\varphi\,\Delta t^2}\left(\mathrm{id} -\left(\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,-i\,\varphi(t)\right)\,\Delta t^2\right)$ and use this to estimate your error, by working out $\left(\mathrm{id} -\left(\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,-i\,\varphi(t)\right)\,\Delta t^2\right) \,\psi$ for your currently evolving solution $\psi(x,\,t)$ and using this to set your $\Delta t$ on-the-fly after each cycle of the algorithm. You can of course make these ideas the basis for an adaptive stepsize controller for your simulation. Here $\varphi$ is a global phase pulled out of the dataset to minimise the norm of $\left(\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,-i\,\varphi(t)\right)\,\Delta t^2$; you can of course often throw such a global phase out: depending on what you're doing with the simulation results often we're not bothered by a constant phase global $\exp\left(\int \varphi\,\mathrm{d}t\right)$.

A relevant paper about errors in the SSFM/BPM is:

Lars Thylén. "The Beam Propagation Method: An Analysis of its Applicability", Optical and Quantum Electronics 15 (1983) pp433-439.

Lars Thylén thinks about the errors in non-Lie theoretic terms (Lie groups are my bent, so I like to look for interpretations of them) but his ideas are essentially the same as the above.

Answer 4

I can recommend using the finite-difference time-domain (FDTD) method. I even wrote a tutorial some time back that should answer most of your questions:

J. R. Nagel, "A review and application of the finite-difference time-domain algorithm applied to the Schrödinger equation," ACES Journal, Vol. 24, No. 1, February 2009

I have some Matlab codes that run nicely for 1D systems. If you have experience with FDTD doing electromagnetics, it works great for quantum mechanics as well. I can post my codes if you're interested.

Basically, it just operates on the wavefunctions directly by splitting the derivatives up into finite differences. It is kind of similar to the Crank-Nicholson scheme, but not exactly. If you are familiar with FDTD from electromagnetic wave theory, then FDTD will be very intuitive when solving the Schrodinger equation.

Answer 5

The most straightforward finite difference method is fast and easy to understand, but is not unitary in time - so probability is not conserved. Crank-Nicholson-Crout averages the forward and backward finite difference methods to produce a hybrid implicit/explicit method that is still pretty easy to understand and to implement and is unitary in time. This site explains the method well, provides pseudocode, and gives the relevant properties:

http://www.physics.utah.edu/~detar/phycs6730/handouts/crank_nicholson/crank_nicholson/ Note: There is a - sign missing from the LHS of equation one of this link, which propagates throughout the page.

Where does the nonunitarity come from?

In a nut shell, solving the TDSE comes down to figuring out how to deal with

$| \psi(x,t)\rangle = e^{-iHt}|\psi(x,0)\rangle$

which contains a differential operator in an exponential.

Applying a forward finite difference turns the differential operator into a tridiagonal matrix (converting the Reals to a grid) and the exponential into the first two terms of its Taylor series

$e^{-iHt}\approx1-iHt $

This discretization and linearization is what gives rise to the nonunitarity. (You can show that the tridiagonal matrix is not unitary by direct computation.) Combining the forward finite difference with the backward finite difference produces the approximation

$e^{-iHt}\approx \frac {1-\frac{1}{2} iHt} {1+\frac{1}{2} iHt} $

which, kindly, happens to be unitary (again you can show it by direct computation).

Answer 6

A few answers and comments here conflate confusingly the TDSE with a wave equation; perhaps a semantics issue, to some extent. The TDSE is the quantized version of the classical non-relativistic hamiltonian $$H=\frac{p^2}{2m} + V(x)= E.$$ With the rules $$p \rightarrow i\hbar\partial_x,\ \ E\rightarrow i\hbar\partial_t, \ \ x\rightarrow x,$$ (as discussed in chapter 1 of d'Espagnat, Conceptual foundations of quantum mechanics, https://philpapers.org/rec/ESPCFO), it therefore reads $$\left[-\frac{\hbar^2}{2m}\partial_{xx} + V(x)\right]\psi = i\hbar\partial_t\psi,$$ so it is clearly a diffusion-like equation. If one used the relativistic energy, which contains a E$^2$ term, then a wave-like equation such as $$ \partial_{xx}\psi = \partial_{tt}\psi +\dots $$ would obtain (for V=0 for simplicity), such as the Pauli or Klein-Gordon equations. But that is, of course, a completely different matter.

Now, back to the TDSE, the obvious method is Crank-Nicolson as has been mentioned, because it is a small-time expansion that conserves the unitarity of the evolution (FTCS, e.g., does not). For the 1-space-D case, it can be treated as a matrix iteration, reading $$ \vec{\psi}^{n+1}=({\cal I} + \frac{i\tau}{2\hbar}\tilde{H})^{-1} ({\cal I} - \frac{i\tau}{2\hbar}\tilde{H})\vec{\psi}^{n} $$ with ${\cal I}$ the identity matrix and $$ H_{jk}=(\tilde{H})_{jk}=\frac{-\hbar^{2}}{2m}\left[\frac{\delta_{j+1,k}+\delta_{j-1,k}-2\delta_{jk}}{h^{2}}\right]+V_{j}\delta_{jk}. $$ (Details e.g. in Numerical methods for physics, http://algarcia.org/nummeth/nummeth.html, by A. L. Garcia). As most clearly seen in periodic boundary conditions, a space-localized $\psi$ spreads out in time: this is expected, because the initial localized $\psi$ is not an eigenstate of the stationary Schroedinger equation, but a superposition thereof. The (classical massive free particle) eigenstate with fixed momentum (for the non-relativistic kinetic operator) is simply $\psi_s=e^{ikx}/\sqrt{L}$, i.e. fully delocalized as per Heisenberg principle, with constant probability density 1/L everywhere (note that I am avoiding normalization issues with continuum states by having my particle live on a finite, periodically repeated line). Using C-N, the norm $$\int |\psi|^2 dx$$ is conserved, thanks to unitarity (this is not the case in other schemes, such as FTCS, e.g.). Incidentally, notice that starting from an energy espression such as $$cp=E$$ with $c$ fixed, you'd get $$ic\hbar\partial_x=i\hbar\partial_t$$ i.e. the advection equation, which has no dispersion (if integrated properly with Lax-Wendroff methods), and your wavepacket will not spread in time in that case. The quantum analog is the massless-particle Dirac equation.