Operator splitting to solve time dependent Schrödinger equation

Question

I encountered the split operator method to solve the time dependent Schrödinger equation during a lecture. I understand the method on a theoretical basis (I think at least), but I'm struggling to understand one specific aspect of the implementation.

Take a look at the following Pseudo-code (taken from arxiv:1306.3247, page 7): $$ \begin{array}{l} \tilde{\psi}(\mathbf{p}) \leftarrow \mathcal{F} \psi(\mathbf{x}) \\ \text {Multiply } \tilde{\psi}(\mathbf{p}) \text { by } \exp \left[-\frac{i \Delta t}{2 \hbar} \frac{\mathbf{p}^{2}}{2 m}\right] \\ \psi(\mathbf{x}) \leftarrow \mathcal{F}^{-1} \tilde{\psi}(\mathbf{p})\\ \begin{array}{l} \text {for } i \leftarrow 1 \text { to } n-1 \text { do } \\ \quad \begin{array}{l} \text {Multiply } \psi(\mathbf{x}) \text { by } \exp \left[-\frac{i \Delta t}{\hbar} V(\mathbf{x}, t)\right] \\ \tilde{\psi}(\mathbf{p}) \leftarrow \mathcal{F} \psi(\mathbf{x}) \\ \text {Multiply } \tilde{\psi}(\mathbf{p}) \text { by } \exp \left[-\frac{i \Delta t}{\hbar} \frac{\mathbf{p}^{2}}{2 m}\right] \\ \psi(\mathbf{x}) \leftarrow \mathcal{F}^{-1} \tilde{\psi}(\mathbf{p})\end{array} \\ \text { end }\end{array}\\ \text {Multiply } \psi(\mathbf{x}) \text { by } \exp \left[-\frac{i \Delta t}{\hbar} V(\mathbf{x}, t)\right] \\ \tilde{\psi}(\mathbf{p}) \leftarrow \mathcal{F} \psi(\mathbf{x}) \\ \text {Multiply } \tilde{\psi}(\mathbf{p}) \text { by } \exp \left[-\frac{i \Delta t}{2 \hbar} \frac{\mathbf{p}^{2}}{2 m}\right] \\ \psi(\mathbf{x}) \leftarrow \mathcal{F}^{-1} \tilde{\psi}(\mathbf{p})\end{array} $$ For the sake of simplicity let us assume that we are looking at a one-dimensional problem. The algorithm makes use of the Fourier transformation to switch between momentum space and the normal coordinate space. When implementing this one naturally has to switch from the continuous Fourier transformation to a discrete one (ideally FFT). Let $\Delta x$ be a resolution of space, assume $\Delta t$ to be finite and small enough and let $\hbar=m=1$.

So if we define $\psi_{x_n}(t) :=\psi(n\Delta x,t)$ we get $\psi_{p_n}(t):=\mathcal{F}\psi_n(t)$. The algorithm now demands that we multiply this by $\exp(-\frac{i \Delta t}{2} \frac{p_n^2}{2})$, but how exactly do I determine $p_n$ (in each step of the computation)?

Answer 1

I can describe the simplified version of this algorithm which might help you to understand what's going on here. Let's you have this Time-dependent Schrodinger equation:

$$i \hbar \partial_{t} \Psi = -\frac{\hbar^{2}}{2m} \nabla^{2} \Psi$$

For a moment forget about potential cause that complicate things here, which might be necessarily helpful for educational purposes. Now, take Fourier transform from both sides and call $\mathcal{F} \{ \Psi(\mathbf{x},t) \} = \tilde{\Psi}(\mathbf{k},t)$:

$$i \hbar \frac{d \tilde{\Psi}}{dt} = \frac{\hbar^{2}}{2m} \mathbf{k} \cdot \mathbf{k} \tilde{\Psi}$$

The Fourier transform of momentum operator $\mathbf{p} = -i\hbar \nabla$ is: $$\tilde{\mathbf{p}} = \hbar \mathbf{k}$$

So:

$$i \hbar \frac{d \tilde{\Psi}}{d t} = \frac{\tilde{\mathbf{p}}^{2}}{2m} \tilde{\Psi}$$

If you discretize the above equation for different values of wave vector ($\mathbf{k}$) or equivalently the momentum in Fourier space, you have a matrix called $\tilde{P}^{2}$ which represent the discretized version of momentum operator in Fourier space. So:

$$i \hbar \frac{d |\tilde{\Psi} \rangle}{dt} = \frac{\tilde{P}^{2}}{2m} | \tilde{\Psi} \rangle$$

Where $|\tilde{\Psi} \rangle$ is the discretized wave function vector in Fourier space. Now you have this algorithm:

1. $t = 0$: Find Fourier transform of initial condition: $|\tilde{\Psi}_{0} \rangle = \mathcal{F} \{ |\Psi_{0} \rangle \}$.
2. $t = \Delta t$: $|\tilde{\Psi}_{\Delta t} \rangle = \exp{(-\frac{i\Delta t}{\hbar} \frac{\tilde{P}^{2}}{2m})} | \tilde{\Psi}_{0} \rangle$.
3. Repeat item 2 by replacing $\Delta t \rightarrow \Delta t + \Delta t$ and $0 \rightarrow \Delta t$.

It was much easier than your initial algorithm from that paper cause we don't have potential and we don't to use Strang Splitting method because of that. But imagine if you have this more general Time-dependent Schrodinger equation:

$$i \hbar \partial_{t} \Psi = (-\frac{\hbar^{2}}{2m} \nabla^{2} + V(\mathbf{x}))\Psi$$

Now if you take Fourier transform, you need to deal with $\mathcal{F} \{ V(\mathbf{x}) \Psi(\mathbf{x},t) \} = \mathcal{F} \{ V \} * \mathcal{F} \{ \Psi \} = \tilde{V} * \tilde{\Psi}$. Note that $*$ operator is a convolution operator and differs from simple multiplication:

$$i \hbar \frac{d \tilde{\Psi}}{dt} = \frac{\tilde{\mathbf{p}}^{2}}{2m} \tilde{\Psi} + \tilde{V} * \tilde{\Psi}$$

Now let's have closer look at term $\tilde{V} * \tilde{\Psi}$:

$$\tilde{V} * \tilde{\Psi} = \int_{\Omega_{\tilde{\mathbf{p}}}} \tilde{V}(\tilde{\mathbf{p}}^{'}) \tilde{\Psi} (\tilde{\mathbf{p}} - \tilde{\mathbf{p}}^{'}) d^{3} \tilde{\mathbf{p}}^{'} = \sum_{\tilde{\mathbf{p}}^{'}} \tilde{V}(\tilde{\mathbf{p}}^{'}) \tilde{\Psi} (\tilde{\mathbf{p}} - \tilde{\mathbf{p}}^{'})$$

So:

$$i \hbar \frac{d |\tilde{\Psi} \rangle}{dt} = \frac{\tilde{P}^{2}}{2m} | \tilde{\Psi} \rangle + \tilde{\mathcal{V}} | \tilde{\Psi} \rangle$$

Where $\tilde{\mathcal{V}}$ is the potential matrix in Fourier space. Now you have this algorithm by using Strang Splitting:

1. Find Fourier transform of initial condition: $|\tilde{\Psi_{0}} \rangle = \mathcal{F} \{ |\Psi_{0} \rangle \}$.
2. $|\tilde{\Psi}_{\frac{\Delta t}{2}}^{*} \rangle = \exp{(-\frac{i\Delta t}{2\hbar} \frac{\tilde{P}^{2}}{2m})} | \tilde{\Psi}_{0} \rangle$
3. $|\tilde{\Psi}_{\Delta t}^{*} \rangle = \exp{(-i \frac{\Delta t}{\hbar} \tilde{\mathcal{V}})} |\tilde{\Psi}_{\frac{\Delta t}{2}}^{*} \rangle$
4. $|\tilde{\Psi}_{\Delta t} \rangle = \exp{(-\frac{i\Delta t}{2\hbar} \frac{\tilde{P}^{2}}{2m})} | \tilde{\Psi}_{\Delta t}^{*} \rangle$

I hope now it's a bit clearer that what's going on here.