The Euclidean metric in $\mathbb{R}^3$ induces a metric on a hyperboloid, and the shortest path will be shortest with respect to this distance. Standard coordinates on a hyperboloid $z^2 - x^2 - y^2 = 1$ would be
$$
x = \sinh(t) \sin \phi \qquad
y = \sinh(t) \cos \phi \qquad
z = \cosh(t)
$$
with the interval:
$$
\mathrm{d}s^2 = \cosh(2 t) \mathrm{d} t^2 + \sinh^2(t) \mathrm{d} \phi^2
$$
That means that non-zero Christoffel symbols are
$$
\Gamma^t_{tt}=\tanh(2t) \qquad \Gamma^t_{tt}=-\frac{1}{2}\tanh(2t) \qquad
\Gamma^\phi_{t\phi} = \Gamma^\phi_{\phi t} = \frac{1}{\tanh(t)}
$$
Geodesic equations are readily obtained:
$$
t^{\prime\prime}(s) + \tanh(2 t(s)) \left( (t^\prime(s))^2 - \frac{1}{2} (\phi^\prime(s))^2 \right) = 0 \qquad
\phi^{\prime\prime}(s) + \frac{2}{\tanh(t(s))} t^\prime(s) \phi^\prime(s) = 0
$$
It is not hard to see, that the latter equation admits an integral of motion, i.e. $\phi^\prime(s) \sinh^2(t(s)) = \mathcal{L}$, because
$$
\frac{\mathrm{d}}{\mathrm{d} s} \left( \phi^\prime(s) \sinh^2(t(s)) \right) = \sinh^2(t(s)) \left( \phi^{\prime\prime}(s) + \frac{2}{\tanh(t(s))} t^\prime(s) \phi^\prime(s) \right) \left. = \right|_{\text{eq. of motion}} = 0
$$
Unfortunately this does not get us any closer to the geometric interpretation of the geodesic. Is it an interesection of the hyperboloid with a plane containing two points ? I suspect so, but I do not see how to prove it.
Calculus of Variations
Using latitude, $\phi$ and longitude, $\theta$, we have
$$
1=\cos^2(\phi)\left(\frac{\mathrm{d}\theta}{\mathrm{d}s}\right)^2+\left(\frac{\mathrm{d}\phi}{\mathrm{d}s}\right)^2\tag1
$$
Integrating $(1)$ with respect to $s$, taking the first variation, and integrating by parts, we get
$$
\scriptsize\delta s=\int\left[\left(\color{#C00}{2\sin(\phi)\cos(\phi)\,\frac{\mathrm{d}\theta}{\mathrm{d}s}\,\frac{\mathrm{d}\phi}{\mathrm{d}s}-\cos^2(\phi)\,\frac{\mathrm{d}^2\theta}{\mathrm{d}s^2}}\right)\delta\theta-\left(\color{#090}{\sin(\phi)\cos(\phi)\,\left(\frac{\mathrm{d}\theta}{\mathrm{d}s}\right)^2+\frac{\mathrm{d}^2\phi}{\mathrm{d}s^2}}\right)\delta\phi\right]\mathrm{d}s\tag2
$$
Orthogonality requires that a curve where we have $\delta s=0$ for any $\delta\theta$ and $\delta\phi$ must satisfy
$$
\color{#C00}{\frac{\mathrm{d}^2\theta}{\mathrm{d}s^2}=2\tan(\phi)\,\frac{\mathrm{d}\theta}{\mathrm{d}s}\,\frac{\mathrm{d}\phi}{\mathrm{d}s}}\tag3
$$
and
$$
\color{#090}{\frac{\mathrm{d}^2\phi}{\mathrm{d}s^2}=-\sin(\phi)\cos(\phi)\,\left(\frac{\mathrm{d}\theta}{\mathrm{d}s}\right)^2}\tag4
$$
Solving the Equations
$$
\begin{align}
\log\left(\frac{\mathrm{d}\theta}{\mathrm{d}s}\right)
&=2\log(\sec(\phi))+\log(\cos(\epsilon))\tag{5a}\\
\frac{\mathrm{d}\theta}{\mathrm{d}s}
&=\cos(\epsilon)\sec^2(\phi)\tag{5b}\\
\left(\frac{\mathrm{d}\phi}{\mathrm{d}s}\right)^2
&=1-\cos^2(\epsilon)\sec^2(\phi)\tag{5c}\\
s
&=\int\frac{\mathrm{d}\phi}{\sqrt{1-\cos^2(\epsilon)\sec^2(\phi)}}\tag{5d}\\
&=\int\frac{\mathrm{d}\sin(\phi)}{\sqrt{\cos^2(\phi)-\cos^2(\epsilon)}}\tag{5e}\\
&=\int\frac{\mathrm{d}\sin(\phi)}{\sqrt{\sin^2(\epsilon)-\sin^2(\phi)}}\tag{5f}\\
&=\sin^{-1}\left(\frac{\sin(\phi)}{\sin(\epsilon)}\right)+s_0\tag{5g}\\[9pt]
\end{align}
$$
Explanation:
$\text{(5a)}$: divide $(3)$ by $\frac{\mathrm{d}\theta}{\mathrm{d}s}$ and integrate
$\text{(5b)}$: remove the log
$\text{(5c)}$: apply $(1)$ to $\text{(5b)}$
$\text{(5d)}$: separate the differential equation and integrate
$\text{(5e)}$: multiply the integrand by $\frac{\cos(\phi)}{\cos(\phi)}$
$\text{(5f)}$: $\cos^2(x)=1-\sin^2(x)$
$\text{(5g)}$: evaluate the integral
Solving $\text{(5g)}$ for $\sin(\phi)$ gives
$$
\bbox[5px,border:2px solid #C0A000]{\sin(\phi)=\sin(s-s_0)\sin(\epsilon)}\tag6
$$
Furthermore,
$$
\begin{align}
\theta
&=\int\frac{\cos(\epsilon)\,\mathrm{d}(s-s_0)}{1-\sin^2(s-s_0)\sin^2(\epsilon)}\tag{7a}\\
&=\int\frac{\cos(\epsilon)\,\mathrm{d}\tan(s-s_0)}{\sec^2(s-s_0)-\tan^2(s-s_0)\sin^2(\epsilon)}\tag{7b}\\
&=\int\frac{\cos(\epsilon)\,\mathrm{d}\tan(s-s_0)}{1+\tan^2(s-s_0)\cos^2(\epsilon)}\tag{7c}\\[6pt]
&=\tan^{-1}(\tan(s-s_0)\cos(\epsilon))+\theta_0\tag{7d}
\end{align}
$$
Explanation:
$\text{(7a)}$: apply $\text{(5h)}$ to $\text{(5b)}$, separate the differential equation, and integrate
$\text{(7b)}$: multiply the integrand by $\frac{\sec^2(s-s_0)}{\sec^2(s-s_0)}$
$\text{(7c)}$: $\sec^2(x)=1+\tan^2(x)$
$\text{(7d)}$: evaluate the integral
Solving $\text{(7d)}$ for $\tan(\theta-\theta_0)$ gives
$$
\bbox[5px,border:2px solid #C0A000]{\tan(\theta-\theta_0)=\tan(s-s_0)\cos(\epsilon)}\tag8
$$
Compare with a Great Circle
Given the parameterization in $\mathbb{R}^3$ of a great circle
$$
(x,y,z)=(\cos(s-s_0),\sin(s-s_0)\cos(\epsilon),\sin(s-s_0)\sin(\epsilon))\tag9
$$
we get
$$
\sin(\phi)=z=\sin(s-s_0)\sin(\epsilon)\tag{10}
$$
and
$$
\tan(\theta-\theta_0)=\frac yx=\tan(s-s_0)\cos(\epsilon)\tag{11}
$$
Equation $(10)$ matches equation $(6)$ and equation $(11)$ matches equation $(8)$. Thus, the shortest path between two points on a sphere is an arc of a great circle.
Plotting the Solved Equations
Since $\phi\in\left[-\frac\pi2,\frac\pi2\right]$, we can get $\phi$ by applying $\sin^{-1}$ to $(6)$:
$$
\bbox[5px,border:2px solid #C0A000]{\phi=\sin^{-1}(\sin(\epsilon)\sin(s-s_0))}\tag{12}
$$
However, to get $\theta-\theta_0$, we cannot simply apply $\tan^{-1}$ to $(8)$. Instead, we note that
$$
\begin{align}
\tan((\theta-\theta_0)-(s-s_0))
&=\frac{\tan(\theta-\theta_0)-\tan(s-s_0)}{1+\tan(\theta-\theta_0)\tan(s-s_0)}\tag{13a}\\
&=\frac{\tan(s-s_0)(\cos(\epsilon)-1)}{1+\cos(\epsilon)\tan^2(s-s_0)}\tag{13b}
\end{align}
$$
Explanation:
$\text{(13a)}$: formula for the tangent of a difference
$\text{(13b)}$: apply $(8)$
Thus, we get
$$
\bbox[5px,border:2px solid #C0A000]{\theta-\theta_0=(s-s_0)-\tan^{-1}\left(\frac{(1-\cos(\epsilon))\tan(s-s_0)}{1+\cos(\epsilon)\tan^2(s-s_0)}\right)}\tag{14}
$$
Great circles with $\epsilon=\left\{\color{#DF0000}{0^{\large\circ}},\color{#FF8F00}{15^{\large\circ}},\color{#EFDF00}{30^{\large\circ}},\color{#00BF00}{45^{\large\circ}},\color{#0000FF}{60^{\large\circ}},\color{#8000FF}{75^{\large\circ}},\color{#CC00FF}{90^{\large\circ}}\right\}$ plotted using $(12)$ and $(14)$.
Best Answer
Here is a low-tech approach which does not use derivatives much. Given distinct points $a,b$ on the sphere, let $\gamma$ be the shorter arc of great circle between them. Let $V$ be the surface obtained by rotating $\gamma$ about the line $ab$. Note that $V$ is convex (i.e., bounds a convex body) and lies inside the sphere. The nearest-point projection onto a convex set does not increase Euclidean distance and therefore does not increase length of curves (defined by means of sums of Euclidean distances over partitions). Therefore, it suffices to prove that a curve connecting $a$ to $b$ on surface $V$ cannot be shorter than $\gamma$.
Let $\Gamma$ be any such curve. Pick partition points $p_k$ on $\Gamma$ to approximate its length. For each $k$, draw the plane through $p_k$ that is orthogonal to line $ab$. This plane crosses $\gamma$ at a point $q_k$. It is easy to see that $|p_k-p_{k+1}|\ge |q_k-q_{k+1}|$ for each $k$. (Just write the Euclidean distance in cylindrical coordinates with $ab$ as the axis.) Since $\sum_k |q_k-q_{k+1}|$ is greater than or equal to the length of $\gamma$, the conclusion follows.