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