I know Gauss' divergence theorem, according to which $$\iiint_D\nabla\cdot\boldsymbol{F}\text{d}x\text{d}y\text{d}z=\iint_{\partial D}\boldsymbol{F}\cdot\boldsymbol{N}_e\text{d}\sigma$$ where $D$ is a solid region satisfying various regularity conditions, whose frontier is $\partial D$, having external unit normal vector $\boldsymbol{N}_e$, and where $\boldsymbol{F}:A\subset\mathbb{R}^3\to\mathbb{R}^3$, $\boldsymbol{F}\in C^1(A)$ with $A$ open and such that $\bar{D}\subset A$.
I read a derivation of Gaus's law in differential form $\nabla\cdot\boldsymbol{E}=\rho/\varepsilon_0$, where $\boldsymbol{E}$ is the electric field, $\varepsilon_0$ is vacuum permittivity and $\rho$ the density of electric charge, from the divergence theorem. In fact, Gauss' laws says that the charge contained in the volume $D$ is $Q=\varepsilon_0\iint_{\partial D}\boldsymbol{E}\cdot\boldsymbol{N}_e\text{d}\sigma$ and, by using the aforesaid theorem, $$\iiint_D\rho \text{d}x\text{d}y\text{d}z=Q=\varepsilon_0\iint_{\partial D}\boldsymbol{E}\cdot\boldsymbol{N}_e\text{d}\sigma=\varepsilon_0\iiint_D\nabla\cdot\boldsymbol{E}\text{d}x\text{d}y\text{d}z$$whence, if we chose $D$ as a parallelepiped, by dividing by the sides of $D$ and by letting the diagonal of $D$ approach $0$, we prove that $\varepsilon_0\nabla\cdot\boldsymbol{E}=\rho$.
But, in this reasoning, we assume that $\boldsymbol{E}$ satisfies the conditions upon $\boldsymbol{F}$ necessary for the divergence theorem to hold.
Nevertheless, if $\rho(\boldsymbol{x}_0)\ne 0$ and $\boldsymbol{x}_0\in D$, how can $\boldsymbol{E}(\boldsymbol{x}_0)=k\int_D\rho(\boldsymbol{x})\|\boldsymbol{x}_0-\boldsymbol{x}\|^{-3}(\boldsymbol{x}_0-\boldsymbol{x})\text{d}x\text{d}y\text{d}z$ (where, as in "standard" notation, $(x,y,z)=\boldsymbol{x}$) exist, finite, how can the field exist finite everywhere and be even of $\boldsymbol{E}\in C^1(\bar{D})$?
Is this one of the cases where physics is not as mathematically rigourous as mathematics itself, as I have been expained here to happen, or am I missing something? I heartily thank you for any answer.
Update (Oct 22 '15) with my trial to prove that $\boldsymbol{E}\in C^k$ using what Dominik, whom I thank again, says in his answer:
Let $\rho\in C^{k}(A)$, $k\ge 0$, where $A$ is an open set such that $\bar{D}\subset A$ and $\forall\boldsymbol{x}\notin \bar{D}\quad\rho(\boldsymbol{x})=0$, with $\bar{D}$ closed and bounded. Without loss of generality I think we can assume $A=\mathbb{R}^3$ because $\rho$ can be extended to such a function. Then the function $\varphi:\mathbb{R}^3\times\mathbb{R}^3\to\mathbb{R}$ defined by $\varphi(\boldsymbol{x},\boldsymbol{x}_0)=-k\frac{\rho(\boldsymbol{x}+\boldsymbol{x}_0)}{\|\boldsymbol{x}\|^3}x$ (and all that I am going to say can be repeated with $y$ and $z$ in place of $x$), is Lebesgue integrable on all $\mathbb{R}^3$ :$$(\boldsymbol{E}(\boldsymbol{x}_0))_x=\int_D\frac{k\rho(\boldsymbol{x})}{\|\boldsymbol{x}_0-\boldsymbol{x}\|^3}(x_0-x)d\mu_\boldsymbol{x}=\int_{D-\boldsymbol{x_0}}\varphi(\boldsymbol{x},\boldsymbol{x}_0)d\mu_\boldsymbol{x}=\int_{\mathbb{R}^3}\varphi d\mu_\boldsymbol{x}.$$
We can see, if $k\ge 1$, that the conditions upon $\rho$ guarantee that the derivatives $\frac{\partial\varphi}{\partial x_0}$, $\frac{\partial\varphi}{\partial y_0}$ and $\frac{\partial\varphi}{\partial z_0}$ all exists, are, for almost all $\boldsymbol{x}\in\mathbb{R}^3$, continuous in $\boldsymbol{x}_0$ on all $\mathbb{R}^3$ and Lebesgue integrable on the same domain, therefore, by using a standard corollary about differentiation under the integral sign to Lebesgue's dominated convergence theorem, we have, for ex. for the derivative with respect to $y_0$, but the same applies to $x_0$ and $z_0$, that $$\bigg(\frac{\partial\boldsymbol{E}}{\partial y_0}\bigg)_x=\frac{\partial }{\partial y_0}\int_{\mathbb{R}^3}\varphi d\mu_\boldsymbol{x}=\int_{\mathbb{R}^3}\frac{\partial \varphi}{\partial y_0} d\mu_\boldsymbol{x}$$Moreover, since the derivative is bounded by the Lebsgue summable function $\boldsymbol{x}\mapsto\frac{x}{\|\boldsymbol{x}\|^3}\max_{\bar{D}}\frac{\partial \rho}{\partial y_0}$, Lebesgue's dominated convergence theorem guarantees that, for any sequence $\{\boldsymbol{x}_{0,n}\}$ such that $\boldsymbol{x}_{0,n}\to\boldsymbol{x}_{0}$, $$\bigg(\frac{\partial\boldsymbol{E}}{\partial y_0}(\boldsymbol{x}_{0,n})\bigg)_x=\int_{\mathbb{R}^3}\frac{\partial \varphi}{\partial y_0}(\boldsymbol{x},\boldsymbol{x}_{0,n}) d\mu_\boldsymbol{x}\to\int_{\mathbb{R}^3}\frac{\partial \varphi}{\partial y_0}(\boldsymbol{x},\boldsymbol{x}_{0}) d\mu_\boldsymbol{x}=\bigg(\frac{\partial\boldsymbol{E}}{\partial y_0}(\boldsymbol{x}_0)\bigg)_x$$
If we repeat the same reasonings made for $\frac{\partial\varphi}{\partial y_0}$ for all the derivatives up to the $k$-th, or just for $\varphi$ if $k=0$, we see that $\boldsymbol{E}\in C^k(\mathbb{R}^3)$.
Best Answer
This is one of the places where we can make things perfectly rigorous if we make certain assumptions on the charge density $\rho$ (and $D$). I will rigorously show you in the following that $\Vert \mathbf{E}(\mathbf{x}_{0}) \Vert < \infty$ for all $\mathbf{x}_{0}$ in the interior of $D$, assuming that $\rho \in (L^{1} \cap L^{\infty})(D)$.
The only dangerous place lies in $\mathbf{x}_{0}$, where the integrand is singular.Pick $\epsilon > 0 $ such that the ball of radius $\epsilon$ around $\mathbf{x}_{0}$, let's call it $B_{\epsilon}(\mathbf{x}_{0})$, is contained in $D$. Then $$\Vert \mathbf{E}(\mathbf{x}_{0}) \Vert \leq \int\limits_{D} \frac{ \vert \rho(\mathbf{x})\vert}{\Vert \mathbf{x}-\mathbf{x}_{0} \Vert^{2}} d\mathbf{x} =\int\limits_{D \backslash B_{\epsilon}(\mathbf{x}_{0})} \frac{ \vert \rho(\mathbf{x})\vert}{\Vert \mathbf{x}-\mathbf{x}_{0} \Vert^{2}} d\mathbf{x}+\int\limits_{B_{\epsilon}(\mathbf{x}_{0})} \frac{ \vert \rho(\mathbf{x})\vert}{\Vert \mathbf{x}-\mathbf{x}_{0} \Vert^{2}} d\mathbf{x}.$$
Using standard inequalities we find that $$\int\limits_{D \backslash B_{\epsilon}(\mathbf{x}_{0})} \frac{ \vert \rho(\mathbf{x})\vert}{\Vert \mathbf{x}-\mathbf{x}_{0} \Vert^{2}} d\mathbf{x} \leq \frac{1}{\epsilon^{2}}\int\limits_{D \backslash B_{\epsilon}(\mathbf{x}_{0})} \vert \rho(\mathbf{x})\vert d\mathbf{x} \leq \frac{1}{\epsilon^{2}} \Vert \rho \Vert_{L^{1}(D)}$$ and that $$\int\limits_{B_{\epsilon}(\mathbf{x}_{0})} \frac{ \vert \rho(\mathbf{x})\vert}{\Vert \mathbf{x}-\mathbf{x}_{0} \Vert^{2}} d\mathbf{x} \leq \Vert \rho \Vert_{L^{\infty}(B_{\epsilon}(\mathbf{x}_{0}))} \int\limits_{B_{\epsilon}(\mathbf{0})} \frac{ 1}{\Vert \mathbf{y}\Vert^{2}} d\mathbf{y} = \Vert \rho \Vert_{L^{\infty}(B_{\epsilon}(\mathbf{x}_{0}))} 4 \pi \epsilon, $$ where the volume element of the ball absorbed the singularity of the integrand.
Therefore we have $\Vert \mathbf{E}(\mathbf{x}_{0}) \Vert \leq \frac{1}{\epsilon^{2}} \Vert \rho \Vert_{L^{1}(D)} + 4 \pi \epsilon \Vert \rho \Vert_{L^{\infty}(D)} < \infty$.
The same analysis works out on the boundary $\partial D$ as well (with minor changes).
I will give you some hints in how one could get continuity/differentiability for $\mathbf{E}$: