I would like to solve the Helmholtz equation of the following form (with $U=U(r, z)$ in radial coordinates): $$ \partial_zU-i\varepsilon\partial_z^2U=i\nabla^2_rU$$ When setting $\varepsilon=0$, I can solve the equation step-wise using a Hankel transformation: $$U(r, z+\Delta z) = H^{-1}\left(\exp\left(-2i\pi^2\rho^2\Delta z\right)H(U(r, z))\right)$$ using the knowledge that $$H(i\nabla^2_rU)=-2i\pi^2\rho^2H(U)$$ Now, according to "Banerjee, P. P., Nehmetallah, G., & Chatterjee, M. R. (2005). Numerical modeling of cylindrically symmetric nonlinear self-focusing using an adaptive fast Hankel split-step method. Optics Communications, 249(1–3), 293–300. https://doi.org/10.1016/j.optcom.2004.12.048", the equation can be modified for $\varepsilon\neq0$, such that $$ U(r, z+\Delta z)=H^{-1}\left(\exp\left(-\frac{i}{2\varepsilon}\left(1-\sqrt{1-16\varepsilon\pi^2\rho^2}\right)\Delta z\right)H(U(r, z))\right)$$ How do I get to that step from the equation above for including the second derivative?
Update:
Following the explanation in the answer, I get as initial ODE
$$-i\varepsilon\partial_z^2U+\partial_zU-i\nabla^2_rU=0$$
which then can be solved using
$$U_{1,2}=\frac{-b\pm\sqrt{b^2-4ac}}{2a}$$
or when using the prefactors from the equation above:
$$ \begin{split}U_{1, 2}&=\frac{-1\pm\sqrt{1-4\cdot\left(-i\varepsilon\cdot2i\pi^2\varrho^2\right)}}{-2i\varepsilon}\\
&=-\frac{i}{2\varepsilon}\left(1\mp\sqrt{1-8\varepsilon\pi^2\varrho^2}\right)
\end{split}$$
After I am only interested in the first solution, I get
$$U(\varrho, z)=A\exp\left(\frac{iz}{2\varepsilon}\left(-1+\sqrt{1-8\varepsilon\pi^2\varrho^2}\right)\right)$$
How exactly did you get $16$ as prefactor, instead of $8$?
Thanks!
Taking the Hankel transform of the equation, with $\tilde{U}(\rho,z)=H\left( U(r,z) \right)$, you get \begin{equation} i\varepsilon\partial_z^2\tilde{U}(\rho,z)-\partial_z\tilde{U}(\rho,z)-2i\pi^2\rho^2\tilde{U}(\rho,z)=0 \end{equation} This is a 2nd order ODE on $\tilde{U}(\rho,z)$ considered as a function of $z$. It can be easily solved: \begin{equation} \tilde{U}(\rho,z)=A\exp\left(\frac{z}{2\varepsilon}\left( -i+\sqrt{16\pi^2\rho^2\varepsilon-1} \right) \right)+B\exp\left(-\frac{z}{2\varepsilon}\left( -i+\sqrt{16\pi^2\rho^2\varepsilon-1} \right) \right) \end{equation} where $A$ and $B$ are some constants. If you have some physical reason to suppose that $B=0$ (for example, considering the behaviour of the solution for $z\to\infty$, taking care of the determination of the square root), then \begin{align} \tilde{U}(\rho,z)&=A\exp\left(\frac{z}{2\varepsilon}\left( -i+\sqrt{16\pi^2\rho^2\varepsilon-1} \right) \right)\\ \tilde{U}(\rho,z+\Delta z)&=A\exp\left(\frac{z+\Delta z}{2\varepsilon}\left( -i+\sqrt{16\pi^2\rho^2\varepsilon-1} \right) \right) \end{align} Then forming the ratio of the preceding expressions and taking the Hankel transform you obtain the quoted result.