Operator splitting to solve time dependent Schrödinger equation

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

• Note that in continuous version exponential of momentum operator or any other operator will be translated to exponential of matrices in discrete version, so multiplication means matrix multiplication. Mar 1 '20 at 14:56
• @AloneProgrammer Do I understand correctly that you mean that in the discrete version $\vec{p} ^2 = diag(p_1^2,\dots, p_n^2)$? I think I understand that, but I still don't know how to find the $p_n$.
– Sito
Mar 1 '20 at 15:00
• The continous momentum operator in real space is: $$\mathbf{p}= -\frac{\hbar^{2}}{2m} \nabla$$and in Fourier space: $$\tilde{\mathbf{p}} = -i \frac{\hbar^{2}}{2m} \mathbf{k}$$ when you discretize these operators in real or Fourier spaces, you’ll get a vector and as a result your Hamiltonian would be a matrix. Mar 1 '20 at 15:17
• I'm sorry, but I do understand that the term $\exp(-\frac{i \Delta t}{2} \frac{\mathbf{p}^2}{2})$ and $\exp (-\frac{i \Delta t}{\hbar} V(\mathbf{x}, t))$ are matrices. I also understand how one calculates $\exp (-\frac{i \Delta t}{\hbar} V(\mathbf{x}, t))$ numerically (matrix exponential), but I don't see how one calculates $\exp(-\frac{i \Delta t}{2} \frac{\mathbf{p}^2}{2})$. Do I just calculate $\exp(-\frac{i \Delta t}{2} \frac{\nabla^2}{2})$ using finite differences and then somehow "Fourier transform" it (even tough I have no idea how that would even look like), (cont.)
– Sito
Mar 1 '20 at 15:54
• @AloneProgrammer (cont.) or do we calculate the $p_i$ separately and then just calculate $\exp(-\frac{i \Delta t}{2} \frac{\mathbf{p}^2}{2})$? The questions boils down to understanding if we need the $p_i$ explicitly (and if, how to calculate them), or if we can somehow directly arrive at the quantity $\exp(-\frac{i \Delta t}{2} \frac{\mathbf{p}^2}{2})$.
– Sito
Mar 1 '20 at 15:55

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.