Formulation of question: Consider a complete, simply connected Riemannian manifold $M$ with Riemannian metric $d$. For $x,y \in M$ that are distinct but close enough to each other, there is a geodesic $\gamma(t)$, $0 \le t \le 1$ such that $\gamma(0)=x$ and $\gamma(1)=y$. Assume that $\gamma$ is the unique minimizing geodesic joining $x$ and $y$ and that $y$ is the first conjugate point of $x$ along $\gamma$. Further, assume that $y$ has a suitable normal neighborhood that is contained in $M$. Is $d(x,y)$ differentiable at $y$ when $x$ is fixed?
Background: We know that cut loci carry important information on the topology of $M$ and information on how the distance function $d(x,\cdot)$ behaves at a cut point also reveals important information on $M$. Further, a lot of techniques in Riemannian geometry utilize the differentiability of $d(x,\cdot)$, and in statistics, uniqueness of Fréchet mean, i.e., $$\arg \min_{y \in M} \int_{M} d^2(x,y) \mu(dx)$$ for a probability measure $\mu$ on $M$, depends critically on the cut loci of $d(x,\cdot)$.
Related posts and results: there is a nice post smoothness of $d(\cdot,\cdot)$ that discusses differentiability of $d$ on $M \times M$, and a classic result on page 108, i.e., Proposition 4.8, in Takashi Sakai's classic book "Riemannian Geometry", that says "Suppose that there exist at least two normal minimal geodesics joining" $x$ and $y$, then $d(x,\cdot)$ is not differentiable at $y$. These discussions confirm the differentiability of $d(x,y)$ when $y$ is not a cut point to $x$ when $x$ is fixed.
Recall that a cut point $z$ of $x=\gamma(0)$ is either a first conjugate point to $x$ or there are at least two minimal geodesics joining $x$ and $z$. Recall also that from Theorem 2.1.12 on page 133 of W.A. Klingenberg's book "Riemannian Geometry", within any neighborhood of the first conjugate point $y$, there exists a point $w$ such that there are at least two distinct geodesics that connects $x$ and $w$. So, if $d(x,y)$ with fixed $x$ is indeed differentiable at $y$, then we also know that $d(x,\cdot)$ is not differentiable in any neighborhood of $y$ except at $y$ itself.
However, none of the above discusses whether $d(x,y)$ is differentiable at $y$ in the setting of the question stated at the beginning. Is is trivial that the answer is "No"? If the question is so easy, I must have missed something and would like to be enlightened by someone who knows the subject better. Thank you for your attention to my question.
Update and answer (Aug 22, 2023): I spent some time trying to write a proof based on what I was told in the update on Feb 21, 2023 (i.e., "suppose $d(x,y)$, $x \ne y$ and $x$ fixed, is differentiable at $y$, then $\vert \nabla d(x,y)\vert = 1$ has to hold. But when $y$ is a conjugate point, we have $\nabla d(x,y) = 0$, a contradiction.") but was not successful since I did not arrive at the contradiction. This question has been bugging me continually for at least 1.5 years.
Finally, I found Richard L. Bishop's 1977 paper "Decomposition of cut loci", where he stated in the middle of page 134 of the paper the following:
The smoothness properties of $r$ at a singular cut-point can be more complicated. We can make calculations, as indicated later in this paragraph, on surfaces, which show that $r$ can be $C^1$ but not $C^2$.
Here $r$ is the distance function from a fixed point $m$, a "singular cut-point" is a cut point to a fixed point that is a first conjugate point with a unique minimal geodesic connecting this point and the fixed point, and $C^1$ means "continuously differentiable" and $C^2$ "twice continuously differentiable".
In the paper, Richard L. Bishop stated a conjecture on the generality of the above phenomenon.
In addition, the paper "Cut loci and distance functions" by Jin-Ichi Itoh and Takashi Sakai states the following:
On the other hand, if $q \in C(p)$ then $d_p$ is in general not differentiable at $q$
Here $d_p$ is the distance from $p$, and $C(p)$ is the cut locus of $p$. The above statement implies that $d_p$ can be differentiable under the settings of the question.
So, the general conclusion is the following:
On a Riemannian manifold $(M,d)$, the distance function from a fixed point $p$, i.e., $f_p(x)=d(p,x), p, x \in M$ can still be differentiable at a first conjugate point $q \ne p$ of $p$ when there is a unique minimizing geodesic connecting $p$ and $q$.
Due to this above findings, I have deleted the "answer".
Comment (Aug 22, 2023) Right under Lemma 5.7.9 on page 218 of the 3rd Edition of Peter Petersen's book "Riemannian Geometry", the author says the following "Notice that in the first case (where there are at least 2 minimal geodesics between the fixed point and its cut point) the gradient (of the distance function) is undefined (at this cut point) …; while in the second case (where the exponential map is singular, i.e., at a first conjugate point) the Hessian (of the distance function) becomes undefined (at this conjugate point) …".
This statement has been given full details by the answer to the post Are distance functions necessarily nonsmooth on the cut locus?.
So, in summary, here is the full picture:
The distance function $f_p(x)=d(p,x)$, $p, x \in M$ cannot be smooth at any point $q \in C(p)$, where $C(p)$ is the cut locus of $p$. It is not differentiable at $q \in C(p)$ when there are at least 2 minimal geodesics connecting $p$ and $q$. But it can be differentiable if $q \ne p$ is a first conjugate point with a unique minimal geodesic connecting $p$ and $q$; in this case, it cannot be differentiable up to the 2nd order though.
Here "a smooth function" refers to "a function whose component functions have continuous partial derivatives of all orders", following the definition of John M. Lee's book "Introduction to smooth manifolds". Obviously, smoothness is a much stronger than differentiability, and for the Fréchet function $F(y)= \int_{M} d^2(x,y) \mu(dx)$, usually 2nd order continuous differentiability, i.e., $C^2$ suffices.
Update (Aug 30,2023) Proposition 1.8 of the paper "Optimal transport and curvature" by Alessio Figalli and Cedric Villani states, in the notations here, "$d^2(p,x)$, $p, x \in M$ for a fixed $p$ is differentiable at $x$ if $p$ and $x$ are joined by a unique minimizing geodesic", and right after Proposition 1.9 there, the authors state, in the notations here, "$d^2(p,x)$ is differentiable at $x$ if and only if $p$ and $x$ are joined by a unique minimizing geodesic." Namely, the squared distance function from a base point to another is differentiable at the other point if and only if these two points are connected by a unique minimizing geodesic.
The proof of these two results is based on the first variation formula of the energy functional, and it only depends on the $C^2$ smoothness (i.e., $C^2$ atlas) and completeness (i.e., geodesics extending infinitely) of a Riemannian manifold.
Hopefully this dissolves the wide spread misconception among non-experts in differential geometry that are self-learning that the distance function is not differential inside a cut locus.