[proofplan]
Shrink the neighbourhood only so that the signed-distance description of $C$ and the Hessian hypothesis are valid throughout the region where the geodesic is assumed to lie. For a unit-speed minimizing geodesic $\gamma$ joining two points of $C \cap U$, consider the one-variable function $f=\rho\circ\gamma$. The geodesic equation converts the second derivative of $f$ into the Hessian of $\rho$ along $\dot{\gamma}$, so the assumed nonnegativity makes $f$ convex. Since $f$ is nonpositive at both endpoints, convexity forces $f\le 0$ along the entire segment, which is exactly the statement that $\gamma$ remains in $C$.
[/proofplan]
[step:Choose a neighbourhood where the signed distance defines the region]
Since $\partial C$ is $C^2$ near $x_0$, the signed-distance function is $C^2$ on some neighbourhood of $x_0$ after restricting to a sufficiently small tubular neighbourhood of $\partial C$. Replace $U_0$ by a smaller open neighbourhood, still denoted $U_0$, such that
\begin{align*}
C \cap U_0 = \{x \in U_0 : \rho(x) \le 0\}.
\end{align*}
Let $U \subset U_0$ be any open neighbourhood of $x_0$ for which the geodesic segment appearing in the hypothesis is contained in $U$. The Hessian assumption remains valid on $U$ because $U \subset U_0$.
[/step]
[step:Compose the signed distance with the minimizing geodesic]
Let $p,q \in C \cap U$. If $p=q$, the unique minimizing geodesic is the constant geodesic with image $\{p\} \subset C \cap U$, so assume $p \ne q$.
Let $L := d_g(p,q) > 0$, and let
\begin{align*}
\gamma: [0,L] &\to U \\
t &\mapsto \gamma(t)
\end{align*}
be the unique minimizing geodesic from $p$ to $q$, parametrized by arclength. Thus $\gamma(0)=p$, $\gamma(L)=q$, $|\dot{\gamma}(t)|_g=1$, and $\nabla_{\dot{\gamma}}\dot{\gamma}=0$ for every $t \in [0,L]$.
Define
\begin{align*}
f: [0,L] &\to \mathbb{R} \\
t &\mapsto \rho(\gamma(t)).
\end{align*}
Since $\rho \in C^2(U)$ and $\gamma$ is smooth, $f$ is $C^2$ on $[0,L]$.
[/step]
[step:Use the geodesic equation to prove convexity of $\rho\circ\gamma$]
For each $t \in [0,L]$, the chain rule and the definition of the Riemannian Hessian give
\begin{align*}
f''(t)
&= \frac{d}{dt}\bigl(d\rho_{\gamma(t)}(\dot{\gamma}(t))\bigr) \\
&= \operatorname{Hess}\rho_{\gamma(t)}(\dot{\gamma}(t),\dot{\gamma}(t))
+ d\rho_{\gamma(t)}(\nabla_{\dot{\gamma}}\dot{\gamma}(t)).
\end{align*}
Because $\gamma$ is a geodesic, $\nabla_{\dot{\gamma}}\dot{\gamma}=0$. Therefore
\begin{align*}
f''(t)
= \operatorname{Hess}\rho_{\gamma(t)}(\dot{\gamma}(t),\dot{\gamma}(t))
\ge 0
\end{align*}
for every $t \in [0,L]$. Hence $f$ is convex on $[0,L]$.
[guided]
The point of introducing $f=\rho\circ\gamma$ is that membership in $C$ is detected by the sign of $\rho$. Thus, instead of studying the whole geodesic geometrically, we study the scalar function measuring its signed distance from the boundary.
We compute the second derivative of
\begin{align*}
f: [0,L] &\to \mathbb{R} \\
t &\mapsto \rho(\gamma(t)).
\end{align*}
The first derivative is
\begin{align*}
f'(t)=d\rho_{\gamma(t)}(\dot{\gamma}(t)).
\end{align*}
Differentiating once more along the curve gives
\begin{align*}
f''(t)
&= \frac{d}{dt}\bigl(d\rho_{\gamma(t)}(\dot{\gamma}(t))\bigr) \\
&= \operatorname{Hess}\rho_{\gamma(t)}(\dot{\gamma}(t),\dot{\gamma}(t))
+ d\rho_{\gamma(t)}(\nabla_{\dot{\gamma}}\dot{\gamma}(t)).
\end{align*}
This identity is the standard decomposition of the second derivative of a composition: the Hessian term measures the second variation of $\rho$ in the direction $\dot{\gamma}(t)$, while the final term measures the acceleration of the curve.
Here the acceleration term vanishes because $\gamma$ is a geodesic:
\begin{align*}
\nabla_{\dot{\gamma}}\dot{\gamma}=0.
\end{align*}
Therefore
\begin{align*}
f''(t)
= \operatorname{Hess}\rho_{\gamma(t)}(\dot{\gamma}(t),\dot{\gamma}(t)).
\end{align*}
The Hessian hypothesis applies at the point $y=\gamma(t)\in U$ and to the tangent vector $v=\dot{\gamma}(t)\in T_{\gamma(t)}M$, so
\begin{align*}
f''(t)\ge 0
\end{align*}
for every $t\in[0,L]$. A twice differentiable real-valued function on an interval with nonnegative second derivative is convex, so $f$ is convex on $[0,L]$.
[/guided]
[/step]
[step:Use endpoint nonpositivity to keep the geodesic inside $C$]
Since $p,q \in C \cap U$ and $C \cap U=\{x\in U:\rho(x)\le 0\}$, we have
\begin{align*}
f(0)=\rho(p)\le 0,
\qquad
f(L)=\rho(q)\le 0.
\end{align*}
For any $t \in [0,L]$, set $\lambda := t/L \in [0,1]$. Convexity of $f$ gives
\begin{align*}
f(t)
= f((1-\lambda)0+\lambda L)
\le (1-\lambda)f(0)+\lambda f(L)
\le 0.
\end{align*}
Thus $\rho(\gamma(t))\le 0$ for every $t\in[0,L]$. Since $\gamma(t)\in U$ and $C\cap U=\{x\in U:\rho(x)\le0\}$, it follows that $\gamma(t)\in C\cap U$ for every $t\in[0,L]$. Hence the image of the minimizing geodesic is contained in $C\cap U$.
[/step]
