On a Riemannian Manifold $(M,g)$, for any $p \in M$, there exists a normal coordinate system at p, ${U, x_1, ..., x_n}$. I want to integrate over a normal coordinate system with U being a geodesic ball $B_\epsilon(p)$, where $\epsilon>0$ is small enough.
$$\int_{B_\epsilon(p)} \rho(x,p)^{n-\delta} |det(g)|^{1/2} dx_1 ... dx_n,$$
where $\rho{x,p}$ is the geodesic distance between $x$ and $p$. I know that $det(g) = 1 + O(r^2)$ since we are using normal coordinates. However, I would like to know if one can perform some kind of change of variables via polar coordinates to get something like that
$$ \int_{B_\epsilon(p)} \rho(x,p)^{-n+\delta} det(g)^{1/2} dx_1 ... dx_n \leq C\int_0^\epsilon r^{n-\delta} r^{n-1} (1+O(r^2))dr.$$
Basically, I do not understand what happens to $\rho(x,p)^{-n+\delta}$. I understand that normal coordinates allows us to transform the volume form into the usual Euclidean case + O(r^2), but I do not know what happens to the $\rho(x,p)^{-n+\delta}$.