Gauss's law... if the integral defining $\boldsymbol{E}$ diverges?

Question

I have been told (here) that, under particular conditions, the electric field produced by a charge present in space $D$, defined by $$\boldsymbol{E}(\boldsymbol{x})=k\int_D\frac{\rho(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}}(\boldsymbol{x}-\boldsymbol{y})\text{d}y_1\text{d}y_2\text{d}y_3$$where $\boldsymbol{y}=(y_1,y_2,y_3)$ and $\rho$ is the charge density, can converge even if $\boldsymbol{x}\in\bar{D}$ and can even be differentiable.

Although I admit that I have not been able to understand the answers there explaining to me that it may converge - in particular I do not understand whether Riemann triple integrals or Lebesgue integrals are intended - I see that quite strict hypothesis on $\rho$ must be made in order to have $\boldsymbol{E}\in C^1$.

Nevertheless, the elementary electrostatics coursebooks that I have seen use a density function $\rho$ constant on $D$ and zero outside $D$ very often; for example, if we have a sphere $D$ having radius $R$, charge $Q$, and constant charge density $\rho$, my book, Gettys' Physics 2, calculates the field on its surface (as $\frac{Q}{4\pi\varepsilon_0 R}\hat{\boldsymbol{r}}$ where $\hat{\boldsymbol{r}}$ is the radial unitary vector pointing from the centre of the sphere to the point where $\boldsymbol{E}$ is evaluated) by using Gauss's law although, before using Gauss's law, we must know that $\boldsymbol{E}$ and the surface integral $$\int_{\partial D}\boldsymbol{E}\cdot d\boldsymbol{S}$$ exist, while I do not think that $$k\rho\int_D\|\boldsymbol{x}-\boldsymbol{y}\|^{-3}(\boldsymbol{x}-\boldsymbol{y})\text{d}y_1\text{d}y_2\text{d}y_3$$ exists finite [EDIT: I was wrong: Emilio's answer shows that it does] for $\boldsymbol{x}\in\bar{D}$. Am I missing something?

Let us notice that, in order for $\int_{\partial D}\boldsymbol{E}\cdot d\boldsymbol{S}$ to exist, and therefore for Gauss's law to be applicable, the continuity of $\boldsymbol{E}$ would be a sufficient condition and, since$$\boldsymbol{E}_n(\boldsymbol{x}):=k\int_{D\setminus B(\boldsymbol{x},\varepsilon_n)}\frac{\rho(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}}(\boldsymbol{x}-\boldsymbol{y})\text{d}y_1\text{d}y_2\text{d}y_3$$ (where $B(\boldsymbol{x},\varepsilon_n)$ is a ball of radius $\varepsilon_n$ centred in $\boldsymbol{x}$) is continuous as a function of $\boldsymbol{x}$, the uniform convergence of $\boldsymbol{E}_n$ to $\boldsymbol{E}$ as $\varepsilon_n\to 0$ would imply the continuity of $\boldsymbol{E}$.

Therefore I would like to know how to prove that $\boldsymbol{E}$ is continuous and I would be very grateful to anybody confirming or correcting my own answer, or clarifying the points that I say I do not understand in the other existing answer (see update below), or giving an alternative answer (using tools, of course, that I could understand, like multivariate calculus, while I know nothing of measure theory applied to integration in $\mathbb{R}^3$).

---

Update Oct 23 '15 I think that Emilio's proof, whom I heartily thank, uses the Taylor approximation $$\left(1+2\frac{d}{r^2} \hat{\boldsymbol{ n}}\cdot \boldsymbol{r}+\frac{d^2}{r^2}\right)^{-3/2}=1-\frac{3}{2}\left(2\frac{d}{r^2}\boldsymbol{r}\cdot\hat{\boldsymbol{n}}+\frac{d^2}{r^2}\right)+o\left(2\frac{d}{r^2}\boldsymbol{r}\cdot\hat{\boldsymbol{n}}+\frac{d^2}{r^2}\right)$$ as $2\frac{d}{r^2}\boldsymbol{r}\cdot\hat{\boldsymbol{n}}+\frac{d^2}{r^2}\to 0$. Therefore $$\frac{1}{r^3}\left(1+2\frac{d}{r^2} \hat{\boldsymbol{ n}}\cdot \boldsymbol{r}+\frac{d^2}{r^2}\right)^{-3/2}(\boldsymbol{r}+d\hat{\boldsymbol{n}}) $$ $$=\frac{\boldsymbol{r}}{r^3}+\frac{d}{r^3}\hat{\boldsymbol{n}} -\frac{3d\hat{\boldsymbol{ n}}\cdot \boldsymbol{r}}{r^5}\boldsymbol{r}-\frac{3d^2\hat{\boldsymbol{ n}}\cdot \boldsymbol{r}}{r^5}\hat{\boldsymbol{n}}-\frac{3d^2}{2r^5}\boldsymbol{r}-\frac{3d^3}{2r^5}\hat{\boldsymbol{n}}+\frac{o\left(2\frac{d}{r^2}\boldsymbol{r}\cdot\hat{\boldsymbol{n}}+\frac{d^2}{r^2}\right)}{r^3}(\boldsymbol{r}+d\hat{\boldsymbol{n}})$$ $$=\frac{\boldsymbol{r}}{r^3}+d\left( \frac{\hat{\boldsymbol{n}}}{r^3} -\frac{3 \boldsymbol{r}\cdot\hat{\boldsymbol{n}}}{r^5}\boldsymbol{r} \right)+O\left(\frac{d^2}{r^4}\right)+\frac{o\left(2\frac{d}{r^2}\boldsymbol{r}\cdot\hat{\boldsymbol{n}}+\frac{d^2}{r^2}\right)}{r^3}(\boldsymbol{r}+d\hat{\boldsymbol{n}})$$as $2\frac{d}{r^2}\boldsymbol{r}\cdot\hat{\boldsymbol{n}}+\frac{d^2}{r^2}\to 0$. Sadly, I cannot understand why this expression equates $\frac{\boldsymbol r}{r^3}+d\left(\frac{\hat{\boldsymbol n}}{r^3}-3\boldsymbol r\frac{\boldsymbol r\cdot \hat{\boldsymbol n}}{r^5}\right)+O(d^2/r^4)$ (as $d/r^2\to 0$?), and even more, how the inequality $$\left|\frac{\boldsymbol{x}-\boldsymbol{y}}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}}- \frac{\boldsymbol{x}_0-\boldsymbol{y}}{\|\boldsymbol{x}_0-\boldsymbol{y}\|^{3}}\right|\leq 4\frac{||\boldsymbol x-\boldsymbol x_0||}{||\boldsymbol x-\boldsymbol y||^3}+O\left(\frac{\|\boldsymbol{x}-\boldsymbol{x}_0\|^2}{||\boldsymbol x-\boldsymbol y||^4}\right)<4\frac{\delta}{\varepsilon_n^3}+C\frac{\delta^2}{\varepsilon_n^4}$$ (and the denominator $\|\boldsymbol{x}-\boldsymbol{y}\|$ of the first addendum of the right member of the inequality in particular) is derived, while I only can see that $\left\|\frac{\boldsymbol{x}-\boldsymbol{y}}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}}- \frac{\boldsymbol{x}_0-\boldsymbol{y}}{\|\boldsymbol{x}_0-\boldsymbol{y}\|^{3}}\right\|$ $\leq2\frac{\|\boldsymbol{x}_0-\boldsymbol{y}\|}{\|\boldsymbol{x}_0-\boldsymbol{y}\|^3}+O\left(\frac{\|\boldsymbol{x}_-\boldsymbol{x}_0\|^2}{\|\boldsymbol{x}_0-\boldsymbol{y}\|^4}\right)$ $+\left\|\frac{o\left(2\frac{d}{r^2}\boldsymbol{r}\cdot\hat{\boldsymbol{n}}+\frac{d^2}{r^2}\right)}{r^3}(\boldsymbol{r}+d\hat{\boldsymbol{n}})\right\|$ (where the $O$ is intended as $d^2/r^4\to 0$ and the $o$ is intended as $2\frac{d}{r^2}\boldsymbol{r}\cdot\hat{\boldsymbol{n}}+\frac{d^2}{r^2}\to 0$)...

Thank you a lot again!

---

Edit Nov 16 '15: Since I have been told that I may be right I accept my own answer to mark the question as "solved" per se. Nevertheless, since I am also interested in understanding Emilio's proof (which would be my preferred choice, but which it would be dishonest on my part to accept since I cannot even understand its correctness), if somebody can explain the above passages that I cannot understand, I will accept his/her answer instead.

Answer 1

As Valter Moretti has explained, the integral

$$ \newcommand{\bs}[1]{\boldsymbol{#1}} \bs {E}(\bs {x})=k\int_D\frac{\rho(\bs {y})}{\|\bs {x}-\bs {y}\|^{3}}(\bs {x}-\bs {y})\text{d}y_1\text{d}y_2\text{d}y_3 $$

will in general converge as long as the volume charge density $\rho(\boldsymbol y)$ is piecewise continuous and has a finite support. The thing that bugs you is the fact that the Coulomb kernel, $$ \frac{\bs {x}-\bs {y}}{\|\bs {x}-\bs {y}\|^{3}} $$ is of course divergent, and its three components diverge to $\pm \infty$ as $\bs y\to\bs x$. However, this is cancelled out by the fact that there is only a very small amount of charge in the regions where the field-per-charge is large, and this smallness beats out the largeness of the kernel.

This intuition can be made precise by using spherical coordinates centred on $\bs x$. Because the integrand is singular, the integral only makes sense as an improper Riemann integral, i.e. as the limit of Riemann integrals excluding spheres of radius $\varepsilon$ about $\bs x$, as $\varepsilon\to 0$. (It also makes sense as a Lebesgue integral, but both coincide where they are both defined.) For those integration regions the integrand is continuous and you are allowed to make any smooth coordinate transformation you want.

Having made the switch to spherical coordinates in $\bs y=\bs x+r\bs \Omega$, you can evaluate the contribution to $\bs E$ from a spherical shell $S$ between $r_1<r<r_2$, $$ \int_S\rho(\bs x+r\bs \Omega)\frac{\bs \Omega}{r^2}r^2\text{d}r\text{d}\Omega = \int_S\rho(\bs x+r\bs \Omega)\bs \Omega\:\text{d}r\text{d}\Omega $$ so it is bounded. Moreover, as long as $\rho$ is bounded (say, $|\rho(\bs y)|<\rho_0$ in $S$), then the magnitude of the integral over $S$ is also bounded:$$ \int_S \left|\rho(\bs x+r\bs \Omega)\frac{\bs \Omega}{r^2}\right| r^2\text{d}r\text{d}\Omega = \int_S|\rho(\bs x+r\bs \Omega)|\:\text{d}r\text{d}\Omega <\rho_0 r_2. $$ Both of these facts mean that as the radius of the excluded sphere tends to zero, the integral will converge, and it will converge absolutely. Moreover, as long as $\rho(\bs y)$ is bounded for $\bs y$ in a neighbourhood of $\bs x$, the convergence is also uniform, which implies that $\bs E$ is continuous.

---

Addendum

Regarding your edit: The function

$$\bs{E}_n(\bs{x}):=k\int_{D\setminus B(\bs{x}_0,\varepsilon_n)}\frac{\rho(\bs{y})}{\|\bs{x}-\bs{y}\|^{3}}(\bs{x}-\bs{y})\text{d}y_1\text{d}y_2\text{d}y_3$$

is indeed a continuous function of $\bs x$ at $\bs x_0$ - note the altered domain $D\setminus B(\bs{x}_0,\varepsilon_n)$. This is because the integrand is sufficiently continuous at $\bs x_0$. To see this, calculate \begin{align} ||\bs{E}_n(\bs{x})-\bs{E}_n(\bs{x}_0)|| & = \left\| k\int_{D\setminus B(\bs{x}_0,\varepsilon_n)} \left( \frac{\bs{x}-\bs{y}}{\|\bs{x}-\bs{y}\|^{3}} - \frac{\bs{x}_0-\bs{y}}{\|\bs{x}_0-\bs{y}\|^{3}} \right) \rho(\bs{y}) \text{d}y_1\text{d}y_2\text{d}y_3 \right\| \\ & \leq k\int_{D\setminus B(\bs{x}_0,\varepsilon_n)} \left| \frac{\bs{x}-\bs{y}}{\|\bs{x}-\bs{y}\|^{3}} - \frac{\bs{x}_0-\bs{y}}{\|\bs{x}_0-\bs{y}\|^{3}} \right| \,|\rho(\bs{y})|\,\text{d}y_1\text{d}y_2\text{d}y_3 . \end{align} Here you need to restrict $\bs x$ to a closer neighbourhood of $\bs x_0$ than $B(\bs x_0,\varepsilon_n)$ to ensure that it is sufficiently far from all $\bs y$. To fix something, enforce $\bs x\in B(\bs x_0,\delta)\subset B(\bs x_0,\varepsilon_n/2)$, so that $||\bs x-\bs y||>\varepsilon_n/2$ (in addition to $||\bs x_0-\bs y||>\varepsilon_n$) for all $y\in D\setminus B(\bs x_0,\varepsilon_n)$. Under these conditions, the coulomb kernel $$ \frac{\bs{x}-\bs{y}}{\|\bs{x}-\bs{y}\|^{3}} $$ can be expanded by setting $\bs x=\bs x_0+d \hat{\bs n}$ and $\bs r=\bs x_0-\bs y$, $r=||\bs r||$, for convenience, as \begin{align} \frac{\bs{x}-\bs{y}}{\|\bs{x}-\bs{y}\|^{3}} &= \frac{\bs r+d \hat{\bs n}}{\|\bs r+d \hat{\bs n}\|^{3/2}} = \frac{1}{r^3}(\bs r+d \hat{\bs n})\left(1+2\frac{d}{r^2} \hat{\bs n}\cdot \bs r+\frac{d^2}{r^2} \right)^{-3/2} \\ & = \frac{\bs r}{r^3} +d\left(\frac{\hat{\bs n}}{r^3}-3\bs r\frac{\bs r\cdot \hat{\bs n}}{r^5}\right) +O(d^2/r^4). \end{align} Here $d<\delta$ and $r>\varepsilon_n$, which means that the quadratic term is not a problem, and the derivative term is bounded, $$ \left\|\frac{\hat{\bs n}}{r^3}-3\bs r\frac{\bs r\cdot \hat{\bs n}}{r^5}\right\| \leq \frac{1}{r^3}+3\frac{r^2}{r^5}=\frac{4}{r^3}<\frac{4}{\varepsilon_n^3}, $$ so that $$ \left| \frac{\bs{x}-\bs{y}}{\|\bs{x}-\bs{y}\|^{3}} - \frac{\bs{x}_0-\bs{y}}{\|\bs{x}_0-\bs{y}\|^{3}} \right| \leq 4\frac{||\bs x-\bs x_0||}{||\bs x-\bs y||^3} +O\left(\frac{||\bs x-\bs x_0||^2}{||\bs x-\bs y||^4}\right) <4\frac{\delta}{\varepsilon_n^3}+C\frac{\delta^2}{\varepsilon_n^4}. $$ Ignoring the quadratic term from now on, we get for the electric field that \begin{align} ||\bs{E}_n(\bs{x})-\bs{E}_n(\bs{x}_0)|| & < 4k\frac{\delta}{\varepsilon_n^3} \int_{D\setminus B(\bs{x}_0,\varepsilon_n)} \,|\rho(\bs{y})|\,\text{d}y_1\text{d}y_2\text{d}y_3 \end{align} for all $\bs x\in B(\bs x_0,\delta)$ if $\delta<\varepsilon_n/2$, so $||\bs{E}_n(\bs{x})-\bs{E}_n(\bs{x}_0)||$ can be made small by taking $\delta$ small enough. (If that last integral is too big because of contributions from far away, you can treat those separately by splitting it up before you begin approximating stuff.) This means that $\bs E_n$ is continuous at $\bs x_0$.

Finally, to put that in place, you need to prove that the limit $$ \bs{E}(\bs{x}):= \lim_{\varepsilon_n\to 0}\bs{E}_n(\bs{x}) $$ is uniform: that is, that there exists a neighbourhood $B(\bs x_0,d)$ of $\bs x_0$ such that given an arbitrary $\varepsilon>0$, you can enforce $||\bs{E}(\bs{x})-\bs{E}_n(\bs{x})||<\varepsilon$ for all $\bs x\in B(\bs x,d)$ as long as $n>N_\varepsilon$ for a sufficiently big $N_\varepsilon$.

If I understand you correctly you've already done that, and if not then it's a calculation for another day. The techniques, though, are all stuff I talked about above.

Answer 2

I think we can also prove continuity by seeing that the assumptions on $\rho$'s finite support made by Emilio in his answer allow us to use the closure $\overline{D\setminus B(\boldsymbol{x}_0,\varepsilon_n)}$ of $D\setminus B(\boldsymbol{x}_0,\varepsilon_n)$, closed and bounded, as the domain of integration. The function $\boldsymbol{f}:\overline{B\left(\boldsymbol{x}_0,\frac{\varepsilon_n}{2}\right)}\times\overline{D\setminus B(\boldsymbol{x}_0,\varepsilon_n)}\to\mathbb{R}^3$ defined by $\boldsymbol{f(\boldsymbol{x},\boldsymbol{y})}=\frac{\boldsymbol{x}-\boldsymbol{y}}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}}$ can be straightforwardly proved to be continuous (if I'm not wrong) and therefore uniformly continuous on its compact domain. Therefore, for all $\varepsilon>0$ there is a $\delta$ such that, if $\|(\boldsymbol{x},\boldsymbol{y})-(\boldsymbol{x}_0,\boldsymbol{y})\|<\delta$, then $$\|\boldsymbol{E}_n(\boldsymbol{x})-\boldsymbol{E}_n(\boldsymbol{x}_0)\|\leq k\sup_{D}|\rho|\int_{D\setminus B(\boldsymbol{x}_0,\varepsilon_n)}\left\|\frac{\boldsymbol{x}-\boldsymbol{y}}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}} - \frac{\boldsymbol{x}_0-\boldsymbol{y}}{\|\boldsymbol{x}_0-\boldsymbol{y}\|^{3}} \right\|\,\text{d}y_1\text{d}y_2\text{d}y_3$$ $$\leq k\sup_{D}|\rho|\,\varepsilon\,\mu(D\setminus B(\boldsymbol{x}_0,\varepsilon_n))$$ where $\mu$ is the volume.

Then, I would say that we can see that, as $\boldsymbol{x}\to \boldsymbol{x}_0$ (in the domain of integration defining the two different functions called both $\boldsymbol{E}_n$), the $\boldsymbol{E}_n(\boldsymbol{x})$ defined in my original post approaches the $\boldsymbol{E}_n(\boldsymbol{x})$ defined by Emilio in his answer, by using the fact that $\int_{D\setminus B(\boldsymbol{x},\varepsilon_n)}-\int_{D\setminus B(\boldsymbol{x}_0,\varepsilon_n)}=\int_{(D\cap B(\boldsymbol{x}_0,\varepsilon_n))\setminus B(\boldsymbol{x},\varepsilon_n)}$ $-\int_{(D\cap B(\boldsymbol{x},\varepsilon_n))\setminus B(\boldsymbol{x}_0,\varepsilon_n)}$. In fact, for $\boldsymbol{x}\in B(\boldsymbol{x}_0,\delta)$, $\delta\le\frac{\varepsilon_n}{2}$, $$\left\|\int_{D\setminus B(\boldsymbol{x},\varepsilon_n)}\frac{k\rho(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}}(\boldsymbol{x}-\boldsymbol{y})\text{d}y_1\text{d}y_2\text{d}y_3-\int_{D\setminus B(\boldsymbol{x}_0,\varepsilon_n)}\frac{k\rho(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}}(\boldsymbol{x}-\boldsymbol{y})\text{d}y_1\text{d}y_2\text{d}y_3\right\|$$ $$=\Bigg\| \int_\limits{(D\cap B(\boldsymbol{x}_0,\varepsilon_n))\setminus B(\boldsymbol{x},\varepsilon_n)}^{} \frac{k\rho(\boldsymbol{y})(\boldsymbol{x}-\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}}\text{d}y_1\text{d}y_2\text{d}y_3-\int_\limits{(D\cap B(\boldsymbol{x},\varepsilon_n))\setminus B(\boldsymbol{x}_0,\varepsilon_n)}^{}\frac{k\rho(\boldsymbol{y})(\boldsymbol{x}-\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{3}}\text{d}y_1\text{d}y_2\text{d}y_3 \Bigg\|$$$$\le\int_\limits{(D\cap B(\boldsymbol{x}_0,\varepsilon_n))\setminus B(\boldsymbol{x},\varepsilon_n)}^{} \frac{k\sup|\rho|}{\varepsilon_n^2}\text{d}y_1\text{d}y_2\text{d}y_3+\int_\limits{(D\cap B(\boldsymbol{x},\varepsilon_n))\setminus B(\boldsymbol{x}_0,\varepsilon_n)}^{}\frac{k\sup|\rho|}{(\varepsilon_n-\delta)^2}\text{d}y_1\text{d}y_2\text{d}y_3 $$$$\le k\sup|\rho|\left(\frac{\mu((D\cap B(\boldsymbol{x}_0,\varepsilon_n))\setminus B(\boldsymbol{x},\varepsilon_n))}{\varepsilon_n^2}+\frac{\mu((D\cap B(\boldsymbol{x},\varepsilon_n))\setminus B(\boldsymbol{x}_0,\varepsilon_n))}{(\varepsilon_n -\delta)^2}\right)$$where both volumes $\mu$ approaches 0 as $\delta\to 0$.

---

UPDATE: A more straightforward approach, which need the theory of Lebesgue integration is to identify $-k\rho$ with the measurable bounded $\varphi$ and $D$ with the $V$ of the result (yellow paragraph) used in this answer: the components of $\boldsymbol{E}$ are continuous since they are what is called $\frac{\partial\Phi}{\partial x_k}$ there (and the electric potential is $\Phi$).