I am trying to numerically solve the 3D Schrodinger equation with an extra non-linear term (the Gross-Pitaevskii equation) using the split-step operator method. This is a follow up to a previous question which was answered nicely, but now I would like to do the same in cylindrical coordinates (to exploit a symmetry that is already in my problem).
The equation is of the form $$ \frac{\partial\Psi(\rho,z,t)}{\partial t} = -\frac{i}{\hbar}\hat{H}\;\Psi(\rho,z,t) \tag{1} $$ where $$ \hat{H} = \hat{T} + \hat{V} = -\frac{\hbar^2}{2m}\nabla^2 + V(\rho,z), \tag{2} $$ and $\rho$ and $z$ are the cylindrical coordinates (I am assuming no angular dependence in the equation), and $\nabla^2$ is the 3D laplacian.
If I split the operator $\hat{H}$ symetrically (Strang), then when solving for the $\hat{T}$ part (kinetic) over a full step $\Delta t$, I obtain the following $$ \begin{align} \frac{\partial \psi}{\partial t} &= -\frac{i}{\hbar} \hat{T} \psi \tag{3}\\ &= \frac{i\hbar}{2m}\nabla^2\psi . \tag{4} \end{align} $$ In cartesian coordinates, the equation is easy to solve by transforming the whole equation into Fourier space, as described here. However, if I try to follow that approach here using the expression for the Laplacian in cylindrical coordinates I obtain $$ \begin{align} \frac{\partial \psi}{\partial t} &= \frac{i\hbar}{2m}\nabla^2\psi \tag{5}\\ &= \frac{i\hbar}{2m} \Bigg( \frac{\partial^2\psi}{\partial \rho^2} + \frac{1}{\rho}\frac{\partial\psi}{\partial \rho} + \frac{\partial^2\psi}{\partial z^2} \Bigg) \tag{6}\\ \implies \frac{\partial \hat{\psi}}{\partial t} &= \mathcal{F}\Bigg\{ \frac{i\hbar}{2m} \Bigg( \frac{\partial^2\psi}{\partial \rho^2} + \frac{1}{\rho}\frac{\partial\psi}{\partial \rho} + \frac{\partial^2\psi}{\partial z^2} \Bigg) \Bigg\} \tag{7}\\ &= -\frac{i\hbar}{2m} \Bigg( k_\rho^2 \hat{\psi} + k_z^2\hat{\psi} + \mathcal{F}\Bigg\{ \frac{1}{\rho}\frac{\partial\psi}{\partial \rho} \Bigg\} \Bigg) \tag{8} \end{align} $$ where $\hat{\psi} = \hat{\psi}(k_\rho,k_z,t)$ is the Fourier transform of $\psi(\rho,z,t)$, and $k_\rho$ and $k_z$ are the corresponding k-vectors. The last term here now poses a problem - the equation does not have the simple exponential solution as in the cartesian case.
Is there a known way to handle this situation? Or is the split-step method somehow only suited to cartesian coordinate systems? Any advice would be great!
(As mentioned, the aim is to numerically solve the equation eventually, where in the cartesian case one would usually use the FFT algorithm.)
This is a partial answer that should get you going in the right direction. I have some uncertainty about satisfying the (unstated) boundary conditions and the sign conventions in this field, so there is significant room for clarifying this answer.
The result from taking the Fourier transform of both sides can also be arrived at by substituting the spectral representation of the wave function: $$ \psi(\mathbf{x},t) = \int\hat\psi(\mathbf{k},t)\exp(-i \mathbf{k}\cdot\mathbf{x}) \;d\mathbf{k} $$ I'm not 100% sure about the Fourier/time/phase convention in this field, but the sign convention above is one way to do things. In cartesian, we have $$ \nabla^2 \psi(\mathbf{x},t) = \int\hat\psi(\mathbf{k},t)\nabla^2\Big(\exp(-i \mathbf{k}\cdot\mathbf{x})\Big) \;d\mathbf{k}, $$ and $$ \nabla^2\Big(\exp(-i \mathbf{k}\cdot\mathbf{x})\Big) = -k^2\exp(-i \mathbf{k}\cdot\mathbf{x}), $$ and so the Fourier transform turned all the spatial derivatives into constants, simplifying things enough to directly write a solution for $\hat\psi(\mathbf{k},t)$ as the solution to $$ \frac{\partial \hat\psi}{\partial t}=-\frac{i\hbar}{2 m}k^2\hat\psi $$ Now in cylindrical, I think we need a spectral representation that makes the cylindrical Laplacian into a constant. In other words, we're looking for eigenfunctions of the cylindrical Laplacian (with this extra symmetry that nothing depends on the azimuthal angle). $$ \nabla^2\psi = -k^2\psi, $$ with $$ \nabla^2 = \frac{\partial^2}{\partial\rho^2} + \frac{1}{\rho}\frac{\partial}{\partial\rho} + \frac{\partial^2}{\partial z^2} $$ That's the Helmholtz equation and simple separation of variables gives a general solution to the eigenfunction problem: $$ \psi = \Big(A\exp(-ik_z z)+B\exp(ik_z z)\Big)\Big(C\mathrm{H}_0^{(1)}(k_\rho \rho)+D\mathrm{H}_0^{(2)}(k_\rho \rho)\Big). $$ The H functions are the usual Hankel functions. I've also named my separation constants such that $k_\rho^2 + k_z^2 = k^2=|\mathbf{k}|^2$. Okay, so some superposition of waves of the above form with various values for the four constants will be the general solution for $\psi$ which will be our replacement for the Fourier transform. It's here that we usually start appealing to the physical nature of the problem and boundary conditions. There are probably sources of the wavefunction in this problem, and those sources make waves that propagate away from the sources and off to infinity. In classical EM (the field I'm more familiar with) we appeal to some kind of radiation condition, like a boundary condition at infinity, to say that waves can't enter from infinity, only go out. That effectively excludes one of the exponential functions from the solution. Which one is excluded depends on if $z<0$ or $z>0$. For the Hankel functions, I assume you want finiteness of $\psi$ at $\rho=0$ in general since there is nothing special there, so you have to have $C=D=1/2$ to cancel the Bessel functions of the second kind that have a singularity at $\rho=0$, so taken together we can write our superposition integral of the "physical" eigenfunctions (considering only $z>0$, the upper half-plane) as $$ \psi(\mathbf{x},t) = \int \hat\psi(\mathbf{k},t)\exp(-i k_z z)J_0(k_\rho\rho)\;d\mathbf{k} $$ This is our new spectral transform that replaces the Fourier transform. Note I could be wrong about the sign in the exponential, because ultimately that's a matter of convention. I also could be wrong about the presence of the $J_0$ Bessel function; I think the exact form you want here depends exactly on your boundary conditions. In any case, you can modify accordingly with $J_0$, $Y_0$, $H_0^{(1)}$, and $H_0^{(2)}$. The same goes for the exponential, I guess it could be replaced with $\exp(i k_z z)$, $\cos(k_z z)$, or $\sin(k_z z)$ depending on exactly what you're doing and your radiation convention / boundary condition at infinity. For a source of waves sitting at the origin, you can even replace $z$ with $|z|$ for certain kinds of Green's function problems with an impulse at the origin. But lets plug my form into the equation of interest: $$ \frac{\partial \psi}{\partial t}=\frac{i\hbar}{2 m}\nabla^2\psi $$ Term by term, we have: $$ \frac{\partial \psi}{\partial t} = \int \frac{\partial\hat\psi(\mathbf{k},t)}{\partial t}\exp(-i k_z z)J_0(k_\rho\rho)\;d\mathbf{k} $$ $$ \nabla^2\psi = \int (-k^2)\hat\psi(\mathbf{k},t)\exp(-i k_z z)J_0(k_\rho\rho)\;d\mathbf{k} $$ This last one is by construction and was the point of finding the eigenfuctions, so that this Laplacian turned into multiplication by a constant. Any of the alternate forms with Hankel functions or sines and cosines also of course have this property. Putting those two together: $$ \int \frac{\partial\hat\psi(\mathbf{k},t)}{\partial t}\exp(-i k_z z)J_0(k_\rho\rho)\;d\mathbf{k} = \int \left(-k^2\frac{i\hbar}{2 m}\hat\psi(\mathbf{k},t)\right)\hat\psi(\mathbf{k},t)\exp(-i k_z z)J_0(k_\rho\rho)\;d\mathbf{k} $$ and collecting on one side $$ \int \left(\frac{\partial\hat\psi(\mathbf{k},t)}{\partial t}+k^2\frac{i\hbar}{2 m}\hat\psi(\mathbf{k},t)\right)\exp(-i k_z z)J_0(k_\rho\rho)\;d\mathbf{k} = 0 $$ Its sufficient for the term in parentheses to be identically zero $$ \frac{\partial\hat\psi(\mathbf{k},t)}{\partial t}+k^2\frac{i\hbar}{2 m}\hat\psi(\mathbf{k},t) = 0 $$ So we have transformed the equation into a spectral domain where the solutions are indeed simple exponentials in time; we just had to use Hankel AND Fourier transforms. Now when you go to implement this numerically, that's a 1D Fourier transform (can use the FFT) followed by a 1D Hankel transform (you will have to artificially divide the thing you're transforming by a factor of $k_\rho$ to get it in the right form). It turns out there are indeed "fast" Hankel transforms (FHT) buts good codes are much more niche than the FFT.
You can continue to generalize to all four terms instead of retaining the single term as I did. Then you'd have a transform that maps a single function of space and time to a vector function on four dimensions of wavenumber and time. I'm not really sure here but wanted to give you some ideas.