How to find the Riemannian Distance on the Sphere $S^n\subset \Bbb{R}^{n+1}$

I am trying to prove that, if $\rho:S^n\times S^n\to \Bbb{R}$ is the distance induced by the Riemannian metric on the sphere $S^n\subset \Bbb{R}^{n+1}$, then $$\rho(p,q)=\arccos \langle p,q\rangle,\,\,\,\forall p,q\in S^n$$ with $\arccos$ defined from $[-1,1]$ to $[0,\pi]$.

If $p=q$, then the formula is trivial. We have then two cases: $p=-q$ and $p\neq- q$. I am having trouble even with the apparently easier case, $p=-q$.

If $p=-q$, then let $v\in S^n$ be any vector orthogonal to $p$. Then $\alpha:[0,\pi]\to S^n$, $\alpha(t)=(\cos t)p+(\sin t)v$ is a well defined differentiable path from $p$ to $q$ such that $\ell_0^\pi(\alpha)=\pi$ ("length of $\alpha$"). This guarantees that $\rho(p,q)\leq \pi (=\arccos \langle p,q\rangle$, in this case).

In order to show that $\rho(p,q)=\pi$, I must consider an arbitrary differentiable by parts path $\beta:[a,b]\to S^n$ from $p$ to $q$ and show that $$\pi\leq \ell_a^b(\beta)=\int_a^b|\beta'(t)|\,dt.$$

I've done some geometric observations and computations, but without success. How can I do this? I wish I could do it without using any facts about geodesics (they appear later in the book I'm studying with).

Let $p=(1,0,\dots,0)$ and $q=(-1,0,\dots,0)$. Then any path from $p$ to $q$ is of the form $\gamma(t)=(t,h(t))$, where $h(t)\in \Bbb{R}^n$. Then the length of the path is $$ L=\int_{-1}^{1}|\gamma'(t)|dt. $$ But $\gamma'(t)=(1,h'(t))$, and so $|\gamma'(t)|=\sqrt{1+|h'(t)|^2}$.

Now $t^2+\langle h(t),h(t)\rangle=\langle \gamma(t),\gamma(t)\rangle=1$, so $$ 2t+2\langle h(t),h'(t)\rangle=0. $$ By Cauchy Schwartz we have $|h'(t)|^2|h(t)|^2\ge |\langle h(t),h'(t)\rangle|^2=|t|^2$, and so $$ 1+|h'(t)|^2\ge \frac{t^2}{|h(t)|^2}+1, $$ but $|h(t)|^2=1-t^2$, hence $$ \sqrt{1+|h'(t)|^2}\ge \sqrt{1+\frac{t^2}{1-t^2}}=\sqrt{\frac{1}{1-t^2}}, $$ and so $$ L=\int_{-1}^{1}\sqrt{1+|h'(t)|^2}dt\ge\int_{-1}^1 \sqrt{\frac{1}{1-t^2}}dt=\pi. $$

A similar reasoning applies for any other point: Assume $p=(1,0,\dots,0)$, then you can assume that the other point $q$ is of the form $(a,s,0,\dots,0)$. The path $\gamma(t)=(t,\sqrt{1-t^2},0,\dots,0)$ has length $$ L=\int_{a}^{1}|\gamma'(t)|dt=\int_{a}^1 \sqrt{\frac{1}{1-t^2}}dt =-arccos(t)|_{a}^1=arccos(a)=arccos(\langle p,q\rangle). $$ By a similar argument as above, any other path has the form $(t,h(t))$ with $\sqrt{1+|h'(t)|^2}\ge \sqrt{\frac{1}{1-t^2}}$, so the result follows.

To get a lower bound for the length of a path joining two points, you may place one point at the north pole, introduce one angular coordinate and mimic the standard proof from Euclidean space. It is convenient to use a slightly singular coordinate system (essentially the same as you are discussing):

Let $x=(x_0,x_1,\ldots,x_n)\equiv (x_0,\vec{x})\in S^n$ be coordinates on the $n$-sphere and introduce the angular coordinate $\theta=\theta(x)\in[0,\pi]$ by setting $x_0=\cos \theta$. This is a smooth coordinate when $\theta\in(0,\pi)$, i.e. a part from the two poles. When $\theta\in(0,\pi)$ we may also define a complementary set of smooth coordinates using the $(n-1)$-sphere by defining $\vec{z}=(z_1,\ldots,z_n)\in S^{n-1}$ to be the vector that verifies $\vec{x}=\vec{z}\sin \theta$. We may then write: $$(x_0,\vec{x})=(\cos \theta, \vec{z} \sin \theta)$$ Now suppose that $\gamma(t)=(x_0(t),\vec{x}(t))$ is a $C^1$-path in $S^n$ avoiding the poles. Then in our (partially) spherical coordinate system $(\theta,\vec{z})$ we may calculate the derivative: $$ \dot\gamma(t) = (-\sin \theta, \vec{z} \cos\theta) \;\dot\theta(t) \; + \; (0,\dot{\vec{z}} \sin\theta)$$ The last two vectors are orthogonal, because $\vec{z}\cdot \vec{z}=1$ implies $\dot{\vec{z}}\cdot \vec{z}=0$. For the length we then get: $$ |\dot\gamma|^2= \dot\theta^2 + |\dot{\vec{z}}|^2 \sin^2\theta \geq \dot\theta^2 $$ whence the simple lower bound: $$ |\dot\gamma | \geq |\dot\theta|$$ with equality iff $\dot{\vec{z}}=0$. Integrating we see that the distance (minimal path length) between $p$ and $q$ satisfies: $d(p,q)\geq |\theta(p)-\theta(q)|$ a bound which by taking limits also holds when $p$ and/or $q$ is one (or both) of the poles. In particular, if $p=(1,0\ldots,0)$ we get: $d(p,q) \geq \theta(q) \in [0,\pi]$ with equality iff $\dot{\vec{z}}=0$, i.e. when $\vec{z}(t) = \vec{z}_0$ is a constant vector along the path. This extremal path runs along a great circle passing through the north pole. So in particular, the distance between the north and the south pole is exactly $\pi$ with extremal path being any great circle passing through the poles. We also have $\cos(\theta(q))= q_0 = p\cdot q$, so $d(p,q)= \arccos (p\cdot q)$ is the distance, first between the north pole and some other point, but because of rotational invariance of the scalar product and distance the formula holds for any two points on the sphere.