The integral form of Gauss' law is:

${\displaystyle {\scriptstyle _{S}}}$ ${\displaystyle \mathbf {E} \cdot \mathrm {d} \mathbf {A} }$${\displaystyle ={\frac {Q}{\varepsilon _{0}}}}$

for any closed surface S containing charge Q. By the divergence theorem, this equation is equivalent to:

$${\displaystyle \iiint _{V}\nabla \cdot \mathbf {E} \,\mathrm {d} V={\frac {Q}{\varepsilon _{0}}}}$$

for any volume V containing charge Q. By the relation between charge and charge density, this equation is equivalent to:

$${\displaystyle \iiint _{V}\nabla \cdot \mathbf {E} \,\mathrm {d} V=\iiint _{V}{\frac {\rho }{\varepsilon _{0}}}\,\mathrm {d} V}$$

for any volume V. In order for this equation to be simultaneously true for every possible volume V, it is necessary (and sufficient) for the integrands to be equal everywhere. Therefore, this equation is equivalent to:

$${\displaystyle \nabla \cdot \mathbf {E} ={\frac {\rho }{\varepsilon _{0}}}.}$$

Thus the integral and differential forms are equivalent.

Coulomb's law states that the electric field due to a stationary point charge is:

$${\displaystyle \mathbf {E} (\mathbf {r} )={\frac {q}{4\pi \varepsilon _{0}}}{\frac {\mathbf {e} _{r}}{r^{2}}}}$$

where

• e_r is the radial unit vector,
• r is the radius, |r|,
• ε₀ is the electric constant,
• q is the charge of the particle, which is assumed to be located at the origin.

Using the expression from Coulomb's law, we get the total field at r by using an integral to sum the field at r due to the infinitesimal charge at each other point s in space, to give

$${\displaystyle \mathbf {E} (\mathbf {r} )={\frac {1}{4\pi \varepsilon _{0}}}\int {\frac {\rho (\mathbf {s} )(\mathbf {r} -\mathbf {s} )}{|\mathbf {r} -\mathbf {s} |^{3}}}\,\mathrm {d} ^{3}\mathbf {s} }$$

where ρ is the charge density. If we take the divergence of both sides of this equation with respect to r, and use the known theorem^[9]

$${\displaystyle \nabla \cdot \left({\frac {\mathbf {r} }{|\mathbf {r} |^{3}}}\right)=4\pi \delta (\mathbf {r} )}$$

where δ(r) is the Dirac delta function, the result is

$${\displaystyle \nabla \cdot \mathbf {E} (\mathbf {r} )={\frac {1}{\varepsilon _{0}}}\int \rho (\mathbf {s} )\,\delta (\mathbf {r} -\mathbf {s} )\,\mathrm {d} ^{3}\mathbf {s} }$$

Using the "sifting property" of the Dirac delta function, we arrive at

$${\displaystyle \nabla \cdot \mathbf {E} (\mathbf {r} )={\frac {\rho (\mathbf {r} )}{\varepsilon _{0}}},}$$

which is the differential form of Gauss' law, as desired.

${\displaystyle \Omega \subseteq R^{3}}$

$${\displaystyle \mathbf {E} _{0}(\mathbf {r} )={\frac {1}{4\pi \varepsilon _{0}}}\int _{\Omega }\rho (\mathbf {r} '){\frac {\mathbf {r} -\mathbf {r} '}{\left\|\mathbf {r} -\mathbf {r} '\right\|^{3}}}\mathrm {d} \mathbf {r} '\equiv {\frac {1}{4\pi \varepsilon _{0}}}\int _{\Omega }e(\mathbf {r,\mathbf {r} '} ){\mathrm {d} \mathbf {r} '}}$$

${\displaystyle \rho (\mathbf {r} ')}$

It is true for all ${\displaystyle \mathbf {r} \neq \mathbf {r'} }$ that ${\displaystyle \nabla _{\mathbf {r} }\cdot \mathbf {e} (\mathbf {r,r'} )=0}$.

Consider now a compact set ${\displaystyle V\subseteq R^{3}}$ having a piecewise smooth boundary ${\displaystyle \partial V}$ such that ${\displaystyle \Omega \cap V=\emptyset }$. It follows that ${\displaystyle e(\mathbf {r,\mathbf {r} '} )\in C^{1}(V\times \Omega )}$ and so, for the divergence theorem:

$${\displaystyle \oint _{\partial V}\mathbf {E} _{0}\cdot d\mathbf {S} =\int _{V}\mathbf {\nabla } \cdot \mathbf {E} _{0}\,dV}$$

But because ${\displaystyle e(\mathbf {r,\mathbf {r} '} )\in C^{1}(V\times \Omega )}$,

$${\displaystyle \mathbf {\nabla } \cdot \mathbf {E} _{0}(\mathbf {r} )={\frac {1}{4\pi \varepsilon _{0}}}\int _{\Omega }\nabla _{\mathbf {r} }\cdot e(\mathbf {r,\mathbf {r} '} ){\mathrm {d} \mathbf {r} '}=0}$$

${\displaystyle \Omega \cap V=\emptyset \implies \forall \mathbf {r} \in V\ \ \forall \mathbf {r'} \in \Omega \ \ \ \mathbf {r} \neq \mathbf {r'} }$

${\displaystyle \nabla _{\mathbf {r} }\cdot \mathbf {e} (\mathbf {r,r'} )=0}$

Therefore the flux through a closed surface generated by some charge density outside (the surface) is null.

Now consider ${\displaystyle \mathbf {r} _{0}\in \Omega }$, and ${\displaystyle B_{R}(\mathbf {r} _{0})\subseteq \Omega }$ as the sphere centered in ${\displaystyle \mathbf {r} _{0}}$ having ${\displaystyle R}$ as radius (it exists because ${\displaystyle \Omega }$ is an open set).

Let ${\displaystyle \mathbf {E} _{B_{R}}}$ and ${\displaystyle \mathbf {E} _{C}}$ be the electric field created inside and outside the sphere respectively. Then,

${\displaystyle \mathbf {E} _{B_{R}}={\frac {1}{4\pi \varepsilon _{0}}}\int _{B_{R}(\mathbf {r} _{0})}e(\mathbf {r,\mathbf {r} '} ){\mathrm {d} \mathbf {r} '}}$, ${\displaystyle \mathbf {E} _{C}={\frac {1}{4\pi \varepsilon _{0}}}\int _{\Omega \setminus B_{R}(\mathbf {r} _{0})}e(\mathbf {r,\mathbf {r} '} ){\mathrm {d} \mathbf {r} '}}$ and ${\displaystyle \mathbf {E} _{B_{R}}+\mathbf {E} _{C}=\mathbf {E} _{0}}$

$${\displaystyle \Phi (R)=\oint _{\partial B_{R}(\mathbf {r} _{0})}\mathbf {E} _{0}\cdot d\mathbf {S} =\oint _{\partial B_{R}(\mathbf {r} _{0})}\mathbf {E} _{B_{R}}\cdot d\mathbf {S} +\oint _{\partial B_{R}(\mathbf {r} _{0})}\mathbf {E} _{C}\cdot d\mathbf {S} =\oint _{\partial B_{R}(\mathbf {r} _{0})}\mathbf {E} _{B_{R}}\cdot d\mathbf {S} }$$

The last equality follows by observing that ${\displaystyle (\Omega \setminus B_{R}(\mathbf {r} _{0}))\cap B_{R}(\mathbf {r} _{0})=\emptyset }$, and the argument above.

The RHS is the electric flux generated by a charged sphere, and so:

$${\displaystyle \Phi (R)={\frac {Q(R)}{\varepsilon _{0}}}={\frac {1}{\varepsilon _{0}}}\int _{B_{R}(\mathbf {r} _{0})}\rho (\mathbf {r} '){\mathrm {d} \mathbf {r} '}={\frac {1}{\varepsilon _{0}}}\rho (\mathbf {r} '_{c})|B_{R}(\mathbf {r} _{0})|}$$

${\displaystyle r'_{c}\in \ B_{R}(\mathbf {r} _{0})}$

Where the last equality follows by the mean value theorem for integrals. Using the squeeze theorem and the continuity of ${\displaystyle \rho }$, one arrives at:

$${\displaystyle \mathbf {\nabla } \cdot \mathbf {E} _{0}(\mathbf {r} _{0})=\lim _{R\to 0}{\frac {1}{|B_{R}(\mathbf {r} _{0})|}}\Phi (R)={\frac {1}{\varepsilon _{0}}}\rho (\mathbf {r} _{0})}$$

Taking S in the integral form of Gauss' law to be a spherical surface of radius r, centered at the point charge Q, we have

$${\displaystyle \oint _{S}\mathbf {E} \cdot d\mathbf {A} ={\frac {Q}{\varepsilon _{0}}}}$$

By the assumption of spherical symmetry, the integrand is a constant which can be taken out of the integral. The result is

$${\displaystyle 4\pi r^{2}{\hat {\mathbf {r} }}\cdot \mathbf {E} (\mathbf {r} )={\frac {Q}{\varepsilon _{0}}}}$$

where r̂ is a unit vector pointing radially away from the charge. Again by spherical symmetry, E points in the radial direction, and so we get

$${\displaystyle \mathbf {E} (\mathbf {r} )={\frac {Q}{4\pi \varepsilon _{0}}}{\frac {\hat {\mathbf {r} }}{r^{2}}}}$$

which is essentially equivalent to Coulomb's law. Thus the inverse-square law dependence of the electric field in Coulomb's law follows from Gauss' law.