[proofplan]
The proof is the rigidity case of the Bonnet--Myers diameter comparison. We choose points $p,q \in M$ at maximal distance and study the two distance functions from $p$ and $q$. Laplacian comparison and the Calabi barrier maximum principle force the sum of these distance functions to be identically equal to the model value
\begin{align*}
\frac{\pi}{\sqrt{k}}.
\end{align*}
Equality in Laplacian comparison then forces equality in the [Riccati equation for geodesic spheres](/theorems/5357), so the metric is the spherical warped product with radial expression
\begin{align*}
dt^2+k^{-1}\sin^2(\sqrt{k}t)\,g_{S^{n-1}},
\end{align*}
which is the round metric of sectional curvature $k$.
[/proofplan]
[step:Choose antipodal points realizing the Bonnet--Myers diameter]
Let
\begin{align*}
d_g:M\times M&\to [0,\infty)
\end{align*}
denote the Riemannian distance function induced by $g$. Set
\begin{align*}
D := \frac{\pi}{\sqrt{k}}.
\end{align*}
By the [Bonnet--Myers theorem](/page/Bonnet--Myers%20Theorem), the Ricci lower bound $\operatorname{Ric}_g \geq (n-1)k\,g$ and the completeness of $(M,g)$ imply that $M$ is compact and
\begin{align*}
\operatorname{diam}(M,g) \leq D.
\end{align*}
The hypothesis gives equality. Since $M$ is compact, there exist points $p,q \in M$ such that
\begin{align*}
d_g(p,q)=D.
\end{align*}
Define the distance functions
\begin{align*}
r_p:M &\to [0,D], & x &\mapsto d_g(p,x),\\
r_q:M &\to [0,D], & x &\mapsto d_g(q,x).
\end{align*}
The triangle inequality gives, for every $x \in M$,
\begin{align*}
D=d_g(p,q)\leq d_g(p,x)+d_g(x,q)=r_p(x)+r_q(x).
\end{align*}
Thus the function
\begin{align*}
u:M &\to [0,\infty), & x &\mapsto r_p(x)+r_q(x)-D
\end{align*}
is non-negative.
[/step]
[step:Apply Laplacian comparison to make the distance sum superharmonic]
Let
\begin{align*}
\Omega := M \setminus \{p,q\}.
\end{align*}
On the [open set](/page/Open%20Set) where $r_p$ and $r_q$ are smooth, define the auxiliary angle functions
\begin{align*}
a:M\setminus\{p\} &\to (0,\pi], & x &\mapsto \sqrt{k}\,r_p(x),\\
b:M\setminus\{q\} &\to (0,\pi], & x &\mapsto \sqrt{k}\,r_q(x).
\end{align*}
Since $0<r_p(x)<D$ and $0<r_q(x)<D$ on $\Omega$, we have
\begin{align*}
0<a(x)<\pi, \qquad 0<b(x)<\pi.
\end{align*}
The inequality $r_p(x)+r_q(x)\geq D$ gives
\begin{align*}
a(x)+b(x)\geq \pi.
\end{align*}
Since also $a(x)+b(x)\leq 2\pi$, the identity
\begin{align*}
\cot a+\cot b=\frac{\sin(a+b)}{\sin a\,\sin b}
\end{align*}
implies
\begin{align*}
\cot a(x)+\cot b(x)\leq 0.
\end{align*}
The [Laplacian comparison theorem for distance functions](/page/Laplacian%20Comparison%20Theorem) under $\operatorname{Ric}_g\geq (n-1)k\,g$ gives, in the barrier sense across the cut loci and classically where the distance functions are smooth,
\begin{align*}
\Delta r_p(x) &\leq (n-1)\sqrt{k}\,\cot(\sqrt{k}\,r_p(x)),\\
\Delta r_q(x) &\leq (n-1)\sqrt{k}\,\cot(\sqrt{k}\,r_q(x)).
\end{align*}
Adding the two inequalities gives, in the barrier sense on $\Omega$,
\begin{align*}
\Delta u(x)
&=\Delta r_p(x)+\Delta r_q(x)\\
&\leq (n-1)\sqrt{k}\left(\cot a(x)+\cot b(x)\right)\\
&\leq 0.
\end{align*}
Thus $u$ is a non-negative superharmonic function on $\Omega$ in the barrier sense.
[guided]
The purpose of this step is to turn the geometric diameter condition into a maximum-principle statement. The function
\begin{align*}
u:M &\to [0,\infty), & x &\mapsto r_p(x)+r_q(x)-D
\end{align*}
measures the excess in the triangle inequality for the triple $(p,x,q)$. Because $d_g(p,q)=D$, the triangle inequality gives
\begin{align*}
u(x)=r_p(x)+r_q(x)-D\geq 0.
\end{align*}
To prove rigidity, we want to show that this excess is always zero.
We now compute the sign of $\Delta u$. On the set where both distance functions are smooth, define the maps
\begin{align*}
a:M\setminus\{p\} &\to (0,\pi], & x &\mapsto \sqrt{k}\,r_p(x),\\
b:M\setminus\{q\} &\to (0,\pi], & x &\mapsto \sqrt{k}\,r_q(x).
\end{align*}
The diameter bound gives $r_p(x)\leq D$ and $r_q(x)\leq D$, while $x \neq p,q$ gives $r_p(x)>0$ and $r_q(x)>0$. Hence
\begin{align*}
0<a(x)<\pi, \qquad 0<b(x)<\pi.
\end{align*}
The triangle inequality gives
\begin{align*}
a(x)+b(x)=\sqrt{k}\,(r_p(x)+r_q(x))\geq \sqrt{k}\,D=\pi.
\end{align*}
Since each of $a(x)$ and $b(x)$ is at most $\pi$, we also have $a(x)+b(x)\leq 2\pi$. Therefore
\begin{align*}
\sin(a(x)+b(x))\leq 0.
\end{align*}
Because $\sin a(x)>0$ and $\sin b(x)>0$, the trigonometric identity
\begin{align*}
\cot a+\cot b=\frac{\sin(a+b)}{\sin a\,\sin b}
\end{align*}
gives
\begin{align*}
\cot a(x)+\cot b(x)\leq 0.
\end{align*}
Now apply the [Laplacian comparison theorem for distance functions](/page/Laplacian%20Comparison%20Theorem) under the Ricci lower bound $\operatorname{Ric}_g\geq (n-1)k\,g$. The theorem applies to $r_p$ and $r_q$ because $(M,g)$ is complete and has the stated Ricci curvature lower bound. It gives, in the barrier sense and classically away from the cut loci,
\begin{align*}
\Delta r_p(x) &\leq (n-1)\sqrt{k}\,\cot(\sqrt{k}\,r_p(x)),\\
\Delta r_q(x) &\leq (n-1)\sqrt{k}\,\cot(\sqrt{k}\,r_q(x)).
\end{align*}
Adding these inequalities yields
\begin{align*}
\Delta u(x)
&=\Delta r_p(x)+\Delta r_q(x)\\
&\leq (n-1)\sqrt{k}\left(\cot a(x)+\cot b(x)\right)\\
&\leq 0.
\end{align*}
Thus the triangle-excess function $u$ is superharmonic in the barrier sense on $\Omega=M\setminus\{p,q\}$. This is the key analytic input: a non-negative superharmonic function that attains an interior minimum must be constant.
[/guided]
[/step]
[step:Use the Calabi maximum principle to force equality in the triangle inequality]
Let
\begin{align*}
\gamma:[0,D] &\to M
\end{align*}
be a unit-speed minimizing geodesic from $p$ to $q$, so $\gamma(0)=p$ and $\gamma(D)=q$. For every $t \in (0,D)$,
\begin{align*}
r_p(\gamma(t))=t,\qquad r_q(\gamma(t))=D-t,
\end{align*}
and hence
\begin{align*}
u(\gamma(t))=0.
\end{align*}
Since $n\geq 2$, the punctured manifold $\Omega=M\setminus\{p,q\}$ is connected. The [Calabi barrier maximum principle for superharmonic functions](/page/Calabi%20Maximum%20Principle) applies because $u$ is continuous on $\Omega$, non-negative on $\Omega$, superharmonic in the barrier sense on $\Omega$, and attains its minimum value $0$ at the interior point $\gamma(t)$ for every $t\in(0,D)$. The principle gives
\begin{align*}
u \equiv 0 \quad \text{on } \Omega.
\end{align*}
By continuity of $r_p$ and $r_q$, the same identity holds on all of $M$:
\begin{align*}
r_p(x)+r_q(x)=D
\end{align*}
for every $x \in M$.
[/step]
[step:Extract equality in the Laplacian comparison inequalities]
On the regular set where both $r_p$ and $r_q$ are smooth, the identity $r_p+r_q=D$ gives
\begin{align*}
\Delta r_p+\Delta r_q=0.
\end{align*}
The comparison inequalities and the trigonometric identity give
\begin{align*}
0
&=\Delta r_p+\Delta r_q\\
&\leq (n-1)\sqrt{k}\left(\cot(\sqrt{k}\,r_p)+\cot(\sqrt{k}\,r_q)\right)\\
&=(n-1)\sqrt{k}\,
\frac{\sin(\sqrt{k}(r_p+r_q))}
{\sin(\sqrt{k}\,r_p)\sin(\sqrt{k}\,r_q)}\\
&=0.
\end{align*}
Therefore equality holds throughout this chain. In particular,
\begin{align*}
\Delta r_p=(n-1)\sqrt{k}\,\cot(\sqrt{k}\,r_p)
\end{align*}
where $r_p$ is smooth.
[claim:Every point between $p$ and $q$ is regular for the distance from $p$]
Every point $x\in M\setminus\{p,q\}$ is not in the cut locus of $p$. Consequently $r_p$ is smooth on $M\setminus\{p,q\}$, and the displayed equality for $\Delta r_p$ holds on all of $M\setminus\{p,q\}$.
[/claim]
[proof]
Fix $x\in M\setminus\{p,q\}$. Let
\begin{align*}
\alpha:[0,r_p(x)]&\to M, & s&\mapsto \alpha(s)
\end{align*}
be a unit-speed minimizing geodesic from $p$ to $x$, and let
\begin{align*}
\beta:[0,r_q(x)]&\to M, & s&\mapsto \beta(s)
\end{align*}
be a unit-speed minimizing geodesic from $x$ to $q$. Since $r_p(x)+r_q(x)=D=d_g(p,q)$, the concatenation of $\alpha$ and $\beta$ has length $D$ and is therefore a minimizing curve from $p$ to $q$. A length-minimizing curve is a geodesic on each interior subinterval and has no corner at the joining point, so the terminal velocity of $\alpha$ equals the initial velocity of $\beta$ at $x$.
This proves uniqueness of the minimizing geodesic from $p$ to $x$: if $\alpha_1$ and $\alpha_2$ are two such geodesics, then concatenating each with the same $\beta$ gives two minimizing geodesics from $p$ to $q$ that agree on the non-empty interval traced by $\beta$ after $x$. By uniqueness for the geodesic initial-value problem, the two geodesics agree after reversing time from $x$ to $p$, hence $\alpha_1=\alpha_2$.
It remains to exclude conjugacy. If $x$ were conjugate to $p$ along the unique minimizing geodesic from $p$ to $x$, then $x$ would be an interior conjugate point along the minimizing geodesic from $p$ to $q$ obtained by continuing through $x$ to $q$. By the [second variation theorem for geodesics](/page/Second%20Variation%20Formula), an interior conjugate point produces a non-zero variation field along the geodesic vanishing at the endpoints and with non-positive index form. Perturbing this field by a compactly supported field near the conjugate point gives a proper variation through curves from $p$ to $q$ with strictly negative second variation of length; hence the original geodesic cannot be length-minimizing between $p$ and $q$. This contradicts $d_g(p,q)=D$. Therefore $x$ is neither a branching cut point nor a conjugate cut point of $p$, so $x$ is not in the cut locus of $p$.
[/proof]
[/step]
[step:Convert equality in Laplacian comparison into the model shape operator]
Let
\begin{align*}
S_pM:=\{v\in T_pM:g_p(v,v)=1\}
\end{align*}
be the unit tangent sphere at $p$, where $g_p:T_pM\times T_pM\to \mathbb{R}$ is the [inner product](/page/Inner%20Product) induced by $g$ on $T_pM$. Let
\begin{align*}
\exp_p:\mathcal{D}_p\subset T_pM&\to M
\end{align*}
denote the Riemannian exponential map at $p$, defined on its maximal star-shaped domain $\mathcal{D}_p$. For each $v\in S_pM$, define the radial geodesic
\begin{align*}
\gamma_v:(0,D)&\to M, & t&\mapsto \exp_p(tv).
\end{align*}
For $t\in(0,D)$, let $A_t:T_{\gamma_v(t)}\{r_p=t\}\to T_{\gamma_v(t)}\{r_p=t\}$ denote the shape operator of the geodesic sphere $\{r_p=t\}$ with respect to the outward unit normal $\nabla r_p$.
The [Riccati equation for the shape operator of geodesic spheres](/page/Riccati%20Equation%20For%20Geodesic%20Spheres) gives
\begin{align*}
\frac{d}{dt}\operatorname{tr}A_t+\operatorname{tr}(A_t^2)+\operatorname{Ric}_g(\dot{\gamma}_v(t),\dot{\gamma}_v(t))=0.
\end{align*}
Since
\begin{align*}
\operatorname{tr}A_t=\Delta r_p(\gamma_v(t))=(n-1)\sqrt{k}\,\cot(\sqrt{k}t),
\end{align*}
we compute
\begin{align*}
\frac{d}{dt}\operatorname{tr}A_t
=-(n-1)k\,\csc^2(\sqrt{k}t).
\end{align*}
The Cauchy--Schwarz inequality for the eigenvalues of the symmetric endomorphism $A_t$ gives
\begin{align*}
\operatorname{tr}(A_t^2)\geq \frac{(\operatorname{tr}A_t)^2}{n-1}
=(n-1)k\,\cot^2(\sqrt{k}t).
\end{align*}
The Ricci lower bound gives
\begin{align*}
\operatorname{Ric}_g(\dot{\gamma}_v(t),\dot{\gamma}_v(t))\geq (n-1)k.
\end{align*}
Substituting these estimates into the Riccati equation yields equality in both estimates, because
\begin{align*}
-(n-1)k\,\csc^2(\sqrt{k}t)
+(n-1)k\,\cot^2(\sqrt{k}t)
+(n-1)k=0.
\end{align*}
Equality in the Cauchy--Schwarz inequality for the eigenvalues of $A_t$ means all principal curvatures are equal. Therefore
\begin{align*}
A_t=\sqrt{k}\,\cot(\sqrt{k}t)\,\operatorname{id}_{T_{\gamma_v(t)}\{r_p=t\}}
\end{align*}
for every $t\in(0,D)$ and every $v\in S_pM$.
[/step]
[step:Integrate the shape operator equation to obtain the spherical warped product metric]
Define the polar-coordinate map
\begin{align*}
F:(0,D)\times S_pM &\to M\setminus\{p,q\}, & (t,v)&\mapsto \exp_p(tv).
\end{align*}
The preceding claim implies that $\exp_p$ has no cut point at distance $t\in(0,D)$. Hence $F$ is locally a diffeomorphism everywhere. It is surjective because every $x\in M\setminus\{p,q\}$ lies on the unique minimizing geodesic from $p$ to $x$, and it is injective because that minimizing geodesic is unique. Therefore $F$ is a diffeomorphism. Write
\begin{align*}
F^*g=dt^2+g_t,
\end{align*}
where $g_t$ is the Riemannian metric induced on the slice $\{t\}\times S_pM$.
The relation between the slice metric and the shape operator is
\begin{align*}
\frac{1}{2}\partial_t g_t(X,Y)=g_t(A_tX,Y)
\end{align*}
for smooth vector fields $X,Y$ tangent to $S_pM$ and independent of $t$. Using
\begin{align*}
A_t=\sqrt{k}\,\cot(\sqrt{k}t)\operatorname{id},
\end{align*}
we get
\begin{align*}
\partial_t g_t(X,Y)=2\sqrt{k}\,\cot(\sqrt{k}t)\,g_t(X,Y).
\end{align*}
Solving this scalar differential equation in $t$ gives
\begin{align*}
g_t=\frac{\sin^2(\sqrt{k}t)}{\sin^2(\sqrt{k}t_0)}\,g_{t_0}
\end{align*}
for every fixed $t_0\in(0,D)$.
In geodesic [normal coordinates](/theorems/2713) at $p$, the slice metrics satisfy
\begin{align*}
\lim_{t\downarrow 0}\frac{1}{t^2}g_t=g_{S_pM},
\end{align*}
where $g_{S_pM}$ is the standard round metric on the unit sphere in the Euclidean [inner product space](/page/Inner%20Product%20Space) $(T_pM,g_p)$. Since
\begin{align*}
\lim_{t\downarrow 0}\frac{\sin^2(\sqrt{k}t)}{t^2}=k,
\end{align*}
the preceding identity gives
\begin{align*}
g_t=\frac{1}{k}\sin^2(\sqrt{k}t)\,g_{S_pM}.
\end{align*}
Consequently
\begin{align*}
F^*g=dt^2+\frac{1}{k}\sin^2(\sqrt{k}t)\,g_{S_pM}
\end{align*}
on $(0,D)\times S_pM$.
[/step]
[step:Identify the warped product with the round sphere of curvature $k$]
Let $g_{\mathrm{round}}$ denote the unit round metric on $S^n$. The round metric of constant sectional curvature $k$ is $k^{-1}g_{\mathrm{round}}$. In geodesic polar coordinates around a pole, this metric has the form
\begin{align*}
dt^2+\frac{1}{k}\sin^2(\sqrt{k}t)\,g_{S^{n-1}},
\qquad 0<t<\frac{\pi}{\sqrt{k}}.
\end{align*}
Since $(S_pM,g_{S_pM})$ is isometric to the unit round sphere $S^{n-1}$, the formula obtained above identifies $(M\setminus\{p,q\},g)$ with the complement of the two poles in the round sphere of sectional curvature $k$.
The metric expression extends smoothly across $t=0$ because
\begin{align*}
\frac{1}{k}\sin^2(\sqrt{k}t)=t^2+O(t^4)
\end{align*}
as $t\downarrow 0$, which is exactly the first-order radial behaviour required in geodesic normal coordinates at $p$. At the other endpoint, put $s:=D-t$. Since $\sqrt{k}D=\pi$, we have
\begin{align*}
\frac{1}{k}\sin^2(\sqrt{k}t)=\frac{1}{k}\sin^2(\sqrt{k}s)=s^2+O(s^4)
\end{align*}
as $s\downarrow 0$. Together with the identity $r_p+r_q=D$, this is the corresponding geodesic normal-coordinate behaviour at $q$. Therefore the isometry on the punctured manifolds extends continuously and smoothly across the two endpoints. Hence $(M,g)$ is isometric to the round sphere $S^n$ with constant sectional curvature $k$.
[/step]