# 11.3: The Split-Step Fourier Method

- Page ID
- 34866

As an example of the usefulness of the DFT, let us discuss a DFT-based method for performing numerical integration of a *partial* differential equation, known as the **split-step Fourier method**. Here, the method will be presented in the context of the time-dependent Schrödinger equation in 1D space:

\[i \frac{d\psi(x,t)}{dt} = \left[-\frac{1}{2}\frac{\partial^2}{\partial x^2} + V(x,t)\right]\psi(x,t)\]

We have taken \(\hbar = m = 1\) for simplicity. At each time, the wavefunction is a continuous function of \(x\). Let us truncate and discretize this spatial coordinate, by defining a computational domain of length \(L\) containing \(N\) discretization points:

\[\psi_n(t) = \psi(x_n,t), \;\;\;\textrm{where}\;\; x_n = -\frac{L}{2} + n \Delta x, \;\; \Delta x = \frac{L}{N}.\]

Thus, the wavefunction at each time is represented by a complex vector, which we call a "state vector":

\[\vec{\psi}(t) = \begin{pmatrix} \psi_0(t) \\ \psi_1(t) \\ \vdots \\ \psi_{N-1}(t)\end{pmatrix}.\]

Given an initial state vector \(\vec{\psi}(t_a)\), the problem is to compute \(\vec{\psi}(t_b)\) at a later time \(t_{b}\). Note that this differs from previously-studied numerical ODE problems in one important respect: evolving in time involves taking second-order *spatial* derivatives. We won't go into the details, but it turns out that standard methods for time-stepping and discretizing space don't work very well here, because the errors from time-stepping and spatial discretization interact badly with one another. The split-step Fourier method provides a better way to solve the problem.

**11.3.1 Factorizing the Time-Evolution Operator**

The split-step Fourier method is based on the concept of the *time-evolution operator*. Given a wavefunction \(\psi(x,t_j)\), the wavefunction after a small time step \(\tau\) is

\[\psi(x,t_j + \tau) = U(t_j + \tau|t_j)\, \psi(x,t_j), \quad\mathrm{where}\;\;\; U(t_j + \tau|t_j) \approx \exp\left\{-i \tau \left[-\frac{1}{2}\frac{\partial^2}{\partial x^2} + V(x,t_j+\tau/2)\right]\right\}\;.\]

Here, the \(\exp(\cdots)\) refers to the exponential of an operator (one involving spatial derivatives). We call \(U(t_b|t_a)\) the time-evolution operator, which evolves the system from time \(t_{a}\) to \(t_{b}\). The exponential of any operator \(\mathcal{A}\) is defined as the infinite series

\[\exp(\mathcal{A}) = I + \mathcal{A} + \frac{1}{2} \mathcal{A}^2 + \frac{1}{6} \mathcal{A}^3 + \cdots\]

In this case, the exponential contains the Hamiltonian, which consists of a kinetic energy term and a potential energy term that do not generally commute. Due to this non-commutivity, the exponential cannot be simplified by factorization:

\[\exp\left\{-i \tau \left[-\frac{1}{2}\frac{\partial^2}{\partial x^2} + V(x,t_j+\tau/2)\right]\right\} \;\; \ne\;\;\exp\Bigg\{\frac{i \tau}{2} \left[\frac{\partial^2}{\partial x^2}\right]\Bigg\} \;\exp\Bigg\{-i \tau V(x,t_j+\tau/2)\Bigg\}.\]

However, we can obtain an approximate factorization by making use of the series definition of the exponential of an operator. One can show that

\[\exp(\mathcal{A} + \mathcal{B}) = \exp(\mathcal{A}/2) \exp(\mathcal{B}) \exp(\mathcal{A}/2) + O\Big((\mathcal{A},\mathcal{B})^3\Big),\]

which is a variant of an important formula known as the Baker–Campbell–Hausdorff formula. On the right-hand side, note that \(\exp(\mathcal{B})\) is sandwiched "symmetrically" between two copies of \(\exp(\mathcal{A}/2)\). This symmetric arrangement reduces the approximation error to third order, by the cancellation of lower-order errors (in a manner similar to the mid-point formula for the discretized derivative). Applying this factorization to the time-evolution operator gives

\[\boxed{U(t_j + \tau|t_j) \;\; \approx\;\; U_\mathcal{K} \,\cdot\, U_\mathcal{V}(t_j+\tau/2) \,\cdot\, U_\mathcal{K}, \quad \mathrm{where}\;\; \begin{cases} \;U_\mathcal{K} &= \exp\left[\frac{i\tau}{4} \frac{d^2}{dx^2}\right] \\ \; U_\mathcal{V}(t) &= \exp\big[-i\tau\, V(x,\,t)\big]\end{cases}}\]

In other words, the time-evolution operator decomposes into three pieces. That's why we call this a "split-step" algorithm: each time step from \(t_{j}\) to \(t_{j}+\tau\) consists of applying a kinetic step, then applying a potential step, then applying another kinetic step, in sequence. As previously noted, we'll be working with state vectors (complex \(N\)-component vectors), defined through spatial discretization of the wavefunction. So we need to figure out how the above stepping operators act on these state vectors:

\[U_\mathcal{K}\; \vec{\psi} \;= \;\;??, \qquad U_\mathcal{V}\, \vec{\psi} \;=\;\; ??\]

The potential stepping operator is simple to deal with. Since the state vector represents the wavefunction at different points in space, the potential operator is represented by a diagonal matrix, and its exponential is also diagonal:

\[U_\mathcal{V} \begin{pmatrix}\psi_0\\ \vdots \\ \psi_{N-1}\end{pmatrix} \;=\; \begin{pmatrix}\exp\left[-i\tau V(x_0,t+\frac{\tau}{2})\right] & & \\ &\ddots&\\ && \exp\left[-i\tau V(x_{N-1},t+\frac{\tau}{2})\right] \end{pmatrix} \begin{pmatrix}\psi_0\\ \vdots \\ \psi_{N-1}\end{pmatrix} \]

## 11.3.2 Kinetic Step

The kinetic stepping operator, \(U_\mathcal{K}\), is less obvious. It contains spatial derivatives and is thus *not* diagonal in the current basis. The key thing to realize, however, is that this operator is diagonal in *wavenumber space*. Let us return to the continuous wavefunction, and write its Fourier representation:

\[\psi(x) = \int_{-\infty}^\infty \frac{dk}{2\pi} \; e^{ikx}\; \Psi_k.\]

Then

\[U_\mathcal{K} \; \psi(x) = \int_{-\infty}^\infty \frac{dk}{2\pi} \; \exp\left(-\frac{i\tau}{4}k^2\right) \; e^{ikx}\; \Psi_k.\]

Let us discretize space in steps of \(\Delta x\), as discussed earlier, and also discretize the Fourier integrals by steps of \(\Delta k\):

\[\begin{align}x_n &= -\frac{L}{2} + n \Delta x, \\k_n &= -\frac{K}{2} + n\Delta k\end{align}\]

The values of \(K\) and \(\Delta k\) will be chosen shortly. The discretized integrals become

\[\begin{align}\psi_m &\approx \sum_{n=0}^{N-1} \frac{\Delta k}{2\pi} \; e^{ik_n x_m}\; \Psi_{k_n}, \\\left(U_\mathcal{K} \psi\right)_m &\approx \sum_{n=0}^{N-1} \frac{\Delta k}{2\pi} \; e^{-\frac{i\tau}{4}k_n^2} \;e^{ik_n x_m}\; \Psi_{k_n}. \end{align}\]

Let us now choose the \(k\)-space discretization parameters to be

\[\Delta k = \frac{2\pi}{N\Delta x} = \frac{2\pi}{L} , \quad K = N\Delta k = \frac{2\pi}{\Delta x} \;\;\;\Rightarrow \;\;\; k_n = -\frac{\pi}{\Delta x} + \frac{2\pi n}{N\Delta x}.\]

With this choice, we can show with a bit of algebra that the integral for \(\psi_{m}\) reduces to an IDFT:

\[\begin{align}\psi_m &= (-1)^m \frac{1}{N} \sum_{n=0}^{N-1} \left(\frac{1}{\Delta x} e^{iN\pi/2}\;e^{-in\pi}\; \Psi_{k_n}\right) e^{\frac{2\pi i mn}{N}}\\ &= (-1)^m\; \mathrm{IDFT}\left\{\frac{1}{\Delta x} e^{iN\pi/2}\;(-1)^n \;\Psi_{k_n}\right\}_m \\ \Rightarrow \Psi_{k_n} &= \Delta x\, e^{-iN\pi/2}\; (-1)^n \; \;\mathrm{DFT}\Big\{(-1)^m\; \psi_m \Big\}_n \end{align}\]

Likewise,

\[(U_\mathcal{K}\; \psi)_m = (-1)^m \; \mathrm{IDFT}\left\{\frac{1}{\Delta x} e^{iN\pi/2}\;(-1)^n e^{-\frac{i\tau}{4}k_n^2}\; \Psi_{k_n}\right\}_m\]

Putting these results together, we get

\[\begin{align}(U_\mathcal{K}\; \psi)_m = (-1)^m \;\mathrm{IDFT}\Bigg\{ e^{-\frac{i\tau}{4}k_n^2}\; \mathrm{DFT}\Big\{(-1)^{p}\, \psi_{p}\Big\}_n \Bigg\}_m \\ \mathrm{where}\;\; k_n = - \frac{\pi N}{L} + \frac{2\pi n}{L} \end{align}\]

Hence, the \(U_{\mathcal{K}}\) kinetic stepping operator can be implemented by taking a DFT, multiplying the resulting vector elements by \(\exp(-i\tau k_n^2/4)\) phase factors, and taking an IDFT. The runtime of the stepping process is \(O(N\log(N))\). The \(m\), \(n\), and \(p\) indices all run over the range \([0, 1, \cdots, N-1]\), consistent with the standard definition of the DFT and IDFT.