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)$$?

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.

The (inverse) einbein $$e$$ is undetermined for at least 2 reasons:

1. 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}

2. 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.