What is the difference between Dirac delta vs removing point from space approach? Let us take for instance E&M, in it when we deal with the failure of divergence theorem to give us the right expression when evaluating the volume integral of divergence of electric field over space, we introduce the concept of divergence of unit radial vector by radius squared being actually equal to dirac delta at origin.
Another approach I saw to solving this discrepancy is to say that we are integrating over a volume with the points which electric charge exist as removed. In such a space with holes, the divergence theorem doesn't hold. This second approach seems conceptually much simpler to me because we don't have to deal with the issue of defining a distribution and such.
What are the rigorous differences between these approaches and are they equivalent?
Refer
 A: The approach you describe appeared in a recent answer of mine, so I'll explain what I meant a bit more thoroughly.
The typical statement of the divergence theorem in $\mathbb R^3$ goes as follows.  Let $B$ be a simple (i.e. non self-intersecting) closed surface which is continuously differentiable (so it has a well-defined outward pointing normal vector $\hat n$ at each point). By the Jordan-Brouwer separation theorem, this surface is the boundary of a bounded, connected interior region $U$.  Given a vector field $\mathbf F$ which is continuously differentiable on $U\cup B$, we have that
$$\int_B \big(\mathbf F \cdot \hat n\big) \mathrm dS = \int_U \mathrm{div}(\mathbf F) \mathrm dV\tag{$\star$}$$
Now consider the vector field $\mathbf F = \hat r/r^2$, with $B$ a simple closed surface enclosing the origin and $U$ the region enclosed.  Naively, one might compute the divergence of $\mathbf F$ and find that $\mathrm{div}(\mathbf F)= 0$, suggesting that both sides of $(\star)$ vanish. However, this is not so, because $\mathbf F$ is not continuously differentiable (or indeed, even defined) at $r=0$.
As a result, $\mathbf F$ does not satisfy the requirements for the divergence theorem to apply.  We could try to fix this problem by removing the origin from $U$ since $\mathbf F$ is continuously differentiable (with vanishing divergence) on $U-\{\mathbf 0\}$, but that still fails to satisfy the conditions of the divergence theorem because $U-\{\mathbf 0\}$ is not the interior of $B$.

This certainly resolves the apparent paradox - that $\mathrm{div}(\mathbf F)$ vanishes on the domain of $\mathbf F$ but the left-hand side of $(\star)$ does not - but it doesn't really tell us how to move forward with our analysis. We're doing physics here, and we'd like to talk about point charges after all.  By extending the divergence theorem to include distribution-valued quantities, we can simply say that the divergence of the electric field due to a point charge is proportional to $\delta$, which matches nicely with the intuition from $\nabla \cdot \mathbf E \propto \rho$.
A: The second paragraph isn't right. The divergence theorem holds provided all its conditions are satisfied. The problem with the point charges is not that there is a hole in the region. Actually, there is no hole in the region (by region I mean where the point charge lives, and that is at the origin of $\Bbb{R}^3$); maybe you meant there is a hole in the domain ($\Bbb{R}^3\setminus\{0\}$) of the fields involved. However, the divergence theorem doesn't really care about the topology of your ambient space, nor whether or not the domain of the fields has any 'holes'. Rather that the fields blow up near the origin.
As for why textbooks introduce $\delta$, there are several reasons

*

*It is convenient mathematically (though depending on the student it may or may not be helpful, depending on how it is presented). Once you learn some basic mechanics of manipulating $\delta$, you can easily arrive at correct answers and usually, ones first intuition is right.

*More significant I think is it allows us to characterize the nature of the singularity. For example, consider on the real line the functions $f_1,f_2,f_3:\Bbb{R}\setminus\{0\}\to\Bbb{R}$, $f_1(x)=\frac{1}{\sqrt{|x|}}$, $f_2(x)=\frac{1}{x^2}$ and $f_3(x)=\frac{1}{x^4}$. These functions satisfy $\lim\limits_{x\to 0}f_1(x)=\lim\limits_{x\to\infty}f_2(x)= \lim\limits_{x\to 0}f_3(x)=\infty$, but clearly $f_2$ is 'more singular' than $f_1$ (e.g $\int_{-1}^1f_1<\infty$ but $\int_{-1}^1f_2=\infty$), and $f_3$ is more singular than $f_2$. Similarly, if we consider a function like $f(x)=x^3+\frac{1}{x^5}$, then we'd say $f$ has a $4^{th}$ order singularity at the origin (more precisely a pole of order $4$). Knowing the nature of singularities of the objects under consideration is very valuable, and since point charges appear all over the theory, it is convenient to be able to characterize singularities by saying "the electric field diverges like that of a point charge".

A: Your premise is wrong: Excluding points of charge from the space(time) might seem "simpler" to you (it's certainly mathematically less involved!), but it's really a non-sensical thing to do: Sure, you could say "well, I guess my electric field just isn't defined on the locations of point charges", but that's physically an extremely strange procedure.
What you're saying in that case is that we have $\nabla \cdot E \propto \rho$ as the rule to determine $E$ from $\rho$ as long as the charge density is some nice differentiable function, and then for point charges we pull out an entirely different rule where $E$ is the superposition of  $\frac{1}{r^2}$ fields centered around the locations of the point charges but not defined there. That's a really bad physical theory: Why are there two sets of rules for two situations when one seems to be the limit of the other - literally, in the sense that you might think about nascent $\delta$-functions as charge densities $\rho_\epsilon$?
This isn't "simple", it's more complicated - instead of a single equation governing the physics you have a single equation for one kind of situations and a weird ad-hoc rule for the other kind! What matters is the simplicity of the physical theory in terms of special cases and assumptions, not how simple you find the math that's going on.
Physical theories are supposed to be parsimonious in their explanations and special cases (cf. Occam's razor) and it's really much simpler in this case to look at nascent $\delta$-functions $\rho_\epsilon$ and say that the equation $\nabla \cdot  E \propto \rho_\epsilon$ still holds as $\nabla \cdot E\propto \delta$ in the distributional sense when $\epsilon\to 0$ and isn't suddenly replaced by a completely different rule.
