The last question is too general, this is a modification.

Let $M$ be an $n$ dimensional Riemannian manifold. Fix $m$ points $x_1,...,x_m$. Suppose $y$ is not in the cut locus of $x_i$ for $1 \leqslant i \leqslant m$. Then the function $g(x)=\sum_{i=1}^m d^2(x_i,x)$ is $C^2$ near $y$. Thus there is a small neighborhood $B_r(y)$ of $y$ and a constant $C$ such that $|\nabla^2 g(x)|\leqslant C$ for any $x \in B_r(y)$.

My question: Is the map $x \mapsto \exp_ x(\nabla g(x))$ Lipschitz in $B_r(y)$?

My understanding (I add an condition $sec(M) \geqslant 0$ to simplify my explanation, and I am only interested in manifolds with lower curvature bound): To prove $\exp_x(\nabla g(x))$ is Lipschitz, we just need to show that for any point $x\in B_r(y)$, any unit geodesic $\gamma(t)$ with $\gamma(0)=x$, we have $$ \frac{d}{dt}|_{t=0} dist(\exp_{\gamma(t)}(\nabla g(\gamma(t))), exp_x(\nabla g(x)) \leqslant C_1, $$ and $C_1$ is the same for all the points in the ball. If this is done, then $\exp_z(\nabla g(x))$ is $C_1$-Lipschitz.

Since $|\nabla^2 g| \leqslant C$, we have $|\nabla_{\dot{\gamma}(0)}\nabla g|\leqslant C$, that is, $$ \frac{d}{dt}[{P_t}^{-1}(\nabla g(\gamma(t))]\leqslant C $$ Where $P_t: T_xM \mapsto T_{\gamma(t)}M$ is the parallel transportation. Then \begin{equation*} \begin{array}{ll} dist(\exp_x (\nabla g(x)),\exp_{\gamma(t)} (\nabla g(\gamma(t)))&\leqslant dist\left(\exp_x (\nabla g(x)),\exp_x ({P_t}^{-1}(\nabla g(\gamma(t))\right)\\ &+dist\left(\exp_x ({P_t}^{-1}(\nabla g(\gamma(t)), \exp_{\gamma(t)} (\nabla g(\gamma(t)\right) \end{array} \end{equation*} If the curvature $secM \geqslant 0$, then the exponential map is a contracting map (i.e. 1-Lipschitz), the first term on the RHS is less than $$ |\nabla g(x), {P_t}^{-1}(\nabla g(\gamma(t))|\leqslant Ct+o(t). $$

But for the second term on the RHS, I don't know how to estimate. This is what I want to ask.