Null geodesic correct Lagrangian I am aware that the Lagrangian for a relativistic massless particle is different than that for a massive particle, as the usual action (dots denote derivative w.r.t. $\lambda$)
$$S =m\int d\lambda \sqrt{\dot{x}^\mu\dot{x}^\nu g_{\mu\nu}}\tag{1}$$ is not differentiable for massless particles (also, the prefactor $m$ is 0). What we usually do is introduce an einbein (as outlined in this answer) and then we can take the limit $m\rightarrow 0$.
That gives us the massless Action of the form
$$S = \int d\lambda~ e(\lambda)\dot{x}^\mu\dot{x}^\nu g_{\mu\nu}\tag{2}$$ where $e(\lambda)$ is the einbein which can be chosen at will by changing the parametrisation.
My problem is the following: how do we determine the right einbein given we chose a parametrisation? As an example, say I have the Schwarzschild metric
$$ds^2 = -(1-\frac{2}{r})dt^2+\frac{dr^2}{1-\frac{2}{r}}+r^2d\theta^2+r^2\sin(\theta^2)d\phi^2\tag{3}$$
and I decide to parametrise my lightlike curves as $(t,r(t),\theta(t),\phi(t))$, how can I determine the right einbein $e(t)$?
 A: The (inverse) einbein $e$ is undetermined for at least 2 reasons:

*

*OP's action (2) is invariant under a worldline (WL) reparametrization gauge symmetry $$\begin{align}\lambda^{\prime}~=~&f(\lambda), \qquad 
d\lambda^{\prime} ~=~ d\lambda\frac{df}{d\lambda},\cr 
\dot{x}^{\mu}~=~&\dot{x}^{\prime\mu}\frac{df}{d\lambda},\qquad 
e^{\prime}~=~e\frac{df}{d\lambda}. \end{align}\tag{A}$$


*If we multiply the action (2) with a non-zero constant (which we can absorb into $e$), it does not change the EOMs.
To get a unique $e$, one would therefore have to make a gauge choice. OP considers the static gauge $x^0=\lambda$. In that case $e$ is determined by the EOMs up to an overall constant, cf. the second reason.
A: The einbein is just a Lagrange multiplier that enforces the mass-shell condition. This is most clearly seen in Hamiltonian formulation:
\begin{align}
S[x,p,e] &= \int d\sigma\,\, p_\mu \dot{x}^\mu - e\,\phi(x,p)
\\
0=\frac{\delta S}{\delta e} &\implies \phi=0,
\\
0=\frac{\delta S}{\delta x},\,\, 0=\frac{\delta S}{\delta p} &\implies \dot{f}=e\{f,\phi\},\,\,\{x^\mu,x^\nu\}=0,\,\,\{x^\mu,p_\nu\}=\delta^\mu_\nu,\,\,\{p_\mu,p_\nu\}=0.
\end{align}
$\phi(x,p)=\frac{1}{2}(g^{\mu\nu}(x) p_\mu p_\nu +m^2)$ is a Hamiltonian constraint that functions as both the time-evolution generator and a constraint on the phase space. It is also the generator of reparametrization gauge transformation, so we also call it by "gauge generator." When $m=0$, we get a massless relativistic particle.
Besides, phase space is just another name (physicists' name) for symplectic manifold. Constrained mechanics on phase space is a matter of describing a symplectic submanifold in a symplectic manifold. Hence we not only consider the constraint $\phi(x,p)$ but its pair constraint $\chi(x,p)$. For instance it is customary to let the gauge-fixing condition depend on the worldline time parameter and set $\chi(x,p,\sigma) = -u_\mu x^\mu - \sigma$ for a constant timelike one-form $u_\mu$. That is, we parameterize the worldline with the coordinate time in the co-moving frame of $u^\mu$. (In the "covariant phase space" setting, we set $\chi(x,p)=-u_\mu x^\mu$ and let $\chi(x,p)=0$ describe a time slice that the particle's initial position $x^\mu(0)$ takes its values.) By introducing this "gauge-fixing condition," the total number of constraints becomes an even number; recall that symplectic manifolds are even-dimensional.
In general, Lagrange multipliers can be determined by studying the following three conditions.
\begin{equation}
(1)\,\,\, \phi=0
\qquad
(2)\,\,\, \dot{\phi}=0
\qquad
(3)\,\,\, \dot{\chi}=0
\end{equation}
And in the case of $(2)\, \dot{\phi}=0$, it holds identically because $\dot{\phi} = e \{\phi,\phi\}=0$. So it does not give new information. Usually, $e(\sigma)$ is determined by enforcing (3). This gives the celebrated Dirac brackets. For example you may find it helpful to have a look at appendix A of arXiv:2102.07063.
The method (3) works for both massive and massless particles, so my answer to your question will be: use not only the mass-shell constraint but also a "gauge-fixing constraint" (chosen with your taste)!
The thing is, for massive particles, the einbein can also be determined by using (1) (which gives $e(\sigma) = \frac{1}{m} \sqrt{-\dot{x}^2(\sigma)}$), but for massless particles enforcing the mass-shell constraint gives nothing. So I think there is some kind of subtle difference between massive and massless particles, and this might be the confusing point addressed by your question.
