[proofplan]
We prove the identity separately in the three constant-curvature models. In the Euclidean case, we place the vertex $\bar{p}$ at the origin and compute the squared distance between the two endpoint vectors. In positive curvature, we use the round sphere of radius $1/\sqrt{k}$ embedded in $\mathbb{R}^3$, where geodesic distance is encoded by the ambient Euclidean [inner product](/page/Inner%20Product). In negative curvature, we use the hyperboloid model in Minkowski space, where geodesic distance is encoded by the Lorentzian inner product; evaluating the two endpoint vectors gives the hyperbolic cosine law.
[/proofplan]
[step:Compute the Euclidean case by expanding the endpoint difference]
Assume $k = 0$. Then $M_0^2$ is isometric to $\mathbb{R}^2$ with the Euclidean metric. Use an isometry to place $\bar{p}$ at $0 \in \mathbb{R}^2$, and choose the standard [orthonormal basis](/page/Orthonormal%20Basis) vectors
\begin{align*}
e_1 &= (1,0),&
e_2 &= (0,1).
\end{align*}
Since $c = d(\bar{p},\bar{q})$ and $b = d(\bar{p},\bar{r})$, and the angle at $\bar{p}$ is $\alpha$, we may write
\begin{align*}
\bar{q} &= c e_1,\\
\bar{r} &= b(\cos \alpha\, e_1 + \sin \alpha\, e_2).
\end{align*}
Therefore
\begin{align*}
a^2
&= |\bar{r}-\bar{q}|^2\\
&= |b(\cos \alpha\, e_1 + \sin \alpha\, e_2) - c e_1|^2\\
&= |(b\cos\alpha-c)e_1 + b\sin\alpha\,e_2|^2\\
&= (b\cos\alpha-c)^2 + b^2\sin^2\alpha\\
&= b^2\cos^2\alpha - 2bc\cos\alpha + c^2 + b^2\sin^2\alpha\\
&= b^2 + c^2 - 2bc\cos\alpha.
\end{align*}
This proves the stated formula when $k=0$.
[/step]
[step:Represent the positive curvature model inside Euclidean space]
Assume $k>0$, and define the radius
\begin{align*}
\rho := \frac{1}{\sqrt{k}}.
\end{align*}
Realize $M_k^2$ as the round sphere
\begin{align*}
S_\rho^2 := \{x \in \mathbb{R}^3 : x \cdot x = \rho^2\}
\end{align*}
with the metric induced by the Euclidean inner product $x \cdot y$ on $\mathbb{R}^3$. The geodesic distance on $S_\rho^2$ satisfies
\begin{align*}
x \cdot y = \rho^2 \cos\left(\frac{d(x,y)}{\rho}\right)
\end{align*}
for any two points $x,y \in S_\rho^2$ joined by a minimizing geodesic segment.
Use an isometry of $S_\rho^2$ to place
\begin{align*}
\bar{p} = (\rho,0,0).
\end{align*}
Let
\begin{align*}
u &:= (0,1,0),&
v &:= (0,\cos\alpha,\sin\alpha).
\end{align*}
Then $u,v \in T_{\bar{p}}S_\rho^2$, $|u|=|v|=1$, and $u \cdot v = \cos\alpha$. The unit-speed geodesic on $S_\rho^2$ starting at $\bar{p}$ with initial velocity $w \in T_{\bar{p}}S_\rho^2$ and $|w|=1$ is the map
\begin{align*}
\gamma_w: \mathbb{R} &\to S_\rho^2\\
t &\mapsto \cos\left(\frac{t}{\rho}\right)\bar{p}
+ \rho \sin\left(\frac{t}{\rho}\right)w.
\end{align*}
Thus, choosing the geodesic to $\bar{q}$ in direction $u$ and the geodesic to $\bar{r}$ in direction $v$, we have
\begin{align*}
\bar{q} &= \cos\left(\frac{c}{\rho}\right)\bar{p}
+ \rho \sin\left(\frac{c}{\rho}\right)u,\\
\bar{r} &= \cos\left(\frac{b}{\rho}\right)\bar{p}
+ \rho \sin\left(\frac{b}{\rho}\right)v.
\end{align*}
[guided]
The useful feature of the sphere model is that geodesic distance can be recovered from the ordinary Euclidean inner product in $\mathbb{R}^3$. We set
\begin{align*}
\rho := \frac{1}{\sqrt{k}},
\end{align*}
so the model space of curvature $k$ is the round sphere
\begin{align*}
S_\rho^2 := \{x \in \mathbb{R}^3 : x \cdot x = \rho^2\}.
\end{align*}
For two points $x,y \in S_\rho^2$ connected by a minimizing geodesic segment, the central angle between the radius vectors $x$ and $y$ is $d(x,y)/\rho$, hence
\begin{align*}
x \cdot y = \rho^2 \cos\left(\frac{d(x,y)}{\rho}\right).
\end{align*}
We may rotate the sphere without changing distances or angles, so place
\begin{align*}
\bar{p} = (\rho,0,0).
\end{align*}
The tangent plane at $\bar{p}$ is the plane perpendicular to $\bar{p}$, namely
\begin{align*}
T_{\bar{p}}S_\rho^2 = \{0\}\times \mathbb{R}^2.
\end{align*}
Choose unit tangent vectors
\begin{align*}
u &:= (0,1,0),&
v &:= (0,\cos\alpha,\sin\alpha).
\end{align*}
They satisfy
\begin{align*}
u \cdot v = \cos\alpha,
\end{align*}
so the angle between them is exactly $\alpha$.
The unit-speed geodesic starting at $\bar{p}$ in a unit tangent direction $w \in T_{\bar{p}}S_\rho^2$ is
\begin{align*}
\gamma_w: \mathbb{R} &\to S_\rho^2\\
t &\mapsto \cos\left(\frac{t}{\rho}\right)\bar{p}
+ \rho \sin\left(\frac{t}{\rho}\right)w.
\end{align*}
This formula stays on $S_\rho^2$ because $\bar{p}\cdot w=0$, $|\bar{p}|=\rho$, and $|w|=1$. It has speed one because differentiating gives
\begin{align*}
\gamma_w'(t)
= -\frac{1}{\rho}\sin\left(\frac{t}{\rho}\right)\bar{p}
+ \cos\left(\frac{t}{\rho}\right)w,
\end{align*}
and therefore
\begin{align*}
|\gamma_w'(t)|^2
= \sin^2\left(\frac{t}{\rho}\right)
+ \cos^2\left(\frac{t}{\rho}\right)
= 1.
\end{align*}
Hence the points at distances $c$ and $b$ from $\bar{p}$ in the two chosen directions are
\begin{align*}
\bar{q} &= \cos\left(\frac{c}{\rho}\right)\bar{p}
+ \rho \sin\left(\frac{c}{\rho}\right)u,\\
\bar{r} &= \cos\left(\frac{b}{\rho}\right)\bar{p}
+ \rho \sin\left(\frac{b}{\rho}\right)v.
\end{align*}
[/guided]
[/step]
[step:Evaluate the Euclidean inner product to obtain the spherical formula]
Using $\bar{p}\cdot u=\bar{p}\cdot v=0$, $\bar{p}\cdot\bar{p}=\rho^2$, and $u\cdot v=\cos\alpha$, we compute
\begin{align*}
\bar{q}\cdot\bar{r}
&= \rho^2\cos\left(\frac{c}{\rho}\right)\cos\left(\frac{b}{\rho}\right)
+ \rho^2\sin\left(\frac{c}{\rho}\right)\sin\left(\frac{b}{\rho}\right)\cos\alpha.
\end{align*}
On the other hand, since $a=d(\bar{q},\bar{r})$,
\begin{align*}
\bar{q}\cdot\bar{r}
= \rho^2\cos\left(\frac{a}{\rho}\right).
\end{align*}
Dividing by $\rho^2>0$ gives
\begin{align*}
\cos\left(\frac{a}{\rho}\right)
=
\cos\left(\frac{b}{\rho}\right)\cos\left(\frac{c}{\rho}\right)
+
\sin\left(\frac{b}{\rho}\right)\sin\left(\frac{c}{\rho}\right)\cos\alpha.
\end{align*}
Since $1/\rho=\sqrt{k}$, this is
\begin{align*}
\cos(\sqrt{k}\,a)
=
\cos(\sqrt{k}\,b)\cos(\sqrt{k}\,c)
+
\sin(\sqrt{k}\,b)\sin(\sqrt{k}\,c)\cos\alpha.
\end{align*}
This proves the stated formula when $k>0$.
[/step]
[step:Represent the negative curvature model inside Minkowski space]
Assume $k<0$, and define
\begin{align*}
\rho := \frac{1}{\sqrt{-k}}.
\end{align*}
Let $\mathbb{R}^{2,1}$ denote $\mathbb{R}^3$ equipped with the Lorentzian [bilinear form](/page/Bilinear%20Form)
\begin{align*}
\langle x,y\rangle_L := -x_0y_0 + x_1y_1 + x_2y_2
\end{align*}
for $x=(x_0,x_1,x_2)$ and $y=(y_0,y_1,y_2)$. Realize $M_k^2$ as the upper hyperboloid
\begin{align*}
H_\rho^2 := \{x \in \mathbb{R}^{2,1} : \langle x,x\rangle_L = -\rho^2,\ x_0>0\}
\end{align*}
with its induced Riemannian metric. The distance on $H_\rho^2$ satisfies
\begin{align*}
-\langle x,y\rangle_L = \rho^2\cosh\left(\frac{d(x,y)}{\rho}\right)
\end{align*}
for all $x,y \in H_\rho^2$.
Use an isometry of $H_\rho^2$ to place
\begin{align*}
\bar{p} = (\rho,0,0).
\end{align*}
Let
\begin{align*}
u &:= (0,1,0),&
v &:= (0,\cos\alpha,\sin\alpha).
\end{align*}
Then $u,v \in T_{\bar{p}}H_\rho^2$, $\langle u,u\rangle_L=\langle v,v\rangle_L=1$, and $\langle u,v\rangle_L=\cos\alpha$. The unit-speed geodesic on $H_\rho^2$ starting at $\bar{p}$ with initial velocity $w \in T_{\bar{p}}H_\rho^2$ and $\langle w,w\rangle_L=1$ is the map
\begin{align*}
\eta_w: \mathbb{R} &\to H_\rho^2\\
t &\mapsto \cosh\left(\frac{t}{\rho}\right)\bar{p}
+ \rho\sinh\left(\frac{t}{\rho}\right)w.
\end{align*}
Thus
\begin{align*}
\bar{q} &= \cosh\left(\frac{c}{\rho}\right)\bar{p}
+ \rho\sinh\left(\frac{c}{\rho}\right)u,\\
\bar{r} &= \cosh\left(\frac{b}{\rho}\right)\bar{p}
+ \rho\sinh\left(\frac{b}{\rho}\right)v.
\end{align*}
[/step]
[step:Evaluate the Lorentzian inner product to obtain the hyperbolic formula]
Using $\langle \bar{p},u\rangle_L=\langle \bar{p},v\rangle_L=0$, $\langle \bar{p},\bar{p}\rangle_L=-\rho^2$, and $\langle u,v\rangle_L=\cos\alpha$, we compute
\begin{align*}
\langle \bar{q},\bar{r}\rangle_L
&=
-\rho^2\cosh\left(\frac{c}{\rho}\right)\cosh\left(\frac{b}{\rho}\right)
+
\rho^2\sinh\left(\frac{c}{\rho}\right)\sinh\left(\frac{b}{\rho}\right)\cos\alpha.
\end{align*}
Since $a=d(\bar{q},\bar{r})$, the hyperboloid distance relation gives
\begin{align*}
-\langle \bar{q},\bar{r}\rangle_L
=
\rho^2\cosh\left(\frac{a}{\rho}\right).
\end{align*}
Combining the last two identities and dividing by $\rho^2>0$ yields
\begin{align*}
\cosh\left(\frac{a}{\rho}\right)
=
\cosh\left(\frac{b}{\rho}\right)\cosh\left(\frac{c}{\rho}\right)
-
\sinh\left(\frac{b}{\rho}\right)\sinh\left(\frac{c}{\rho}\right)\cos\alpha.
\end{align*}
Since $1/\rho=\sqrt{-k}$, this is
\begin{align*}
\cosh(\sqrt{-k}\,a)
=
\cosh(\sqrt{-k}\,b)\cosh(\sqrt{-k}\,c)
-
\sinh(\sqrt{-k}\,b)\sinh(\sqrt{-k}\,c)\cos\alpha.
\end{align*}
This proves the stated formula when $k<0$. Together with the Euclidean and spherical cases, the theorem follows for every $k \in \mathbb{R}$.
[/step]