With metric tensors of the unit sphere in normal coordinates, their Taylor series for $p\in S$ near the north pole $N$ can be written as follows.
$$g_{rr}(p) \equiv 1; g_{r\theta}(p) = g_{\theta r}(p) \equiv 0; g_{\theta\theta}(p) = r^2.$$
However, the third expansion does not check out with the following general expression
$$g_{ij} = \delta_{ij} - \frac{1}{3} \sum_{k, l} R_{ikjl}(0) u^k u^l + O(|u|^3). $$
$$ g_{22} = g_{\theta\theta} = 1 + \frac{1}{3} (\sin^2r)r^2 + O(|u|^3) $$
Where did I make mistake, please? Thank you!
The formula
$$ g_{ij} = \delta_{ij} - \frac{1}{3} \sum_{k, l} R_{ikjl}(0) u^k u^l + O(|u|^3) $$
holds only for geodesic normal coordinates - that is, for the standard cartesian coordinate system $(x,y)$ on $T_p M \simeq \mathbb R^n$ transferred to $M$ via the exponential map. It does not hold for the geodesic polar coordinates $(r,\theta)$ - in fact these coordinates are not even well-defined at $p$, where $\theta$ is undefined with the corresponding degeneracy $g_{\theta\theta} = 0$ in the coordinate expression for the metric.
If you transform to geodesic normal coordinates using $x = r \cos \theta, y = r \sin \theta$ then the formula should hold for $g_{xx}$ etc.