Notes¶

Split-Step-Fourier method¶

Defining a time-propagator \(U(t,t_0)\) with its action on the wave function \(\Psi(t) = U(t,t_0)\Psi(t_0)\) and insert it into the Schrödinger equation, we obtain

\[i\hbar \frac{\mathrm{d}}{\mathrm{d}t}U(t,t_0)\Psi(t) = H(t)U(t,t_0)\Psi(t).\]

A solution to this is the well known Dyson series

\[U(t,t_{0})=1-\frac{i}{\hbar}\int _{{t_{0}}}^{{t}}{\mathrm{d}t_{1}H(t_{1})}+\bigg(-\frac{i}{\hbar}\bigg)^{2}\int _{{t_{0}}}^{t}{\mathrm{d}t_{1}\int _{{t_{0}}}^{{t_{1}}}{\mathrm{d}t_{2}H(t_{1})H(t_{2})}}+\cdots\]

which can be simplified using the time-sorting operator \(\hat{T}\):

\[U(t,t_0) = \hat{T}\exp\Big[-\frac{i}{\hbar}\int_{t_0}^tH(t')\mathrm{d}t' \Big].\]

The naive intuition of the time propagator to be \(\exp[-\frac{i}{\hbar}\int_{t_0}^tH(t')\mathrm{d}t']\), neglecting the time-sorting, is generally not true. If this naive ansatz is inserted into the Schrödinger equation one might expect the equality to hold, however, the time derivative of an exponentiated operator is only well defined if it is diagonal. Yet, generally a Hamiltonian does not commute with itself at different times \([H(t_1),H(t_2)] \neq 0\), hence they do not have the same eigenbasis and time-sorting has to be applied. Nonetheless, as an approximation for small time steps \(\Delta t\), we will assume that \([H(t),H(t+\Delta t)] = 0\) holds and the naive ansatz for the time propagator can be used. (That the commutator relation \([H(t),H(t+\Delta t)] = 0\) leads to \(U(t,t_0) \approx\exp[-\frac{i}{\hbar}\int_{t_0}^tH(t_1)\mathrm{d}t_1]\) can be even better understood using the Magnus expansion of the time propagator:

\[U(t,t_0) =1 -\frac{i}{\hbar}\int_{t_0}^tH(t_1)\mathrm{d}t_1+ \frac{1}{2} \bigg(-\frac{i}{\hbar}\bigg)^2 \int_{t_0}^t \mathrm{d}t_1 \int_{t_0}^{t_1} \mathrm{d}t_2 [H(t_1),H(t_2)]+\dots\]

The Dyson series is described, because it is the common way of writing down \(U(t,t_0)\).)

Hence, we define the time evolution of a state to be approximately

\[\psi(\vec{r},t+\Delta t) \approx e^{-\frac{i}{\hbar} H(\vec{r},t) \Delta t } \psi(\vec{r},t).\]

The Hamiltonian consists of a kinetic \(T\) and potential \(V\) part. As it is intended to evolve the state vector incrementally over time steps of length \(\Delta t\), it is necessary to calculate the exponential function of a sum of two operators. Since \(V\) is dependent on \(\vec{r}\) and the kinetic part consists of the Laplacian operator, it follows (because \([V(\vec{r}),\partial_{\vec{r}}^2] \neq 0\)) that for the operation \(\exp[i\int(T+V)\Delta t]\) applied to a wave function, it can generally not be expressed in an appropriate eigenbasis, as both parts are not simultaneously diagonalizable.

Therefore, the Hamiltonian is split into two parts to diagonalize them separately. The potential \(V(\vec{r},t)\) is diagonal in \(\vec{r}\)-space, while the kinetic part is diagonal in \(\vec{k}\)-space as \(T(\vec{k})= -\frac{\hbar^2 \vec{k}^2}{2m}\). However, whereas for ordinary scalar quantities \(e^{a+b}=e^ae^b\) holds, it does not for non-commuting operators. For this, we have to use the Lie-Trotter-Suzuki formula. This formula can be derived from the Baker–Campbell–Hausdorff relation which states

\[e^{-\frac{i}{\hbar}(T+V)\Delta t} = e^{-\frac{i}{\hbar}T\Delta t} \cdot e^{-\frac{i}{\hbar}V\Delta t} + e^{-\hbar^{-2}[T,V]\Delta t^2} + \mathcal{O}(\Delta t^3).\]

This equation would lead to accuracy of order \(\Delta t\) if terms with commutators would be neglected. However, the equation can be brought into the form of the Symmetric-Strang-Splitting.

\[e^{-\frac{i}{\hbar}(T+V)\Delta t} = e^{-\frac{i}{\hbar}T\frac{\Delta t}{2}} \cdot e^{-\frac{i}{\hbar}V\Delta t} \cdot e^{-\frac{i}{\hbar}T\frac{\Delta t}{2}} + \mathcal{O}(\Delta t^3)\]

which improves the accuracy by one order, without significant additional computational effort.\ Having split the Hamilton-operator into three operators, which can all be applied subsequently if the wave function is in an appropriate basis, we can evolve the wave function in time with the following algorithm:

The initial wave function is Fourier-transformed \(\psi(\vec{r},t)\rightarrow \mathcal{F} [ \psi(\vec{r},t) ] = \psi(\vec{k},t)\)
The first operator is applied to the wave function \(\psi(\vec{k},t+\frac{\Delta t}{2}) = e^{i\frac{\hbar}{2m}|\vec{k}|^2\frac{\Delta t}{2}} \psi(\vec{k},t)\)
Fourier transformation of \(\psi\) back into \(\vec{r}\)-space \(\psi(\vec{k},t+\frac{\Delta t}{2}) \rightarrow \mathcal{F}^{-1} [ \psi(\vec{k},t+\frac{\Delta t}{2}) ]= \psi(\vec{r},t+\frac{\Delta t}{2})\)
Phase correction due to potential and nonlinearity over the entire interval \(\psi(\vec{r},t+\frac{\Delta t}{2}) \rightarrow e^{-\frac{i}{\hbar} V(\vec{r},t+\frac{\Delta t}{2}) \Delta t} \psi(\vec{r},t+\frac{\Delta t}{2})\)
Fourier-transformation into \(\vec{k}\) -space \(\psi(\vec{r},t+\frac{\Delta t}{2}) \rightarrow \mathcal{F} [\psi(\vec{r},t+\frac{\Delta t}{2}) ] = \psi(\vec{k},t+\frac{\Delta t}{2})\)
The operator of the second kinetic half-time-step is applied \(\psi(\vec{k},t+\Delta t) = e^{i\frac{\hbar}{2m}|\vec{k}|^2\frac{\Delta t}{2}} \psi(\vec{k},t+\frac{\Delta t}{2})\)
Last Fourier transformation into \(\vec{r}\) -space \(\psi(\vec{k},t+\Delta t) \rightarrow \mathcal{F}^{-1} [ \psi(\vec{k},t+\Delta t) ] = \psi(\vec{r},t+\Delta t)\)

Using this algorithm, the initial wave function is propagated in time \(\psi(\vec{r},t)\rightarrow \psi(\vec{r},t+\Delta t)\), just as described by formula of the symmetric Strang splitting.

Importance of sampling frequency¶

The split-step Fourier method uses position and momentum space to calculate the wave-function’s time-propagation. Therefore, one hast to keep in mind that the momentum state of the wave-function can only be correctly accounted for if the space-grid has enough sampling points.

This means that if \(\Delta x\) is the whole position space consisting of \(N\) samples, we have a sampling rate of \(\delta x = \Delta x / N\). If we Fourier transform any wave-form, we can safely represent momentum states within \(\Delta k = 2\pi / \delta x\). The sampling rate in momentum space is \(\delta k = 2\pi / \Delta x\).

Importance of number of sampling points¶

Fast Fourier Transforms (FFTs) are fastet if the number of grid points of your spatial grid is a multiple of \(2\) e.g.: 256, 2048 etc. As every timestep involves at least two FFT we advice you to do so.