Let $(M, g)$ be a Riemannian manifold (without boundary) and let $p \in M$. Consider a linear $2$-plane $P$ in $T_p M$ and let $C_r$ be a circle in $P$, of radius $r$ and centered at $0$ (for $r$ small enough). How can we prove that $$\mathrm{Length}(\mathrm{exp}_p(C_r)) = 2 \pi \left(r - \frac{\mathrm{sec}(P)}{6} r^3 + O(r^4) \right), $$ where $\mathrm{exp}_p$ is the Riemannian exponential map at $p$ and $\mathrm{sec}(P)$ is the sectional curvature of the plane $P$?
I tried using the definition of length: the curve $\mathrm{exp}_p(C_r)$ can be parametrized as $$\gamma : [0,2\pi] \to M, \quad \gamma(t) = \mathrm{exp}_p (r\cos(t) v + r\sin(t) w), $$ where $\{v, w\}$ is an orthonormal basis for $P$. Then, we know that $$\mathrm{sec}(P) = R_p(v, w, w, v), $$ where $R$ is the $(0, 4)$ version of the Riemann curvature tensor, and we have that, \begin{align} \mathrm{Length} (\gamma) &= \int_0^{2\pi} |\dot{\gamma}(t)|dt \\ &= \int_0^{2\pi} \left| \left( d \left(\mathrm{exp}_p \right) \right)_{r\cos(t) v + r\sin(t) w} \left( -r\sin(t)v + r\cos(t)w \right) \right| dt, \end{align} but I am not sure how to continue further, as I do not exactly know how to compute the differential of the exponential map at those particular points.
I thought about using the fact that, since $r$ is small enough, the image of $\gamma$ is contained in a chart of Riemannian normal coordinates, call them $(x^1, \cdots, x^n)$. Then, in these coordinates, the components of the Riemannian metric satisfy $$g_{ij} = \delta_{ij} - \frac{1}{3} \sum_{k, l} R_{iklj}(p)x^k x^l + O(|x|^3), $$ where $R_{ikjl}$ are the components of the $(0,4)$ version of the Riemann curvature tensor (this is Exercise 10.1 in Lee's Introduction to Riemannian manifolds, 2nd edition). However, I do not know how to use this.
The solution you are looking for can be found there: Theorem 3.68, page 153, in Riemannian geometry, by Gallot, Hulin & Lafontaine, second edition.
(Community wiki, as it is a reference only answer.)