Change of two normal coordinates based on two nearby points? Let $M$ be a manifold and $L(M)$ be the tangent frame bundle on $M$. Let $\Gamma$ be a linear connection on $L(M)$ which induces a covariant derivative $\nabla$ on $TM$.
Let $p, q$ be two distinguished and sufficiently close points on $M$, connected by a smooth curve $\gamma = \{\gamma_t\}_{t\in[0,\epsilon]}$ so that $\gamma(0)=p$, $\gamma(\epsilon)=q$. Suppose that there is an open set $U\subset M$, such that $U$ contains $\gamma$ and is a normal neighborhood of both $p$ and $q$.
Fix a linear frame $u_p \in L_p(M): \mathbf R^d\to T_pM$. Let $u_q = \Gamma(\gamma)_0^\epsilon (u_p) \in L_q(M)$ be the parallel displacement of $u_p$ along $\gamma$, that is, $u_q$ can be joined to $u_p$ by a horizontal curve on $L(M)$ along $\gamma$.
Then we have two normal coordinate systems on $U$:
$$x=u_p^{-1}\circ \exp_p^{-1}: U \to \mathbf R^d,$$
$$y=u_q^{-1}\circ \exp_q^{-1}: U \to \mathbf R^d,$$
so that $(U,(x^i))$ and $(U,(y^j))$ are coordinate charts based on $p$ and $q$ respectively.

Now my question is: how to do change of coordinates between these two normal coordinate charts? More precisely, let $m \in U$, then what is the relation between $x(m)$ and $y(m)$?
If necessary, you can endow more structures to $M$. Say, $M$ is equipped with a Riemannian metric $g$ and $\nabla$ is the Levi-Civita connection, the frame bundle $L(M)$ is replaced by the orthonormal frame bundle $O(M)$...

Some thinking: Clearly, $x(p) = y(q) = 0$. If $\gamma$ is a geodesic, then it is easy to check $x(q) = -y(p)$. Denote $\rho = x(q) \in\mathbf R^d$. Then I believe that $x(m) = y(m) + \rho$ for all $m \in U$. If $\gamma$ is not a geodesic, then I think that a curvature term should appear, since we may use the holonomy. But I do not know how to prove my conjectures...

EDIT: I think there maybe no exact expressions for $x(m)$ and $y(m)$, but there should be an asymptotic expression with infinitesimals $o(\epsilon)$ or $o(d(x,y))$...
 A: Thanks for the instructive comments by @Deane, and by @WillieWong in the crosspost on mathoverflow.

First of all, the shift property $x(m) = y(m) + \rho$ that I conjectured in my question is not true in general. See the comments by @WillieWong for a counterexample.

Secondly, inspired by their comments, I am now trying to give a feasible option instead of a complete solution, as the question is actually quite general.
Denote
\begin{align*}
u_t &= \Gamma(\gamma)_0^t (u_p) \in L_{\gamma(t)}(M), \\
x_t &= u_t^{-1} \circ \exp_{\gamma(t)}^{-1} (m) \in \mathbf R^d, \\
V(t) &= u_t x_t = \exp_{\gamma(t)}^{-1} (m) \in T_{\gamma(t)} M, \\
\gamma_s(t) &=  \exp_{\gamma(t)}(sV(t)) \in M.
\end{align*}
Then for each $t$, $\gamma_\cdot(t)$ is a geodesic satisfying $\gamma_0(t) = \gamma(t)$ and $\gamma_1(t) \equiv m$. This leads to a Jacobi field $J_t$ along $\gamma_\cdot(t)$ for each $t$:
\begin{equation*}
J_t(s) = \frac{\partial}{\partial t} \gamma_s(t) = \dot\gamma_s(t),
\end{equation*}
which satisfies the Jacobi equation:
\begin{equation}\tag{1}
\frac{D^2}{ds^2} J_t(s) + R\left( J_t(s), \frac{\partial}{\partial s} \gamma_s(t) \right) \frac{\partial}{\partial s} \gamma_s(t) = 0,
\end{equation}
subjected to two boundary conditions:
\begin{equation}\tag{2}
J_t(0) = \dot\gamma_0(t) = \dot\gamma(t),
\end{equation}
and
\begin{equation}\tag{3}
J_t(1) = \dot\gamma_1(t) = 0.
\end{equation}
Now the torsion-freeness (which needs to be assumed in the beginning) gives
\begin{equation*}
\begin{split}
\frac{D}{ds} J_t(0) &= \frac{D}{ds} \frac{\partial}{\partial t} \gamma_s(t) \bigg|_{s=0} = \frac{D}{dt} \frac{\partial}{\partial s} \gamma_s(t) \bigg|_{s=0} = \frac{D}{dt} V(t) \\
&= \frac{d}{dh}\bigg|_{h=0} u_t\circ u_{t+h}^{-1} (V(t+h)) = \frac{d}{dh}\bigg|_{h=0} u_t x_{t+h} = u_t \dot x_t,
\end{split}
\end{equation*}
which gives
\begin{equation}\tag{4}
\dot x_t = u_t^{-1} \left( \frac{D}{ds} J_t(0) \right).
\end{equation}
Therefore, in order to obtain $\dot x_t$ which is exactly the purpose of the question, one need first solve the boundary value problem $(1)$-$(3)$ to get $J_t(s)$, and then calculate $\frac{D}{ds} J_t(0)$ and finally use $(4)$ to derive $\dot x_t$.
