[proofplan]
We prove the operator identity by testing it on Jacobi fields coming from variations through radial geodesics from $p$. Since $\gamma(t)$ stays inside the normal domain, these variation Jacobi fields span the tangent space $E_t$ to the geodesic sphere at $\gamma(t)$. The shape operator sends each such Jacobi field $J(t)$ to its covariant derivative $\nabla_{\dot{\gamma}}J(t)$, and differentiating this identity converts the Jacobi equation into the Riccati equation.
[/proofplan]
[step:Construct radial Jacobi fields spanning the tangent spaces of the geodesic spheres]
Fix $t_0 \in (0,a)$. Let $v_0 := \dot{\gamma}(0) \in T_pM$, so $|v_0|_g=1$ and $\gamma(t)=\exp_p(tv_0)$ for $0 \le t<a$. Since $\gamma(t_0)\in\Omega_p$, the exponential map is nonsingular at $t_0v_0$ on the radial normal domain. Hence the differential
\begin{align*}
d(\exp_p)_{t_0v_0}: T_{t_0v_0}(T_pM) \cong T_pM \to T_{\gamma(t_0)}M
\end{align*}
maps the hyperplane $v_0^\perp \subset T_pM$ isomorphically onto
\begin{align*}
E_{t_0}=T_{\gamma(t_0)}S_{t_0}.
\end{align*}
For each $w \in v_0^\perp$, choose a smooth curve
\begin{align*}
\sigma: (-\varepsilon,\varepsilon) \to T_pM
\end{align*}
with $\sigma(0)=v_0$, $\sigma'(0)=w$, and $|\sigma(s)|_g=1$ for all $s$. Define the radial geodesic variation
\begin{align*}
F: (-\varepsilon,\varepsilon)\times [0,a) &\to M \\
(s,t) &\mapsto \exp_p(t\sigma(s)).
\end{align*}
Let
\begin{align*}
J_w: [0,a) &\to TM \\
t &\mapsto \frac{\partial F}{\partial s}(0,t)
\end{align*}
be the associated variation field along $\gamma$. Then $J_w$ is a perpendicular Jacobi field with $J_w(0)=0$, and
\begin{align*}
J_w(t_0)=d(\exp_p)_{t_0v_0}(t_0w).
\end{align*}
As $w$ ranges over $v_0^\perp$, the vectors $J_w(t_0)$ span $E_{t_0}$.
[guided]
The purpose of this step is to produce enough Jacobi fields to test an operator identity on the whole space $E_{t_0}$. We start at the initial point $p$ and write $v_0:=\dot{\gamma}(0)$. Since $\gamma$ is unit-speed and radial from $p$,
\begin{align*}
\gamma(t)=\exp_p(tv_0)
\end{align*}
for $0\le t<a$.
Because $\gamma(t_0)$ lies inside the normal domain $\Omega_p$, the exponential map is nonsingular at $t_0v_0$. The geodesic sphere $S_{t_0}$ is the image under $\exp_p$ of the sphere of radius $t_0$ in $T_pM$, and its tangent space at $t_0v_0$ is the hyperplane $v_0^\perp$. Therefore the differential
\begin{align*}
d(\exp_p)_{t_0v_0}: T_{t_0v_0}(T_pM)\cong T_pM \to T_{\gamma(t_0)}M
\end{align*}
maps $v_0^\perp$ isomorphically onto
\begin{align*}
E_{t_0}=T_{\gamma(t_0)}S_{t_0}.
\end{align*}
Now fix $w\in v_0^\perp$. Choose a smooth curve
\begin{align*}
\sigma: (-\varepsilon,\varepsilon)\to T_pM
\end{align*}
such that $\sigma(0)=v_0$, $\sigma'(0)=w$, and $|\sigma(s)|_g=1$ for all $s$. The condition $\sigma'(0)\in v_0^\perp$ is exactly the first-order condition for staying on the unit sphere in $T_pM$. Define
\begin{align*}
F: (-\varepsilon,\varepsilon)\times [0,a) &\to M \\
(s,t) &\mapsto \exp_p(t\sigma(s)).
\end{align*}
For each fixed $s$, the map $t\mapsto F(s,t)$ is a unit-speed geodesic starting at $p$. Thus the variation field
\begin{align*}
J_w: [0,a) &\to TM \\
t &\mapsto \frac{\partial F}{\partial s}(0,t)
\end{align*}
is a Jacobi field along $\gamma$ by the Jacobi equation for variation fields (citing a result not yet in the wiki: Jacobi equation for geodesic variation fields). Since all geodesics in the variation start at $p$, we have
\begin{align*}
J_w(0)=\frac{\partial F}{\partial s}(0,0)=0.
\end{align*}
Differentiating $F(s,t_0)=\exp_p(t_0\sigma(s))$ at $s=0$ gives
\begin{align*}
J_w(t_0)=d(\exp_p)_{t_0v_0}(t_0w).
\end{align*}
Because $d(\exp_p)_{t_0v_0}$ maps $v_0^\perp$ isomorphically onto $E_{t_0}$, the vectors $J_w(t_0)$ span $E_{t_0}$ as $w$ ranges over $v_0^\perp$.
[/guided]
[/step]
[step:Identify the shape operator with differentiation of radial Jacobi fields]
Let $w\in v_0^\perp$, and let $J:=J_w$ be the radial Jacobi field constructed above. For $0<t<a$, the vector $J(t)$ is tangent to the geodesic sphere $S_t$, hence $J(t)\in E_t$. Since $F$ is a radial variation, its coordinate vector fields satisfy
\begin{align*}
\frac{\partial F}{\partial t}(0,t)=\dot{\gamma}(t)=\nabla r_{\gamma(t)}
\end{align*}
and
\begin{align*}
\frac{\partial F}{\partial s}(0,t)=J(t).
\end{align*}
The Levi-Civita connection is torsion-free, so the coordinate vector fields of the two-parameter map $F$ satisfy
\begin{align*}
\nabla_{\partial_t F}\partial_s F=\nabla_{\partial_s F}\partial_t F.
\end{align*}
Evaluating at $(0,t)$ gives
\begin{align*}
\nabla_{\dot{\gamma}(t)}J(t)
=
\nabla_{J(t)}\nabla r.
\end{align*}
By the definition of $A(t)$, this is
\begin{align*}
\nabla_{\dot{\gamma}(t)}J(t)=A(t)J(t).
\end{align*}
Writing
\begin{align*}
J'(t):=\nabla_{\dot{\gamma}(t)}J(t),
\end{align*}
we have
\begin{align*}
J'(t)=A(t)J(t)
\end{align*}
for every $0<t<a$.
[/step]
[step:Differentiate the shape operator identity along the geodesic]
For the same Jacobi field $J$, use the definition of the covariant derivative $A'(t)$ of the endomorphism family $A(t)$. Since $J'(t)=A(t)J(t)$, we compute
\begin{align*}
J''(t)
&=\nabla_{\dot{\gamma}(t)}J'(t) \\
&=\nabla_{\dot{\gamma}(t)}(A(t)J(t)) \\
&=A'(t)J(t)+A(t)J'(t) \\
&=A'(t)J(t)+A(t)^2J(t).
\end{align*}
Thus
\begin{align*}
J''(t)=\bigl(A'(t)+A(t)^2\bigr)J(t).
\end{align*}
[/step]
[step:Use the Jacobi equation to obtain the Riccati equation on each spanning vector]
The Jacobi equation for $J$ along $\gamma$ is
\begin{align*}
J''(t)+R(J(t),\dot{\gamma}(t))\dot{\gamma}(t)=0
\end{align*}
for $0<t<a$ (citing a result not yet in the wiki: Jacobi equation for geodesic variation fields). By the definition of $R_\gamma(t)$, this becomes
\begin{align*}
J''(t)+R_\gamma(t)J(t)=0.
\end{align*}
Substituting the identity from the previous step gives
\begin{align*}
\bigl(A'(t)+A(t)^2+R_\gamma(t)\bigr)J(t)=0.
\end{align*}
Hence the endomorphism $A'(t)+A(t)^2+R_\gamma(t)$ vanishes on every vector of the form $J_w(t)$ arising from a radial Jacobi field.
[/step]
[step:Conclude the operator identity from the spanning property]
Fix $t_0\in(0,a)$. By the spanning result in the first step, every vector $X\in E_{t_0}$ has the form
\begin{align*}
X=J_w(t_0)
\end{align*}
for some $w\in v_0^\perp$. Applying the previous step at $t=t_0$ gives
\begin{align*}
\bigl(A'(t_0)+A(t_0)^2+R_\gamma(t_0)\bigr)X=0.
\end{align*}
Since $X\in E_{t_0}$ was arbitrary, the endomorphism
\begin{align*}
A'(t_0)+A(t_0)^2+R_\gamma(t_0):E_{t_0}\to E_{t_0}
\end{align*}
is the zero map. Since $t_0\in(0,a)$ was arbitrary, we have
\begin{align*}
A'(t)+A(t)^2+R_\gamma(t)=0
\end{align*}
for every $0<t<a$, as claimed.
[/step]
