Could anyone show me the details on the deduction of a discrete object's zeroth, first, second, third and fourth mass moment?
To motivate the moment definitions I provide below, suppose that you have a mass density $\rho(\vec x)$ and suppose that you want to compute the convolution of an arbitrary function $f(\vec x)$ with respect to the mass density:
$$
\bar f \equiv \int d^3x f(\vec x)\rho(\vec x)\;.
$$
If you expand the function $f(\vec x)$ about $\vec x=0$ you arrive at an expression (n.b., repeated indices are summed):
$$
\int d^3x f(0)\rho(\vec x)
+\int d^3x \left.\frac{\partial f}{\partial x^i}\right|_{0}x^i\rho(\vec x)
+\int d^3x \frac{1}{2!}\left.\frac{\partial^2 f}{\partial x^i x^j}\right|_{0}x^i x^j\rho(\vec x)
+\int d^3x \frac{1}{3!}\left.\frac{\partial^3 f}{\partial x^i x^j x^k}\right|_{0}x^i x^j x^k\rho(\vec x)
+\int d^3x \frac{1}{4!}\left.\frac{\partial^4 f}{\partial x^i x^j x^k x^\ell}\right|_{0}x^i x^j x^k x^\ell\rho(\vec x)
+\ldots
$$
$$
=f(0)M_0
+ \left.\frac{\partial f}{\partial x^i}\right|_{0}M^i_1
+ \frac{1}{2!}\left.\frac{\partial^2 f}{\partial x^i x^j}\right|_{0} M^{ij}_2
+ \frac{1}{3!}\left.\frac{\partial^3 f}{\partial x^i x^j x^k}\right|_{0}M^{ijk}_3
+ \frac{1}{4!}\left.\frac{\partial^4 f}{\partial x^i x^j x^k x^\ell}\right|_{0}M^{ijk\ell}_4
+\ldots\;,
$$
where
$$
M_0 \equiv \int d^3x \rho(\vec x)\;, \tag{0}
$$
$$
M^i_1 \equiv \int d^3x x^i \rho(\vec x)\;, \tag{1}
$$
$$
M^{ij}_2 \equiv \int d^3x x^i x^j \rho(\vec x)\;, \tag{2}
$$
$$
M^{ijk}_3 \equiv \int d^3x x^i x^j x^k\rho(\vec x)\;, \tag{3}
$$
$$
M^{ijk\ell}_4 \equiv \int d^3x x^i x^j x^k x^\ell\rho(\vec x)\;, \tag{4}
$$
$$
\ldots
$$
$$
M^{i_1i_2\ldots i_n}_n\equiv \int d^3x x^{i_1} x^{i_2}\ldots x^{i_n} \rho(\vec x)\;, \tag{n}
$$
We can refer to $M_n$ above as the "n-th" moment.
The n-th moment transforms like a rank-n tensor. To see this, consider the density in a different coordinate frame $\tilde x^i = R^{ij}x_j$. In this frame we have a different density function $\tilde \rho (\tilde x)$, but it is related to the initial frame's density function like:
$$
\tilde \rho(\tilde{\vec x}) = \rho(R^{-1}{\tilde{\vec x}})\;.
$$
So, for example:
$$
\tilde{M^{ij}_2} \equiv \int d^3\tilde x \tilde x^{i}\tilde x^{j}\tilde \rho(\tilde x) = \int d^3x \det(R^{-1}) R^{ik}x^k R^{j\ell}x^{\ell}\rho(\vec x)
$$
$$
= R^{ik}R^{j\ell}M^{k\ell}_2\;.
$$
This transformation property is appealing and is one reason why we like to use these moments.
However, these transformation properties alone do not uniquely define the moments proposed above. For example, instead of considering $M_2^{ij}$ we could consider a different moment like
$$
N_2^{ij}\equiv M_2^{ij} - \alpha\delta^{ij}\int d^3x |\vec x|^2\rho(\vec x)\;,
$$
where $\alpha$ is an arbitrary number. For example, if we wanted to make $N_2$ a traceless rank-2 tensor we could set $\alpha = \frac{1}{3}$.
The upshot is that, although the set of $n$ equations for the $M_n$ above provides a useful definition for moments, it is not the only definition one could choose. One's choice is typically dictated by the problem one is trying to solve. For example, another systematic way to define moments is the multipole expansion in terms of spherical harmonics, which tends to be most useful when considering problems that involve an inverse power-law potential like a gravitation potential or an electrostatic potential.