Derivation of optical Schrodinger equation

In order to study scattering problems of light in scalar approximation it is possible to use quantum mechanics tools through the analogy between the Helmoltz equation and the Schrodinger equation. In fact, if the permittivity profile is $\epsilon=\epsilon_b+\alpha(z)$ if we start with Helmoltz

$$\nabla^2 E + k_0^2 \epsilon_b E + k_0^2 \alpha(z) E= 0$$

Setting $$E(x,z) = E(x) e^{i k_z z}$$

$$- \frac{d^2 E}{dx^2} - k_z^2 E -k_0^2 \epsilon_b E - k_0^2 \alpha(z) E= 0$$

$$- \frac{d^2 E}{dx^2} - k_0^2 \alpha(z) E = (k_0^2 \epsilon_b - k_z^2) E$$

And defining: $$\hat{T} := -\frac{d}{dx} \qquad \hat{V} :=-k_0^2 \alpha(z) \qquad \mathcal{E}:=k_0^2 \epsilon_b -k^2_z$$

led to an equation similar to a time independent Scrodinger equation: $$(\hat{T}+\hat{V})E=\mathcal{E}E$$

where the time coordinate role is played by $z$.

I used this analogy to perform some simulation on permittivity profiles, but I realized that while this formalism works well for transmission and reflections coefficients, it misses some information about the propagation of the beam, in fact the packets travels twice the speed they should. I want to derive this analogy in a more formal and correct way in order to find an equation similar to the time dependent Schrodinger equation.

The reasoning I tried follows:

Setting: \begin{equation} \tilde{E}_z (x,y,z) = U(x,y,z) e^{i k_z z} \end{equation}

The Helmholtz equation become:

\begin{equation} \frac{d^2U}{dx^2}e^{i k_z z} + \frac{d^2U}{dy^2} e^{i k_z z}+ \frac{d^2U}{dz^2} e^{i k_z z} + 2 i k_z \frac{dU}{dz} e^{i k_z z} - U k_z^2 e^{i k_z z} + k_0^2 \epsilon_r U e^{i k_z z}= 0 \end{equation}

In the SVEA approximation:

\begin{equation} 2 i k_z \frac{dU}{dz} = -\frac{d^2U}{dx^2} - \frac{d^2U}{dy^2} + k_z^2 U - k_0^2 \epsilon_r U \end{equation}

\begin{equation} 2 i k_z \frac{dU}{dz} = -\nabla^2_\perp U + k_z^2 U - k_0^2 \epsilon_b U - k_0^2 \alpha(z) U \end{equation}

\begin{equation} i \frac{1}{k_z} \frac{dU}{dz} = -\frac{1}{2 k_z^2} \nabla^2_\perp U - \frac{1}{2} \left[\frac{k_0^2}{k_z^2} \left(\epsilon_b + \alpha\right)+1 \right] U \end{equation}

This equation seems reasonable because it recover the missing factor $\frac{1}{2}$ in the kinetic term but on the other hand the potential operator seems messy because it contains the previously defined $\mathcal{E}$ term.

Is this the correct way to proceed? Is there a better way? I checked several book on photonics and fiber optics physics but they all use the paraxial approximation that I cannot use for the problem I'm studying.