= \left| \boldsymbol{v}\right|$, $r = \left| \boldsymbol{r}\right|$, $\mu$ is the gravitational parameter of the system and $\varepsilon > 0$ is the constant acceleration imparted to the test particle by the force field. Without loss of generality, the coordinate system has been oriented so that the force field is directed towards the positive $z$ axis.

Following the lead of @epstein_zur_1916 and @born_mechanics_1927, we proceed by expressing the Lagrangian in parabolic coordinates $\left( \xi, \eta, \phi \right)$ via the coordinate transformation $$\begin{aligned}
x & =\xi\eta\cos\phi, & \dot{x} & =\left(\dot{\xi}\eta+\xi\dot{\eta}\right)\cos\phi-\xi\eta\dot{\phi}\sin\phi,\label{eq:x_para}\\
y & =\xi\eta\sin\phi, & \dot{y} & =\left(\dot{\xi}\eta+\xi\dot{\eta}\right)\sin\phi+\xi\eta\dot{\phi}\cos\phi,\\
z & =\frac{\xi^{2}-\eta^{2}}{2}, & \dot{z} & =\xi\dot{\xi}-\eta\dot{\eta},\label{eq:z_para}\end{aligned}$$ where $\xi \geq 0$, $\eta \geq 0$ and $\phi \in \left( -\pi,\pi\right]$ is the azimuthal angle. The inverse transformation from cartesian to parabolic coordinates is $$\begin{aligned}
\xi & =\sqrt{r+z}, & \dot{\xi} & =\frac{\dot{r}+\dot{z}}{2\sqrt{r+z}},\label{eq:xi_inv}\\
\eta & =\sqrt{r-z}, & \dot{\eta} & =\frac{\dot{r}-\dot{z}}{2\sqrt{r-z}},\label{eq:eta_inv}\\
\phi & =\arctan\left(y,x\right), & \dot{\phi} & =\frac{\dot{y}x-\dot{x}y}{x^{2}+y^{2}}\label{eq:phi_inv},\end{aligned}$$ where $\dot{r}=\left(\boldsymbol{v}\cdot\boldsymbol{r}\right)/r$ and $\arctan$ is the two-argument inverse tangent function. In the new coordinate system, $$\begin{aligned}
v^{2} & =\left(\xi^{2}+\eta^{2}\right)\left(\dot{\xi}^{2}+\dot{\eta}^{2}\right)+\xi^{2}\eta^{2}\dot{\phi}^{2},\\
r & =\frac{\xi^{2}+\eta^{2}}{2},\end{aligned}$$ and the Lagrangian becomes $$\begin{gathered}
L=\frac{1}{2}\left[\left(\xi^{2}+\eta^{2}\right)\left(\dot{\xi}^{2}+\dot{\eta}^{2}\right)+\xi^{2}\eta^{2}\dot{\phi}^{2}\right]\\
+\frac{2\mu}{\xi^{2}+\eta^{2}}+\varepsilon\frac{\xi^{2}-\eta^{2}}{2}.\end{gathered}$$ Switching now to the Hamiltonian formulation through a Legendre transformation, the momenta are defined as $$\begin{aligned}
p_{\xi} & =\frac{\partial L}{\partial\dot{\xi}}=\left(\xi^{2}+\eta^{2}\right)\dot{\xi},\label{eq:p_xi}\\
p_{\eta} & =\frac{\partial 