Deriving the Euler-Lagrange equations for the arclength of a curve on the unit sphere I'm trying to derive the extremal solutions to the Lagrangian for arclength on the unit sphere by setting up the Euler-Lagrange equations.
Starting from
$$L = \sqrt{\dot \theta^2 + \sin^2 \theta \ \dot \phi^2} \ $$
I see not long after that
$$\frac{\partial L}{\partial \dot\theta}=\frac{\dot\theta}{\sqrt{\dot \theta^2 + \sin^2 \theta \ \dot \phi^2}}\\
\frac{\partial L}{\partial\theta}=\frac{\dot\phi^2\sin(\theta)\cos(\theta)}{\sqrt{\dot \theta^2 + \sin^2 \theta \ \dot \phi^2}}\\
\frac{\partial L}{\partial\dot\phi}=\frac{\sin^2(\theta)\dot\phi}{\sqrt{\dot \theta^2 + \sin^2 \theta \ \dot \phi^2}}\\
\frac{\partial L}{\partial\phi}=0$$
I'm pretty sure I've understood the planar case, and that had taught me to recognize that in the last two equations, I could then say that $$\frac{\partial L}{\partial\dot\phi}=C_1$$ a constant $C_1$.
Then one can rewrite 
$$\frac{\partial L}{\partial\theta}=C_1\dot\phi\cot(\theta)\\
\frac{\partial L}{\partial\dot\theta}=C_1\frac{\dot\theta}{\sin^2(\theta)\dot\phi}$$
but then $$\frac{d}{dt}\frac{\partial L}{\partial\dot\theta}=\frac{\partial L}{\partial \theta}$$ still seems to be a second order differential equation in both variables, and I feel like I've run aground.
How do you go from here? Maybe I'm supposed to do more at the step where I found $C_1$ to make some reductions?  
The final goal is to see concretely that the answer is going to yield a single solution when I'm going from $(0,0)$ and $(\pi/2, \pi/2)$, say, and that there are infinitely many extremals when going between antipodal points.
 A: When minimizing length, there are always infinitely many extremals because reparameterizations do not change the length. For this reason, one usually studies the better behaved extremal problem $\int |\dot \gamma|^2\to \min$ (rather than $\int |\dot \gamma|\to \min$). The extremals are geometrically the same, since the Cauchy-Schwarz inequality 
$$
\left(\int_0^1 |\dot \gamma| \right)^2 \le \int_0^1 |\dot \gamma|^2
$$
turns into equality for constant-speed parameterizations. And without the square root we get a lot simpler equations. So, $L=\dot\theta^2 + \sin^2\theta \,\dot \phi^2$  leads to 
$$
\frac{d}{dt}\frac{\partial L}{\partial\dot\theta} = \frac{d}{dt}(2\dot\theta) = 2\ddot\theta
$$
$$
\frac{\partial L}{\partial\theta} =  2\sin\theta\cos\theta \, \dot\phi^2
$$
hence 
$$
\ddot\theta = \sin\theta\cos\theta \, \dot\phi^2  \tag1
$$
Also, 
$$
\frac{\partial L}{\partial\dot\phi} = 2\sin^2\theta \, \dot\phi
$$
$$
\frac{\partial L}{\partial\phi} =  0
$$
hence 
$$
\sin^2\theta \, \dot\phi = C
\tag2 $$
From (1) and (2) we get 
$$
\ddot\theta = C^2 \frac{\cos\theta}{\sin^3\theta}
\tag 3$$
which is still unpleasant but is somewhat solvable for the inverse function $t(\theta)$. When working with  (3) by hand, it seems advisable to let $y=\cot \theta$ and rewrite the equation in terms of $y$.
I somehow doubt that this approach will shed any light on the (non)-uniqueness of geodesics.
