For relativistic massive particle, the action is
$$\begin{align}S ~=~& -m_0 \int ds \cr
~=~& -m_0 \int d\lambda ~\sqrt{ g_{\mu\nu} \dot{x}^{\mu}\dot{x}^{\nu}} \cr
~=~& \int d\lambda \ L,\end{align}$$ where $ds$ is the proper time of the particle; $\lambda$ is the parameter of the trajectory; and we used Minkowski signature $(+,-,-,-)$. So what is the action for a massless particle?
Special Relativity – Lagrangian for Relativistic Massless Point Particle
actiongeodesicslagrangian-formalismspecial-relativityvariational-principle
Related Solutions
This is what happens when you overload notation and mix up what is the curve in spacetime and what are the coordinates on spacetime, and mixup several functions by giving them the same names.
My first suggestion is that you use the notation $v^a$ for the velocity term, so that $x^a,v^a$ are independent. More formally, the Lagrangian is a function (on the tangent bundle of a manifold, whose representation in terms of a chart is) given by \begin{align} L(x,v)&=-mc\sqrt{g_{ab}(x)v^av^b} \end{align} So, \begin{align} \frac{\partial L}{\partial v^a}(x,v) &= mc\frac{g_{ab}(x)v^b}{\sqrt{g_{kl}(x)v^kv^l}}\tag{$*$} \end{align} Note that here $L$ is (as mentioned above) simply a function of two variables $(x,v)\in\Bbb{R}^n\times\Bbb{R}^n$, so it makes sense to talk about the $2n$ partial derivatives $\frac{\partial L}{\partial x^1},\dots,\frac{\partial L}{\partial x^n},\frac{\partial L}{\partial v^1},\dots,\frac{\partial L}{\partial v^n}$.
Next, what we require in the Euler-Lagrange equations is a curve $\gamma:\Bbb{R}\to TM$, which we lift to $\gamma':\Bbb{R}\to TM$. Anyway without getting too much into the differential geometry, once you choose charts (i.e coordinates), you can think of them as standard curves $\gamma:\Bbb{R}\to\Bbb{R}^n$ and its velocity $\dot{\gamma}:\Bbb{R}\to\Bbb{R}^n$. The 'least' action principle states that a curve $\gamma$ satisfies the equations of motion if and only if for every parameter value $\lambda\in\Bbb{R}$, one has \begin{align} \frac{d}{ds}\bigg|_{s=\lambda}\left(s\mapsto \frac{\partial L}{\partial v^a}\bigg|_{(\gamma(s), \dot{\gamma}(s))}\right)-\frac{\partial L}{\partial x^a}\bigg|_{\gamma(\lambda), \dot{\gamma}(\lambda)}&= 0 \end{align} In other words, YOU PERFORM THE PARTIAL DERIVATIVES FIRST, AND ONLY AFTERWARDS PLUG IN THE CURVE (and do $\frac{d}{d\lambda}$). Excuse the caps, but this is I think a point which is completely glossed over in any treatment of Lagrangian mechanics and is one of the biggest sources of confusion. Also, I'm clearly aware that the above expression is very cumbersome to write, but atleast it is 100% unambiguous and is the explicit/precise way of writing things (and not appreciating this is contributing to your confusion).
In your case, you're getting the mysterious "$0=0$" (assuming $g_{ab}$ doesn't actually depend on $x$) because you plugged in the curve first and then (erroneously) tried to do partial derivatives (because your notation makes it "seem ok" when it is infact not). More explicitly, what you said is we can choose a parametrization of the curve such that for each $\lambda$, $g_{ab}(\gamma(\lambda))\dot{\gamma}^a(\lambda)\dot{\gamma}^b(\lambda)=\text{const}$. Ok, that's certainly true, but this means the function $\ell:\Bbb{R}\to\Bbb{R}$ defined as \begin{align} \ell(\lambda)&=L(\gamma(\lambda),\dot{\gamma(\lambda)})\\ &=-mc\sqrt{g_{ab}(\gamma(\lambda))\dot{\gamma}^a(\lambda)\dot{\gamma}^b(\lambda)}\\ &=-mc^2 \end{align} is constant. HOWEVER, $\ell$ and $L$ are NOT at all the same function, but you're thinking that they are. The correct statement is that $\ell'=0$, NOT that $\frac{\partial L}{\partial v^a}=0$, and by extension, the mapping $\lambda\mapsto \frac{\partial L}{\partial v^a}|_{\gamma(\lambda),\dot{\gamma}(\lambda)}$ is not the zero mapping. It is this abuse of notation which is tripping you up completely (you wrote $\frac{\partial L}{\partial \dot{x}^a}=-\frac{\partial mc^2}{\partial \dot{x}}=0$, when if anything you should have written $\frac{\partial \ell}{\partial \dot{x}^a}=0$. But of course, $\ell$ is a function of one variable so using partial derivatives is just nonsense. Hence, by writing things clearly, you can immediately spot where you're making errors). Just for completeness, the correct formula is (according to $(*)$ above) \begin{align} \frac{\partial L}{\partial v^a}\bigg|_{(\gamma(\lambda),\dot{\gamma}(\lambda))}&= mc\frac{g_{ab}(\gamma(\lambda))\,\,\dot{\gamma}^b(\lambda)}{\sqrt{g_{kl}(\gamma(\lambda))\dot{\gamma}^k(\lambda)\dot{\gamma}^l(\lambda)}}\\ &=mc\frac{g_{ab}(\gamma(\lambda))\,\,\dot{\gamma}^b(\lambda)}{c}\\ &=m\cdot g_{ab}(\gamma(\lambda))\,\, \dot{\gamma}^b(\lambda) \end{align} Therefore, you can clearly see now that there is no reason for this to be identically zero for all $\lambda$.
This is a very basic notation and definition issue, which I beat to death in the following MSE answer About the "ambiguity" of explicit and implicit function of a variable (yes this is a basic issue, but I really think you could benefit from being more careful with the basics).
Other Links
Below are some answers I've written which I think could help you clarify some really basic issues (particularly since I explain where exactly the sloppiness of the physics notation is occurring, and where the 'pitfalls'. Once you know these 'pitfalls', you can of course go back to writing sloppily/succinctly)
Differentiation in the Geodesic Problem. Here OP's question is about how one can think of $x, \dot{x}$ as independent variables for partial differentiation (so not directly related to your question here), but I explain what the notation actually stands for, and emphasize the importance of keeping track of what is the function and where it is being evaluated.
Doubt about ordinary and partial derivative. Here, OP is confused about a simple application of the chain rule (of course not your confusion here), but the mistake the OP makes is the same as the one you make, namely to confuse two functions $\Lambda$ and $L$, where the former is obtained by a certain composition with the latter.
More generally, an infinitesimal gauge quasi-symmetry $$\delta q^k(t)~=~\int\! dt^{\prime}~ R^k(t,t^{\prime})\epsilon(t^{\prime})\tag{3.5a}$$ of an action functional $S[q]$, $$\begin{align} \text{possible}&\text{ boundary terms}\cr ~=~&\delta S\cr ~=~&\sum_k\int\! dt\int\! dt^{\prime}~\frac{\delta S}{\delta q^k(t)} R^k(t,t^{\prime})\epsilon(t^{\prime})\cr ~+~&\text{possible boundary terms}\end{align}\tag{3.6a}$$ leads to an off-shell Noether identity $$ \sum_k\int\! dt~\frac{\delta S}{\delta q^k(t)} R^k(t,t^{\prime})~=~0, \tag{3.6b}$$ cf. Noether's 2nd theorem and e.g. Ref. 1.
The Noether identity (3.6b) can be interpreted as the EOMs are perpendicular to the gauge quasisymmetry.
The Noether identity (3.6b) is a generalization of OP's last eq. (6) with the following change of notation $$t~\to~\tau, \qquad q^k~\to~x^{\mu}, \qquad R^k(t,t^{\prime})~\to~-\dot{x}^{\mu}\delta(\tau\!-\!\tau^{\prime}).$$
References:
- M. Henneaux & C. Teitelboim, Quantization of Gauge Systems, 1994; eqs. (3.5a) + (3.6a) + (3.6b).
Best Answer
OP's square root action is not differentiable along null/light-like directions, which makes it ill-suited for a massless particle. So we have to come up with something else. One equation of motion for a scalar massless relativistic point particle on a Lorentzian manifold $(M,g)$ is that its tangent should be null/light-like $$ \dot{x}^2~:=~g_{\mu\nu}(x)~ \dot{x}^{\mu}\dot{x}^{\nu}~\approx ~0, \tag{A}$$ where dot denotes differentiation wrt. the world-line (WL) parameter $\tau$ (which is not proper time). [Here the $\approx$ symbol means equality modulo EOM.] This suggests that a possible action is $$ S[x,\lambda ]~=~\int\! d\tau ~L, \qquad L~=~\lambda ~\dot{x}^2, \tag{B} $$ where $\lambda(\tau)$ is a Lagrange multiplier. This answer (B) may seem like just a cheap trick. Note however that it is possible by similar methods to give a general action principle that works for both massless and massive point particles in a unified manner, cf. e.g. Ref. 1 and eq. (3) in my Phys.SE answer here.
More importantly, the corresponding Euler-Lagrange (EL) equations for the action (B) are the null/light-like condition $$ 0~\approx ~\frac{\delta S}{\delta\lambda}~=~\dot{x}^2, \tag{C}$$ and the geodesic equations $$ 0~\approx ~-\frac{1}{2}g^{\sigma\mu}\frac{\delta S}{\delta x^{\mu}} ~=~\frac{d(\lambda\dot{x}^{\sigma})}{d\tau} +\lambda\Gamma^{\sigma}_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu},\tag{D} $$ as they should be.
The action (B) is invariant under WL reparametrization $$ \begin{align} \tau^{\prime}~=~&f(\tau), \qquad d\tau^{\prime} ~=~ d\tau\frac{df}{d\tau},\cr \dot{x}^{\mu}~=~&\dot{x}^{\prime\mu}\frac{df}{d\tau},\qquad \lambda^{\prime}~=~\lambda\frac{df}{d\tau}.\end{align}\tag{E} $$ Therefore we can choose the gauge $\lambda=1$. Then eq. (D) reduces to the more familiar affinely parametrized geodesic equations.
References: