Dirac delta function at singularities in spherical coordinates Background Information
Let $\delta^3(\vec x-\vec a)$ represent a point density at $\vec a$. It satisfies
$$
f(\vec a)=\int \delta^3(\vec x-\vec a)f(\vec x)|J(\vec x)|\mathrm d^3\vec x,
$$
where $f(\vec x)$ is an arbitrary function on the space, and $|J(\vec x)|$ is the Jacobian determinant.
On the order hand, the delta function of a vector can be decomposed into the product of several delta functions, of the vector's components,
$$
f(\vec a)=\int f(\vec x)\prod_i\delta^3(x_i-a_i)\mathrm d x_i,\tag{1}
$$
where $x_i,a_i$ are components in the coordinate systems. Comparing the above two equations, we get
$$
\delta^3(\vec x-\vec a) = \frac{1}{|J(\vec a)|}\prod_i\delta(x_i-a_i).
$$
For example,
in Cartesian coordinates, $|J|=1$, so
$$
\delta^3(\vec x-\vec a)=\delta(x-x_a)\delta(y-y_a)\delta(z-z_a).
$$
In spherical coodinates, $|J|=r^2\sin\theta$, and
$$
\delta^3(\vec x-\vec a)=\frac{1}{r^2\sin\theta}\delta(r-r_a)\delta(\theta-\theta_a)
\delta(\varphi-\varphi_a).
$$

Question
The derive above seems illegal when $\vec a$ is a singularity, where $|J|=0$.
The coordinate of a singularity can be indefinite.
For instance, points on $z$-axis (except the origin) in spherical coordinates can be endowed with arbitrary $\varphi$ component. At the origin, even the $\theta$ component is arbitrary, too.
Therefore, when decomposing $\delta^3(\vec r-\vec a)$, there's no corresponding $\delta(\varphi)$ or $\delta(\theta)$.
We must turn back to the initial Equation (1).
For $\vec a$ on the positive half-axis of the $z$-axis,
$$
f(\vec a)=\int f(\vec x)\delta(r-r_a)\delta(\theta-0)\mathrm dr\mathrm d\theta
\mathrm d\varphi
=2\pi\int f(\vec x)\delta(r-r_a)\delta(\theta)\mathrm dr\mathrm d\theta,
$$
so
$$
\delta^3(\vec r-\vec a)=\frac{1}{2\pi r^2\sin\theta}\delta(r-r_a)\delta(\theta).
$$
If $\vec a=0$, then
$$
f(\vec a)=\int f(\vec x)\delta(r-0)\mathrm dr\mathrm d\theta
\mathrm d\varphi
=4\pi\int f(\vec x)\delta(r)\mathrm dr.
$$
Thus
$$
\delta^3(\vec r-\vec a)=\frac{1}{4\pi r^2}\delta(r).
$$
Is there anything wrong? I feel uncomfortable with this result.
This answer (though gained a downvote) states that symmetry results to the difference in the spherical coordinates.
Is there any unified formula similar to Equation (1) that include singularities?
Or do I have to understand this in measure or distribution language? I'm not familiar with that.
 A: The problem that you have is that when you calculate an integral over some volume $\Omega$ in spherical coordinates you are implicitly throwing out the integration region that is not covered by your spherical coordinates. This is no problem because that region is a null-set so it gives a zero contribution to your integral.
What you are using when writing your integral in spherical coordinates is the
$$
\int_{\Omega} f(y) dy = \int_{\Phi^{-1}(\Omega)} f(\Phi(x)) |\det(D\Phi(x))| dx,
$$
where $\Phi$ is the coordinate transformation, that needs to be bijective(!, and differentiable), and $\det(D\Phi(x)) = J(x)$ is the corresponding Jacobi determinant.
Now look at spherical coordinates, where (naively) $\Omega = \mathbb R^3$ and
$$
\Phi((r,\theta,\varphi)) = (r \cos\varphi \sin \theta, r \sin \varphi \sin \theta, r \cos \theta).
$$
Now if want $\Phi$ to be invertible we need to be careful with defining its domain. Naively we have $\theta \in [0,\pi]$, but there we run into problems when $\theta = 0$ or $\pi$, in which case
$$
\Phi^{-1}((0,0,z)) = \{ (z,0,\varphi) ~|~ \varphi \in[0,2\pi) \},
$$
i.e. a point on the $z$-axis has infinite preimages...
What we need to do is say that
$$
\Phi: [0,\infty) \times (0,\pi) \times [0,2\pi ) \rightarrow \mathbb R - \{ (0,0,z ) ~|~ z \in \mathbb R \}.
$$
Now, as you can check, this is bijective, but now $\Omega = \mathbb R - \{ (0,0,z ) ~|~ z \in \mathbb R \}$. And the integral doesn't care that we removed the $z$-axis, since it is a null-set.
But now you see that you get into trouble when you integrate over
$$
f(\vec x) \delta^3(\vec x - \vec a),
$$
where $\vec a \in \{ (0,0,z ) ~|~ z \in \mathbb R \}$. Since $\vec a$ is not in the integration region the integral must give zero
$$
\int_{\Omega} f(\vec x)  \delta^3(\vec x - \vec a) d^3x = 0,
$$
since the delta function vanishes on all of $\Omega$.
But actually it really just depends on how you define the delta function (which tells you that there is something wrong). I like to define it through a Dirac series, for example
$$
\delta_{\epsilon}(\vec x) = \frac{1}{2\pi \epsilon} \exp \Big ( - \frac{\vec x^2}{2 \epsilon} \Big ).
$$
Then
$$
\int_{\Omega} d^3\vec x~ f(\vec x) \delta^3(\vec x - \vec a) := \lim_{\epsilon \rightarrow 0^+} \int_{\Omega} d^3\vec x~ \delta_{\epsilon}(\vec x) f(\vec x) = f(\vec a).
$$
Now the delta function "stretches" into $\Omega$, so you actually get $f(\vec a)$.
If you want to use your formula
$$
\delta( \vec x - \vec a) = \frac{1}{|J(\vec a)|} \delta( \Phi^{-1}(\vec x - \vec a)),
$$
you have to make sure that your cover $\vec a$ with your domain. For this you can just choose different spherical coordinates with respect to some other axis than $z$. Or you simply use symmetry to make $\vec a$ not point in the $z$-direction. In the end it is just an artificial problem that comes from an incorrect choice of coordinates.
