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.

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