Let $(M,g)$ a manifold and $\gamma (t)$ for $t\in [a,b]$ a curve $\mathcal C^1$ the is minimizing the length. Then, if $p=\gamma (t_0)$ and $q=\gamma (t_1)$, then $\gamma $ is also minimizing the length of $p$ and $q$ for all $a\leq t_0<t_1\leq b$. In polar coordinates in $B(p,r)$, we have $$\gamma |_{[t_0,t_1]}=(r(t),\gamma ^1(t),...,\gamma ^{n-1}(t))$$ and
\begin{align*} d(p,q)&=\ell(\Gamma)|_{[t_0,t_1]}\\ &=\int_{t_0}^{t_1}\|\dot \gamma (t)\|\mathrm d t\\ &=\int_{t_0}^{t_1}\sqrt{\dot r(t)+g_{ij}(\gamma (t))\dot y^i\dot y^j}\mathrm d t\\ &\geq \int_{t_0}^{t_1}\sqrt{\dot r(t)^2}\\ &=r(t_1)-r(t_0)\\ &=d(p,q). \end{align*}
We have equality if and only if $\dot y^i=0$ for all $i$, i.e. if $y^i=const$, which proves the claim.
Question
I don't understand the proof, why we have show that $\gamma $ is a geodesic ? To me, we only proved that $d(p,q)\geq d(p,q)$ (which is in fact obvious). Any way, any explanation is welcome.
In polar coordinates $(r,\theta_1, \dots, \theta_n)$, the rays $t\mapsto (t,a_1,a_2, \dots, a_n)$ (where the $a_i$ are constant) are known to be geodesics, and the same applies to reparametrizations $t\mapsto (\phi(t),a_1,a_2, \dots, a_n)$ .
The sketch of proof you replicated shows that, in polar coordinates, the length minimizer is of this form.