Why we can use the divergence theorem for electric/gravitational fields if they have singular point?

The divergence theorem is stated (using vector calculus notation)

$$\iiint_{V} \vec{\nabla} \cdot \vec{F}\ dV = \iint_{\partial V} \vec{F} \cdot \vec{n}\ dS$$

For some $C^1$ vector field $\vec{F}$. If we consider, say, the electric field of a point particle :

$$\vec{E} = \alpha \frac{\vec{e_r}}{r^2}$$

for some constant $\alpha$, then $E$ isn't differentiable (or even defined) at $r = 0$.

So here are two ways of dealing with this :

Remove the point

We can simply consider our volume to be $V = B_{0, R} \setminus B_{0, \varepsilon}$, and then consider the limit $\varepsilon \to 0$. We could take $\varepsilon = 0$ directly, but then we'd be integrating a divergence on a set of measure zero, which is also not a great idea.

The boundary is just the two spheres that delimit those regions. Let's see what happens. First, using the formula for the divergence in spherical coordinates (this is simple enough since there is radial symmetry),

\begin{eqnarray} \vec{\nabla} \cdot \vec{E} &=& 0 \end{eqnarray}

This isn't too surprising : outside of $r = 0$, there are no charges. On the other hand, consider the integral on the boundary. The normal vector is simply $\pm\vec{e}_r$, the sign depending on which sphere we're considering (the normal vector has to point away from the volume).

\begin{eqnarray} \iint_{S_{0,R} \cup S_{0, \varepsilon}} (\vec{E} \cdot \vec{n})\ r^2 \sin\theta d\theta d\varphi &=& \iint_{S_{0,R}} (\vec{E} \cdot \vec{n})\ r^2 \sin\theta d\theta d\varphi + \iint_{S_{0, \varepsilon}} (\vec{E} \cdot \vec{n})\ r^2 \sin\theta d\theta d\varphi\\ &=& \iint_{S_{0,R}} E_r\ r^2 \sin\theta d\theta d\varphi - \iint_{S_{0, \varepsilon}} E_r\ r^2 \sin\theta d\theta d\varphi \end{eqnarray}

The $r^2$ of the integral volume measure cancels out with the $r^2$ from the field, and everything only depends on the solid angle of our spheres. As both of them are full spheres, this is just

\begin{eqnarray} 4 \pi - 4 \pi = 0 \end{eqnarray}

All is well : for any neighbourhood that doesn't include the origin, at an arbitrarily small distance $\varepsilon$, the divergence theorem holds. This doesn't help us much with computing though.

Use distributions

Functions aren't well-equipped to deal with point-charges. What we technically have is that those integrals should be equal to the charge $Q$, but as the charge distribution is on $r = 0$, a set of measure zero, any integral should be zero. So instead, we consider the Dirac delta distribution.

A distribution is a linear map from smooth functions of compact support to $\mathbb{R}$. ie, if you have a distribution $f$ and a smooth function of compact support $\phi$, then $f[\phi] \in \mathbb{R}$. There is a map defined thusly : if $f$ is a function, then you can define a distribution from it using the identity

$$f[\phi] = \int f(x) \phi(x) dx$$

In particular, the Dirac delta function is defined by

$$\delta_{\vec{x}}[\phi] = \phi(\vec{x})$$

Applying the Dirac distribution on a function returns the value of that function at a fixed point $\vec{x}$. In our case, this is just $\delta_0$.

Derivatives on distributions (called weak derivatives) are defined by the fact that, if $f$ is a distribution, $f'$ acts on $\phi$ as

$$f'[\phi] = - f[\phi']$$

This corresponds to integrating by part. Now let's consider $F$ as a distribution (or three distributions for every component, this doesn't matter too much). The left-hand of the divergence theorem would roughly correspond to $(\vec{\nabla} \cdot \vec{E})[\phi]$. Now $\vec{E}$ itself is a function, as divergent as it may be. We can find its weak derivative by considering the formula above. Since $E_r$ is a function (We can drop the other components since they are all zero)

\begin{eqnarray} \frac{1}{r^2} \frac{\partial (r^2 E_r)}{\partial r} [\phi] &=& - E_r [\partial_r \phi]\\ &=& -\int \frac{\alpha}{r^2} (\partial_r \phi ) r^2 \sin\theta dr d\theta d\varphi\\ &=& -\int (\partial_r \phi) \sin\theta dr d\theta d\varphi\\ &=& -4\pi [\phi]^\infty_0\\ &=& 4\pi \phi(0) \end{eqnarray}

This is because 1) the fundamental theorem of calculus 2) as the function is of compact support, $\phi(\infty) = 0$. Therefore, the derivative of $E_r[\phi]$ gives $4\pi \phi(0)$, so it is $4\pi \delta_0$.

Of course, usually physicists pretend that $\delta$ is a function, with the property

$$\int \delta(x) f(x) dx = f(0)$$

While this is Bad, it makes everything much simpler to deal with. There is a proper theory to deal with the divergence theorem with distributions, but I'll leave you to read it if you really want the proper proof of this. For the more casual proof, we simply have that

\begin{eqnarray} \iiint_{B_{0,R}} \vec{\nabla} \cdot \vec{E}\ dV &=& \iiint_{V} 4\pi \alpha \delta(x)\ dV\\ &=& 4\pi \alpha\\ \iint_{\partial S_{0,R}} \vec{E} \cdot \vec{n}\ dS &=& \iint_{\partial S_{0,R}} \frac{\alpha}{R^2} \ R^2 \sin\theta d\theta d\varphi\\ &=& 4\pi \alpha \end{eqnarray}

All is well and good.

A third method somewhere between the two would have been to consider the limit of a charge distribution which approaches that of a point charge, such as a normal distribution in the limit $\sigma \to 0$, as distributions can also be considered as the limit of sequences of functions as well.

Proceed as follows.

First accept that the notion of a point charge (with a finite amount of charge but no volume) is itself a mathematical abstraction which may or may not be useful or appropriate in dealing with any given piece of analysis.

Next, treat a small charged region, say a sphere, with finite charge density and finite fields.

Finally, let the size of that region tend to zero and retain whichever are the leading terms in the quantities you are interested in.