You can probably make your idea work but it won't be easy. The reason is that your formula that recovers the connection from the parallel transport is true not only for the Levi-Civita connection but also for arbitrary connections. This means that in order to identify the right hand side as the Levi-Civita connection, you will need to understand what makes the parallel transport of the Levi-Civita connection special compared with the parallel transport of a general connection. The compatibility with the metric is easy - this implies that parallel transport is an isometry. However, to understand how the symmetry affects the parallel transport is much more delicate (see here for example).
A much less painful way to solve the exercise is to use the notion of a parallel frame along $c$ (which if I remember correctly is introduced in one of the other exercises). Namely, pick some basis of $\xi_1(p), \dots, \xi_n(p)$ of $T_pM$ and extend it by parallel transport to a frame $(\xi_1, \dots, \xi_n)$ of vectors fields along $c$. Now, write the restriction of $Y$ to $c(t)$ as $Y = Y^i(t) \xi_i(c(t))$ (summation convention is in use) and note that
$$ (\nabla_X Y)(c(t)) = \frac{DY(c(t))}{dt} = \dot{Y}^i(t) \xi_i(c(t)) + Y^i(t) \frac{D\xi_i(t)}{dt} = \dot{Y}^i(t) \xi_i (c(t)) $$
which means that the covariant derivative relative to the frame $\xi_i$ is given simply by the regular derivative. Then,
$$ \frac{d}{dt} \left( P_{c,t_0, t}^{-1}(Y(c(t)) \right)|_{t = t_0} = \frac{d}{dt} \left( P_{c,t_0, t}^{-1}(Y^i(t) \xi_i(c(t))) \right)|_{t = t_0} \\
= \frac{d}{dt} \left( Y^i(t) \xi_i(p) \right)|_{t = t_0} = \dot{Y}^i(t_0) \xi_i(p) = (\nabla_X Y)(p).$$
Note: The Einstein summation convention is in force throughout this answer.End of Note.
Our OP walaa's question,
Now are the geodesics always smooth in a Riemannian metric?
appears to be somewhat of a misstatement, insofar as smoothness is technically defined with respect to local coordinate patches of a given atlas; both the smoothness of the metric tensor $g_{ij}$ and of the Christoffel symbols $\Gamma_{ij}^k$ are then affirmed if these functions are smooth in such coordinates. A more carefully stated question along these lines might read:
Are the geodesics of a smooth Riemannian metric smooth?
This is indeed the case, as may be seen by observing that the coefficients $\Gamma_{ij}^k$ occurring in the geodesic equation
$\ddot \gamma = \dfrac{d^2 \gamma^k}{dt^2} \partial_k + \Gamma_{ij}^k \dfrac{d\gamma^i}{dt} \dfrac{d\gamma^j}{dt} \partial_k = 0 \tag 1$
are themselves smooth functions on $M$, being given in terms of the $g_{ij}$ by
$\Gamma_{ij}^k = \dfrac{1}{2}g^{km}(g_{mi, j} + g_{mj, i} - g_{ij, m}), \tag 2$
where
$[g^{ij}] = [g_{ij}]^{-1}; \tag{2.5}$
in these equations we understand that
$\partial_k = \dfrac{\partial}{\partial x_k}, \tag 3$
and
$g_{ij, m} = \dfrac{\partial g_{ij}}{\partial x_m}, \tag 4$
and so forth.
We may cast (1) in the familiar first order form by setting
$\dfrac{d\gamma^i}{dt} = \beta^i, \; 1 \le i \le \dim M; \tag 5$
then
$\dfrac{d^2\gamma^i}{dt^2} = \dfrac{d\beta^i}{dt}; \tag 6$
(1) may then be written
$\dfrac{d\beta^k}{dt}\partial_k + \Gamma_{ij}^k \beta^i \beta^j \partial_k = 0, \tag 7$
which corresponds to the collection of $\dim M$ first order, non-linear ordinary differential equations
$\dfrac{d\beta^k}{dt} + \Gamma_{ij}^k \beta^i \beta^j \ = 0, \; 1 \le k \le \dim M, \tag 8$
that is,
$\dfrac{d\beta^k}{dt} = -\Gamma_{ij}^k \beta^i \beta^j \ = 0, \; 1 \le k \le \dim M; \tag 9$
since the $\Gamma_{ij}^k$ are functions of the position coordinates $x = (x_1, x_2, \ldots, x_{\dim M})$ in $M$, along $\gamma(t)$ we have
$\Gamma_{ij}^k(x) = \Gamma_{ij}^k(\gamma(t)) = \Gamma_{ij}^k(\gamma^1(t), \gamma^2(t), \ldots, \gamma^{\dim M}(t));\tag{10}$
(9) then becomes
$\dfrac{d\beta^k}{dt} = - \Gamma_{ij}^k(\gamma^1(t), \gamma^2(t), \ldots, \gamma^{\dim M}(t))\beta^i \beta^j = 0, \; 1 \le k \le \dim M; \tag{11}$
these together with (5) form a system of $2\dim M$ ordinary, non-linear differential equations for the $\gamma^i(t)$, $\beta^i(t)$, $1 \le i \le \dim M$.
We observe that if the functions $\Gamma_{ij}^k(\gamma^1, \gamma^2, \ldots, \gamma^{\dim M})\beta^i \beta^j$ occurring on the right-hand side of (11) are in fact $C^m$ in the $\gamma^i$, $\beta^i$, then
(5), (11) forms a $C^m$ system of non-linear, ordinary differential equations for the $\gamma^i$, $\beta^i$. It then follows from the standard theorems that the solution functions $\gamma^i(t)$, $\beta^i(t)$, $1 \le i \le \dim M$ are themselves at least $C^m$; since the $\Gamma_{ij}^k$ in fact depend on the first derivatives of the $g_{ij}$ (cf. (2)), we may infer that for $C^m$-smooth $g_{ij}$ the $\Gamma_{ij}^k$ are $C^{m - 1}$-smooth, and hence so the $\gamma_i(t)$, $\beta_i(t)$ are at least $C^{m - 1}$. Since these assertions bind for any $m \ge 1$ it follows that if the metric $g_{ij}$ is $C^\infty$, so are the functions $\gamma_i(t)$, $\beta_i(t)$, $1 \le i \le \dim M$. Thus the geodesics of a smooth Riemannian metric are themselves smooth.
Best Answer
Let $M$ be a manifold such that there is a length function, where each smooth curve $c: [a,b] \rightarrow M$ has a length $\ell[c]$ and $\ell$ satisfies some expected geometric properties. Then you can define a distance function between two points as the infimum of lengths of curves connecting the two points. A geodesic is now a curve such that given any two sufficiently close points, the length of the curve segment connecting the two points is equal to the distance between the two points. Lots of things are hard to do in this setting. Even the existence of any geodesics is in question.
So far, no connection is needed. It is possible to proceed further without a connection. However, the above can be done within the context of Riemannian geometry, and it turns out that a geodesic, with respect to a constant speed parameterization satisfies the geodesic equation $$ \nabla_{c'}c' = 0. $$ Equivalently, the velocity vector $c'$ is parallel along the geodesic $c$.
This equation, however, does not use the Riemannian metric itself, only the Levi-Civita connection. It follows that if you have a connection on a manifold, then you can define geodesics using the geodesic equation, even if the connection is not the Levi-Civita connection of a Riemannian metric. In particular, no distance or length function is used for this definition of a geodesic. This approach is important in general relativity for defining null geodesics.