As discussed here and within, the Lagrangian for a massless particle, using the $(-,+,+,+)$ metric signature, is
$$L = \frac{\dot{x}_\mu \dot{x}^\mu}{2e} - V,\tag{1}$$
where $\dot{x}^\mu := \frac{dx^\mu}{d\lambda}$ is the velocity, $\lambda$ is some worldline parameter, $e$ is the auxiliary einbein and $V$ is the potential.
The EL equations give us the EOMs
$$\dot{x}_\mu \dot{x}^\mu = 0,\tag{2}$$ $$\ddot{x}^\mu + \Gamma^\mu_{\sigma\rho} \dot{x}^\sigma \dot{x}^\rho - \frac{\dot{e}\dot{x}^\mu}{e} + e\partial^\mu V = 0,\tag{3}$$
where $\Gamma^\mu_{\sigma\rho}$ are the Christoffel symbols of the metric $\eta_{\mu\nu}$ for some choice of coordinates. After this, I'm not sure how to proceed, for the following reasons.
In the $V=\text{constant}$ case, the system is underdetermined, and we are free to choose some $e$, such as setting $e=1$. We then get the consistent EOMs
$$\dot{x}_\mu \dot{x}^\mu = 0,\tag{4}$$ $$\ddot{x}^\mu + \Gamma^\mu_{\sigma\rho} \dot{x}^\sigma \dot{x}^\rho = 0.\tag{5}$$
In the general case, however, we seem to end up with the two EOMs being inconsistent.
For example, suppose we use Cartesian coordinates and a $z$-direction potential such as $$V = z,\tag{6}$$ choosing initial conditions satisfying the null velocity condition such as $\dot{x}^\mu \left(0\right) = \left(1,0,0,1\right)$. The EOM for the coordinates becomes
$$\ddot{x}^\mu = \left(0,0,0,-1\right),\tag{7}$$
yielding $\dot{x}^\mu \left(\lambda\right) = \left(1,0,0,1-\lambda\right)$, which fails to satisfy the null velocity condition for $\lambda \neq 0$.
