How close are the exponential maps on $\mathbb{S}^2$ at two nearby points? Consider the two dimensional sphere $\mathbb{S}^2$ and let $p, q \in \mathbb{S}^2$. Let $\text{exp}_{p}$ and $\text{exp}_{q}$ be the exponential maps on $\mathbb{S}^2$ at points $p$ and $q$ respectively. I am interested in the map $\psi := \text{exp}_{p}^{-1} \circ \text {exp}_{q}$ defined on the unit disc $\mathbb{D} \subset \mathbb{R}^2$. I expect that if $p$ and $q$ are nearby points, then the map $\psi$ is close to the identity map. My question is: Is there a way to quantify this closeness? More precisely, is it possible to write the Taylor series expansion of $\psi = (\psi_1, \psi_2)$ around the origin as 
 \begin{align*}
 \psi_1(x,y) &= \psi_1(0) + \frac{\partial \psi_1}{\partial x}(0)~x + \frac{\partial \psi_1}{\partial y}(0)~ y + \ldots \\
 \psi_2(x,y) &= \psi_2(0) + \frac{\partial \psi_2}{\partial x}(0)~x + \frac{\partial \psi_2}{\partial y}(0)~ y + \ldots
 \end{align*}
 and show that the only significant coefficients in the above expansion are $\frac{\partial \psi_1}{\partial x}(0) \approx 1$, $\frac{\partial \psi_2}{\partial y}(0) \approx 1$ and all other partial derivatives $\frac{\partial^{m+n} \psi_i}{\partial^{m} x \partial^{n} y}(0) \approx 0$, for all other choices of $m,n$ and any $i = 1,2$.
Edit: After reading the comments/ answers of user35593 and Gro-Tsen , I realized that the following comments are essential regarding the above question:


*

*First is regrding the comment of user35593 (please read his/her comment below): The question makes sense only when we specify a certain way of identifying the tangent spaces at the points $p$ and $q$. Consider the rotation $R :\mathbb{R}^3 \rightarrow \mathbb{R}^3$ which maps $p$ to $q$, then it also maps $T_{p}(\mathbb{S}^2)$ to $T_{q}(\mathbb{S}^2)$ and it gives a way of identifying the two tangent spaces. The above question I asked is for the tangent spaces identified in this manner.

*Second is regarding the answer of Gro-Tsen : One possible interpretation of  ``$\frac{\partial^{m+n} \psi_i}{\partial^{m} x \partial^{n} y}(0) \approx 0$''  which I am interested in is: Given $\epsilon >0$, is it possible to find $\delta >0$ such that whenever $d_{\mathbb{S}^2}(p,q)\leq \delta$ we have
\begin{align*}
  \left| \frac{\partial^{m+n} \psi_i}{\partial^{m} x \partial^{n} y}(0)  \right| \leq \epsilon
  \end{align*}
whenever $m+n \leq 100$ (except for the cases of $\frac{\partial \psi_1}{\partial x} (0)$ and $\frac{\partial \psi_2}{\partial y} (0)$ which I expect to be $\approx 1$).
Thanks!
 A: I believe the following proves that the partial derivatives of $\psi$ at the origin cannot be bounded by a constant (so they certainly cannot be $\approx 0$ for any reasonable meaning of this symbol).
Assume on the contrary that (for some fixed $p,q$ distinct and not antipodal), the partial derivatives of $\psi := \exp_p^{-1} \circ \exp_q$ (defined in some neighborhood of the origin) are all bounded by a constant.  Then, by summing the Taylor series expansion of $\psi$ at $0$ we see that $\psi$ extends to a real-analytic function $\psi\colon\mathbb{R}^2 \to \mathbb{R}^2$, which by analytic extension must still satisfy $\exp_p \circ \psi = \exp_q$.  Let me argue why this is impossible.
We can assume w.l.o.g. that the coordinates $(u,v)$ on the tangent plane to $\mathbb{S}^2$ at $q$ were chosen so that $\exp_q$ maps the axis $(u,0)$ to the great circle connecting $q$ and $p$, and more precisely, if $0<\delta<\pi$ is the distance between $q$ and $p$ on $\mathbb{S}^2$, that $\exp_q$ takes $(\delta,0)$ to $p$.  Furthermore, we can similarly assume on the coordinates $(u',v')$ of the tangent plane at $p$ that $\exp_p$ maps the axis $(u',0)$ to the same great circle and takes $(-\delta,0)$ to $q$.  Then $\exp_q(u,0)$ is the point obtained by traveling a distance $u$ on $\mathbb{S}^2$ starting from $q$ in the direction of $p$, and $\exp_p(u,0)$ is the point obtained by traveling a distance $u$ on $\mathbb{S}^2$ starting from $p$ in the direction opposite to $p$, thus $\psi(u,0) = (u-\delta,0)$ for $u$ in the neighborhood of $0$, hence everywhere by analytic extension.
On the other hand, if $(u,v)$ lies on the circle $C$ with radius $\pi$ around the origin then $\exp_q$ takes $(u,v)$ to the antipode $\tilde q$ of $q$.  But the inverse image of $\tilde q$ by $\exp_p$ is discrete (since $\exp_p$ is a diffeomorphism outside of circles of radius $k\pi$ around the origin which are mapped to either $p$ or its antipode $\tilde p$, and we are assuming $p,q,\tilde p,\tilde q$ distinct); and $\psi$ must map $C$ (which is connected) inside this inverse image: so $\psi$ must be constant on $C$.  But this contradicts the fact that $\psi(\pi,0) = (\pi-\delta,0)$ and $\psi(-\pi,0) = (-\pi-\delta,0)$ (as per previous paragraph) are not equal.
A: Not a solution but too long for a comment.
Note that the map depends only on the distance $d(p,q)$ of $p,q$, i.e. we have $\phi(t)\colon \mathbb{R}^2\rightarrow \mathbb{R}^2$ where $t=d(p,q)$. Note further that $\phi(t_1+t_2)=\phi(t_1)\circ \phi(t_2)$, i.e. $\phi(t)$ is a one-parameter group. Let $D=d\phi(t)/dt|_{t=0}$. Note that then $d\phi(t)/dt=D(\phi(t))$. Hence $\phi(T)(u)$ can be found by solving the ODE $f'(t)=D(f(t))$ on the interval $[0,T]$ with initial condition $f(0)=u$. Let us now find $D$ for our case of the sphere.
The exponential and its inverse on the sphere are $exp(v)=\begin{pmatrix} sin(|v|)\frac{v}{|v|}\\ cos(|v|)
\end{pmatrix}$ and $log \left(\begin{pmatrix} x\\ y\\z \end{pmatrix}\right)=\frac{arccos(z)}{\sqrt{1-z^2}}\begin{pmatrix} x\\ y \end{pmatrix}$. Applying and infinitesimal rotation by $\alpha$ after the exponential map yields
$$
\begin{pmatrix} sin(|v|)\frac{v_1}{|v|}-\alpha \cdot cos(|v|)\\sin(|v|)\frac{v_2}{|v|} \\ cos(|v|)+\alpha \cdot sin(|v|)\frac{v_1}{|v|}
\end{pmatrix}
$$
Now applying the logarithm using
$$
d\frac{arccos(z)}{\sqrt{1-z^2}}/dz=\frac{arccos(z)z-\sqrt{1-z^2}}{\sqrt{1-z^2}^3}
$$
yields
$$
v+\alpha \cdot sin(|v|)\frac{v_1}{|v|} \frac{(|v|cos(|v|)-sin(|v|))}{sin(|v|)^3} sin(|v|)\frac{v}{|v|} -\begin{pmatrix} \alpha \frac{cos(|v|)}{sin(|v|)}|v|\\0\end{pmatrix}\\
=v+\alpha \cdot \frac{(|v|cos(|v|)-sin(|v|))v_1}{sin(|v|)|v|^2} v -\begin{pmatrix} \alpha \frac{cos(|v|)}{sin(|v|)}|v|\\0\end{pmatrix}
$$
Hence
$$
D(v)=\frac{(|v|cos(|v|)-sin(|v|))v_1}{sin(|v|)|v|^2} v-\begin{pmatrix} \frac{cos(|v|)}{sin(|v|)}|v|\\0\end{pmatrix}.
$$
