I am trying to understand this answer.
The answer is the following:
I'll assume that $V'$ is bounded and that $\mathbf{M'}$ is continuously differentiable. To your primary question:
Is [the] divergence theorem applicable to πΏπ»π of equation (4)?
The answer is, strictly speaking, no. The hypotheses of the divergence theorem include that the vector field in question is continuously differentiable, and $\frac{\mathbf{M'(r)}}{r}$ is not (in general). However, one can still rigorously employ the divergence theorem to yield your expression by removing a ball $B_\epsilon(P)$ from $V'$. Writing $B_\epsilon '(P) \equiv B_\epsilon(P) \cap V' $,
$$\iiint_{V'} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' = \iiint_{V' \backslash B_\epsilon'(P)} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' + \iiint_{B_\epsilon'(P)} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' $$ The divergence Theorem may be applied to the first term for any $\epsilon>0$, and one can show that the second term goes to $0$ in the limit $\epsilon \to 0$ via your chain of equalities in equation $(4)$, which hold $a.e.$ in $B_\epsilon'(P)$, since $\mathbf{M'}$, being continuous on the compact set $\overline{V'}$, is bounded, say by $M>0$.
With some computation, one can show that the integral over $\partial(V' \backslash B_\epsilon '(P))$ arising from the first term approaches the integral over $S' = \partial V'$ in the same limit, assuming $S'$ is smooth in a neighborhood of $P$. To be explicit,
$$ \left( \unicode{x222F}_{\partial V'} - \unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' = \left( \iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)} - \iint_{\partial B_\epsilon '(P) \backslash \partial V'} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' $$
That is, the integrals are the same up to the difference of the integrals over the piece of $\partial B_\epsilon '(P)$ from $\partial V'$ and that from $\partial B_\epsilon(P)$. Since $S'$ is smooth near $P$, it is locally the graph of a smooth function $f(x,y)$ on the tangent plane at $P$ with standard linear coordinates $(x,y)$. So, for $\epsilon$ sufficiently small, we may compute the first integral on the RHS above as an integral in the tangent plane, with area form $\sqrt{1 + f_x^2 + f_y^2}dxdy=\sqrt{1 + f_x^2 + f_y^2}sdsd\theta$ in polar coordinates $(s,\theta)$. The term $\sqrt{1 + f_x^2 + f_y^2}$ is continuous and therefore bounded, say by $C$ (independent of $\epsilon$ for $\epsilon$ sufficiently small), and further $r'(s,\theta)=\sqrt{s^2+f(s,\theta)^2} \geq s$. Hence
$$\left|\iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)}\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS'\right| \leq C \int_0^\epsilon \int_0^{2\pi}|\mathbf{M}(s,\theta)| \frac{s}{\sqrt{s^2+f^2}}dsd\theta \leq 2\pi MC \epsilon$$
Note that the exact integration here may not range over all of $0 \leq s \leq \epsilon$, but the bound still holds. Finally, in the second term on the RHS above, $r'=\epsilon$, so that integral is bounded by $4\pi M\epsilon$, and we have
$$\left| \left( \unicode{x222F}_{\partial V'} - \unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' \right| \leq 2\pi M(2+C)\epsilon $$
which goes to $0$ as $\epsilon \to 0$. Since my first equation holds for all $\epsilon >0$, we may take the limit of both sides as $\epsilon \to 0$ to obtain
$$\iiint_{V'} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' = \unicode{x222F}_{\partial V'} \left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS'$$ As a naive application of the divergence theorem would have suggested.
Note that that "singularity" in the RHS integration does not cause the integral to diverge for essentially the same reason that it didn't in the volume integral: by passing to an appropriately defined set of coordinates in a neighborhood of $P$, we see that the area form gets small fast enough in this neighborhood for the integral to converge, just as the volume form did in your computations.
My misunderstandings are the following:
(1) What does $$ \left( \unicode{x222F}_{\partial V'} - \unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' = \left( \iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)} - \iint_{\partial B_\epsilon '(P) \backslash \partial V'} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' $$ mean? Why are they equal?
My interpretation:
$\unicode{x222F}_{\partial V'}$ means the integral over the whole of the surface of $V'$.
$\text{ }$
$\unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))}$ means the integral over the surface of "$B_\epsilon '(P)$ removed from $V'$.
I have no idea about $\displaystyle\iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)}$ and $\displaystyle\iint_{\partial B_\epsilon '(P) \backslash \partial V'}$ Please explain.