I want to show that, on a $2$-dimensional Riemannian manifold $(M,g)$, under the conformal change of metric $\tilde g= e^{2f} g$ for some $f\in C^\infty(M)$, we have the following identity about scalar curvatures,
$$s_{\tilde g} = e^{-2f}(2\Delta_g(f)+s_g). $$
I know that we have
$$s_g = \frac{2R_{1221}}{\det(g_{ij})},$$
in local coordinates $(x^1,x^2)$, and
$$\tilde \nabla_XY = \nabla_X Y + \left(Y(f)X +X(f)Y-g(X,Y)grad_g(f)\right), $$
for any $X,Y \in \mathfrak{X}(M)$. The hint that I was given is that to use geodesic normal coordinates, for which we have
$$ \frac{\partial g_{ij}}{\partial x^k}(p) = 0,\quad \Gamma_{ij}^k(p)= 0, $$
where $p$ is the point where the coordinates are centered at.
The rest should be a straight forward computation, but things really turn out to be excruciating for computing $\tilde{R}_{1221}$. It is easy to obtain that
$$\frac{2\tilde{R}_{1221}} {\det(\tilde{g}_{ij})} = e^{-2f} \left(\frac{2R_{1221}+ \text{something}}{\det(g_{ij})} \right),$$
so it remains only to compute that the "something" is equal to $2\Delta_g(f)\det(g_{ij}) = 2\det(g_{ij})\sum_{i,j} g^{ij}\frac{\partial^2 f}{\partial x^i\partial x^j}$. Using the two identities given above by geodesic normal coordinates, I've expanded the "something" into $g(A,\frac{\partial}{\partial x^1})$, where $A$ is given by
$$\begin{aligned} A = &f_{12} \frac{\partial}{\partial x^2} - g_{22} \left( \sum_{i,j} g^{ij} f_{1i}\frac{\partial}{\partial x^j}\right)+ f_1\left(2f_2 \frac{\partial}{\partial x^2} - g_{22} \sum_{i,j}g^{ij} f_i\frac{\partial}{\partial x^j}\right) +\left( (f_2)^2 - g_{22} \sum_{i,j}g^{ij} f_if_j\right)\frac{\partial}{\partial x^1} \\ &-g\left(\frac{\partial}{\partial x^1},2f_2 \frac{\partial}{\partial x^2} - g_{22} \sum_{i,j} g^{ij} f_i \frac{\partial}{\partial x^j}\right)\left(\sum_{i,j} g^{ij} f_i\frac{\partial}{\partial x^j}\right) - f_{22} \frac{\partial}{\partial x^1} +g_{12} \left(\sum_{i,j} g^{ij} f_{2i} \frac{\partial}{\partial x^j} \right)\\ &- f_2\left(f_1\frac{\partial}{\partial x^2} - g_{12} \sum_{i,j}g^{ij}f_i\frac{\partial}{\partial x^j}\right) - \left(2f_1f_2 - g_{12}\sum_{i,j}g^{ij}f_if_j\right)\frac{\partial}{\partial x^2} +g\left(\frac{\partial}{\partial x^2} , f_1\frac{\partial}{\partial x^2} +f_2 \frac{\partial}{\partial x^1} - g_{12}\sum_{i,j}g^{ij}f_i\frac{\partial}{\partial x^j}\right)\left(\sum_{i,j}g^{ij}f_i\frac{\partial}{\partial x^j}\right), \end{aligned} $$
and the convention that $f_i = \frac{\partial f}{\partial x^i} $ and $f_{ij} = \frac{\partial^2 f}{\partial x^i \partial x^j}$ is adapted.
I know that using the matrix relation $(g_{ij})(g^{ij}) = I$, after expanding fully and rearranging wisely, the result should follow. However, this is so tedious that I don't know if I could manage to find the clue to rearrange it after full expansion. And even if I manage to do so, the process would be too time costly.
Is there any smarter way around? I would really appreciate any help.
On a Riemann Surface equipped with a Riemannian metric $g$ the scalar curvature can always be expressed (up to a constant which I now can't recall) locally as $$s(g)=\Delta_g(\log \sqrt{\det g})$$ where $\det g$ is the local expression of the volume form of the Riemannian metric, if you use isothermal coordinates on the surface.
So now if we perform the conformal change we get $$s(e^f g)=e^{-f}\Delta_g(\log\sqrt{\det g}+f)=e^{-f}(s(g)+\Delta(f)).$$ I guess that you can obtain the same result using geodesic normal coordinates, it would be interesting to see it!