I was told that if we consider a circle of a small radius $\varepsilon$ on a two-dimensional Riamannian manifold then it's length and area are $$ L(\varepsilon) = 2 \pi \varepsilon + A \varepsilon^3 + o(\varepsilon^3),\qquad S(\varepsilon) = \pi \varepsilon^2 + B \varepsilon^4 + o(\varepsilon^4), $$ where $A$ and $B$ are some functions of Gaussian curvature $K$. I can directly verify it for a sphere to get $A = - \frac{\pi}{3} K$, $B = - \frac{\pi}{12} K$. I proved the identities above (without explicitly finding $A$ and $B$) and want my proof to be checked.
Consider orthonormal basis $e_a$ at tangent space at point $p$ and consider Riemann normal coordinates $x^a$ around it. In this coordinates geodesic starting at $p$ along vector $v = v^a e_a$ has coordinates $$ x^a(t) = v^a t $$ and Christoffel symbols vanish at $p$ so $g_{ab,c}(p) = 0$, so the metric is $$ g_{ab}(x) = \delta_{ab} + \frac12 g_{ab,cd}(p) x^c x^d + o(x^2). \qquad (*) $$ Consider geodesic that starts at $p$, has initial velocity $v$, $|v| = 1$, which has length $\varepsilon$. Then its end has coordinates $x^i_v = v^i \varepsilon + o(\varepsilon)$.
Consider unit vectors at $p$ $v(\varphi) = \cos \varphi e_1 + \sin \varphi e_2$. Then the length of a circle centered at $p$ with radius $\varepsilon$ is $$ L(\varepsilon) = \int_0^{2\pi} \left|\frac{dx^a}{d\varphi}\right| d \varphi = \int_0^{2\pi} \sqrt{\varepsilon^2 \delta_{ij}v^i{}'v^j{}' + \varepsilon^4 \frac12 g_{ij,kl}(p)v^i{}'v^j{}'v^k{}'v^l{}' + o(\varepsilon^4)} d\varphi\\ = \int_0^{2\pi} \left( \varepsilon + \varepsilon^3 \frac14 g_{ij,kl}(p)v^i{}'v^j{}'v^k{}'v^l{}' + o(\varepsilon^3) \right) d\varphi\\ = 2 \pi \varepsilon + \varepsilon^3 \frac14 g_{ij,kl}(p) \int_0^{2\pi} v^i{}'v^j{}'v^k{}'v^l{}' d\varphi + o(\varepsilon^3), $$ where $v' = \frac{dv}{d\varphi} = -\sin \varphi e_1 + \cos \varphi e_2$. I didn't check whether $$ g_{ij,kl}(p) \int_0^{2\pi} v^i{}'v^j{}'v^k{}'v^l{}' d\varphi = \frac{\pi}{4} \left( 3 (g_{11,11} + g_{22,22}) + 4 g_{12,12} \right). $$ is proportional to scalar curvature, it's too cumbersome and I don't see any way to check it in not computationally awful way.
Consider ''polar coordinates'' $(r, \varphi)$ s.t. $r = \sqrt{(x^1)^2 + (x^2)^2}$, $\tan \varphi = \frac{x^2}{x^1}$. Note that for a point on a circle of radius $\varepsilon$ we have $r = \varepsilon + o(\varepsilon)$.
The volume form is $\omega = \sqrt{g} dx^1 \wedge dx^2$. Using coordinate expansion for metric $(*)$ we have $$ \omega = \sqrt{g} dx^1 \wedge dx^2 = (1 + r^2 C + o(r^2)) dx^1 \wedge dx^2 = (1 + r^2 C + o(r^2)) r dr \wedge d\varphi, $$ where $C$ is a function of $v$, that is a function of $\varphi$. So finally $$ S(\varepsilon) = \int_0^{2\pi} d\varphi \int_0^\varepsilon r dr + \int_0^{2\pi} C(\varphi) d\varphi \int_0^\varepsilon r^3 dr + o(\varepsilon^4) = \pi \varepsilon^2 + \varepsilon^4 \frac14 \int_0^{2\pi} C(\varphi) d\varphi + o(\varepsilon^4). $$