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

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?

• – Kyle Kanos Feb 7 '14 at 21:37
• @EmilioPisanty I've added a discussion of errors to my writeup of the SSFM: I notice (after writing my answer, sorry) you're not keen on spectral methods, but just in case ... – WetSavannaAnimal aka Rod Vance Feb 12 '14 at 1:52
• Cleaned up thread; removing the discussion of topicality from Physics. – Geoff Oxberry Feb 23 '14 at 4:09
• I recommend starting here, though there is a lot of more recent work: A comparison of different propagation schemes for the time dependent Schrödinger equation. – David Ketcheson Feb 23 '14 at 5:16
• @GeoffOxberry could you make available a screenshot of those comments? – Emilio Pisanty Feb 25 '14 at 18:39

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.

• Great answer, only complaint is that you beat me to it! – Kyle Feb 8 '14 at 18:20
• @ChrisWhite: I was thinking about how that could be done earlier this morning and the only thing I came up with is be doing it once for $\mathbb R$ and once for $\mathbb I$. I'll take a look at that paper (and more importantly the free code they give out) and see how they suggest to do it. – Kyle Kanos Feb 8 '14 at 19:00
• @ChrisWhite Maybe the sneakiness is for calculating eigenfunctions, which I have seen calculated by imaginary timestepping: you arrange the step direction so that the lowest energy eigenfunction has the least negative value of $-h\,\nu$ and thus the slowest decay. On iterating on a random input, very swiftly only the shape of the lowest energy eigenfunction is left. Then you subtract this from the random input and do the process again: now the next lowest energy eigenfunction is the dominant one. And so on. It sounds a bit dodgy (especially getting higher $-h\,\nu$ eigenfuncs) but it works! – WetSavannaAnimal aka Rod Vance Feb 9 '14 at 12:39
• @DavidKetcheson: A reaction-diffusion equation takes the form $\partial_tu=D\partial_x^2u+R(u)$. In the case of the Schrodinger equation, $R(u)=Vu$; may I then ask how that it's not a RD-type equation? And, curiously, the Schrodinger equation actually appears in the reaction-diffusion wiki article I referenced. This equivocation I made also appears in many published journals and texts (go ahead and search it). Perhaps it would have been better for me to advise using standard libraries (as is the common MO here), such as PETSc, deal.ii, or pyCLAW? – Kyle Kanos Feb 23 '14 at 18:09
• @KyleKanos: Your post is good. In fact, in the article posted by DavidKetcheson, Crank-Nicolson is advocated by the first reference. The comparison to reaction-diffusion is fine; as you note, the comparison to reaction-diffusion does appear in many published sources. I think DavidKetcheson was looking for something like "dispersive wave equation" mentioned earlier. – Geoff Oxberry Feb 23 '14 at 20:44

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.

• That is a very neat trick for solving the real & imaginary components! Note also that you can get large, centered equations by using $$...$$. I've taken the liberty of doing this for you, I hope you don't mind! – Kyle Kanos Feb 10 '14 at 18:42
• We were delighted to find the algorithm - it was easy to program and ran fast. The hardest part was getting the initial conditions right, R at t=0 and I at 0.5dt... I don't mind the edit, I was happy to get equations at all. – Wally Feb 10 '14 at 19:05
• @user40172 We were doing the same thing for waveguides at the about the same time, and we settled on the BPM described in my answer. The reason was that at the time we could run the FFTs separately from the main CPU using a DSP board. We thought we were oh so clever, but I must say coming up with essentially a hardware solution to a software problem looks pretty naff in 2014 though! The latest version of Visual Studio C++ automatically vectorises code over CPUs and it does a beautiful job with the FFT. – WetSavannaAnimal aka Rod Vance Feb 11 '14 at 0:37
• @user40172 How did you get the initial conditions at $0.5dt$ finally? Just propagating solution to that time using another method? – Ruslan Feb 11 '14 at 18:56
• @Rusian Since we were doing scattering, we used a standard free-particle Gaussian wave packet but made sure to start it "far enough" away from any region where the potential was non-zero. See, for example: demonstrations.wolfram.com/EvolutionOfAGaussianWavePacket – Wally Feb 11 '14 at 19:20

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 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.

• Rod, you are probably aware that you can do better if you use the so-called split-operator approximation, where $\exp[\Delta t({\cal D} + {\cal V})] \approx \exp[\Delta t {\cal V}/2] \exp[\Delta t {\cal D}] \exp[\Delta t {\cal V}/2]$. In fact you can do some further splitting to carry the error to higher $\Delta t$ powers. See for instance Bandrauk and Shen, Chem. Phys. Lett. 176, 428 (1991). Obviously your kinetic term cannot depend on the coordinates, that is, it doesn't work nicely in curvilinear coordinates. – perplexity Feb 17 '14 at 11:01
• Otherwise, this split-operator thing coupled to the FFT evaluation of the kinetic energy operator is one of the standard procedures to solve the TDSE on a grid-based representation in Molecular Physics. – perplexity Feb 17 '14 at 11:05
• @perplexity Many thanks. It's good to know what different fields use. The 1991 date on your reference is interesting: I was always pretty sure the split step idea came out of waveguide simulation in the late 1970s - so maybe I'm wrong. – WetSavannaAnimal aka Rod Vance Feb 17 '14 at 12:12
• You are not wrong at all. That was the inspiration indeed. The first work translating these ideas to QM that I am aware of is Feit, Fleck and Steiger, J. Comput. Phys. 47, 412 (1982) where, if I recall correctly, they essentially use the same tricks with the advantage that the operator here is unitary by construction (unlike in classical waves). The FFT-grid based approach to these type of simulations was first proposed by Ronnie Kosloff, I believe. He has a very nice review about this subject on his web page. – perplexity Feb 17 '14 at 16:27
• Another good reference in my field is David Tannor's book on Quantum Mechanics: A time-dependent perspective. Cheers. – perplexity Feb 17 '14 at 16:28

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.

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).

• Thanks for the quick response. Could you provide more details on both those methods? How do they work, and why? Where does the nonunitarity come from? – Emilio Pisanty Feb 7 '14 at 22:51
• I would be happy to provide more detail, but to avoid missing my target audience, it would be useful to know how much education and experience you've had in each of the following background fields: Calculus, Differential Equations, Linear Algebra, Quantum Mechanics, and Numerical Methods (specifically Finite Difference Methods). – Wally Feb 8 '14 at 14:26
• Please assume as much as you need from standard physics and math (though references to the more complicated parts would probably help). My numerical methods are a bit rusty, though. – Emilio Pisanty Feb 8 '14 at 16:10
• Are there any differences between this and Kyle Kanos's answer? I mean, it's not obvious how to implement your last equation - as you've written it involves inverting a full operator - are you simply saying that the CN method is simply, through the solution of its tridiagonal equation, working out $(1+\frac{i}{2}\,H\,t)^{-1}\,(1+\frac{i}{2}\,H\,t)\,\psi$? Or is there a subtlety that I've missed? Actually you last equation is a good rendering insofar that it makes unitarity explicit for CN, a fact which is unclear in many descriptions of CN. – WetSavannaAnimal aka Rod Vance Feb 9 '14 at 12:17
• No, it is the same algorithm as given by Kyle Kanos. I just wrote it this way to give a different way of looking at it. I hoped easier to conceptualize - whereas his is easier to implement. Yes, you are ultimately just solving a tridiagonal equation. There was an old (1967) paper in AJP that I couldn't find earlier that describes it very well: ergodic.ugr.es/cphys/lecciones/SCHROEDINGER/ajp.pdf They used CN to produce 8mm film loops of gaussian wave packets scattering off various potentials. You can still find those film loops in many university physics demo libraries. – Wally Feb 9 '14 at 20:48

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.