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

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.

• ...using polar spherical coordinates centered on $x$, $||x-y||^{-3}(x-y) dy_1dy_2dy_3= r^{-2} \hat{r} r^2 dr d\Omega$, the integral converges! Commented Oct 20, 2015 at 19:33
• if the integral exists finite in polar coordinates (which it does) you can switch back into Cartesian coordinates and are back at your expression. Simply take the polar coordinate version as your definition. Commented Oct 20, 2015 at 21:02
• The integral is of Lebesgue type and it does not depend on the used coordinates. It can also be considered as an improper Riemann integral. Commented Oct 21, 2015 at 6:00
• Your latest edit is almost right but not quite. The domains in x and y need to not intersect, or the function will not be continuous. Shrink the ball in x and leave a buffer shell in the middle. After that the argument works. Commented Oct 25, 2015 at 11:42
• This post has been automatically flagged for numerous edits. The thing is that this is a question & answer site, not an interactive help portal. Please do not continuously modify and extend a question. Edit for clarity yes, but not to increase the scope of the question. You ask a question you get answers. If you have a new question you should use a new post for that. Commented Oct 25, 2015 at 22:31

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.

$$\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.

• Thank you for your very useful answer (and for introducing me to the nice and useful compact notation $\boldsymbol{\Omega}=(\sin\phi\cos\theta,\sin\phi\sin\theta,\cos\phi)$, $\text{d}\Omega=\sin\phi\text{d}\phi\text{d}\theta$)!!! One thing isn't clear to me: how can we see uniform convergence and continuity (as a function of $\boldsymbol{x}$) of the integral if $\rho$ is only piecewise continuous? I $\infty$-ly thank you again! Commented Oct 22, 2015 at 6:54
• If $\rho$ is bounded then $\boldsymbol E$ will be continuous, but if $\rho$ has discontinuities then $\boldsymbol E$ will in general cease to be differentiable there. (For an example, take a uniformly charged sphere.) For these sorts of things, even if you want to be fully rigorous, it is usually better to ask for forgiveness than for permission: that is, do the calculation first and worry about convergence later. You need the calculation, anyways, for the convergence proof. Commented Oct 22, 2015 at 11:55
• If you want precise statements it's probably more fun for you to try it yourself, both defining what you want to prove, deciding on the necessary conditions, and proving it. The standard uniform convergence toolkit is what you need here. Try also your favourite analysis textbook. Commented Oct 22, 2015 at 11:57
• Sorry, it's not quite clear what you're trying to prove - provide a concise statement and I can probably help. Commented Oct 22, 2015 at 14:16
• Don't hold me too tightly to my big-O/little-o notation, which is essentially handwaving it away here. If you really care about those terms it's possible, though. Write a suitably string statement and it should be enough to get you whatever bounds you need for the continuity proof. (That's not where the real action is, though. That's in the uniformity bit.) I'm gonna have to leave you to it on filing in the details though. Commented Oct 23, 2015 at 23:17

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$).

• Since I have been told that I may be right I accept this answer to mark the question as "solved". Nevertheless, if somebody can explain the passage in Emilio's answer that I cannot understand, I will accept his/her answer instead. Commented Nov 16, 2015 at 12:56