¿Por qué se cumple la ley de Gauss (en forma integral) para cualquier superficie cerrada? Una forma de verlo es que simplemente se cumple: es una de las ecuaciones de Maxwell, y no hay que "demostrarla". En otras palabras, hay pruebas experimentales al respecto. Mucha gente encuentra más natural expresarla en la forma diferencial ( $\vec{\nabla}·\vec{E} = \rho/\epsilon_0$ ), que es equivalente a la forma integral (mediante el teorema de la divergencia) pero no se refiere a una superficie concreta. La ley de Gauss puede ser motivada por la ley de Coulomb, o incluso derivada de ella en la electrostática, pero la ley de Gauss es más fundamental en el sentido de que se mantiene en la electrodinámica mientras que la ley de Coulomb no.
Aquí hay una derivación algo perezosa de la ley de Gauss a partir de la ley de Coulomb utilizando el teorema de la divergencia que también ilustra cómo $4\pi$ se produce. En general se sale de las matemáticas cuando se trata de situaciones de simetría esférica.
De la ley de Coulomb, $$\vec{E}(\vec{r})=\frac{1}{4\pi\epsilon_0}\iiint\limits_\infty dV' \rho(\vec{r}') \frac{\vec{r} - \vec{r}'}{\lVert \vec{r} - \vec{r}' \rVert^3}.$$ Considere una región $V$ delimitada por la superficie cerrada $\partial V$ . Entonces,
$$ \iint\limits_{\partial V}{\vec{E}(\vec{r})} · d\vec{S} =\frac{1}{4\pi\epsilon_0} \iint\limits_{\partial V}{ \iiint\limits_\infty dV' \rho(\vec{r}') \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3} }.$$
Cambiando el orden de las integrales (que no intentaré justificar), $$ \iint\limits_{\partial V}{\vec{E}(\vec{r})} · d\vec{S} =\frac{1}{4\pi\epsilon_0} \iiint\limits_{\infty}{ dV' \rho(\vec{r}') \iint\limits_{\partial V} { \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3} } }$$
Consideremos ahora la integral $$ I(\vec{r}') = \iint\limits_{\partial V} { \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3}} .$$ Supongamos que $\vec{r}'$ no está encerrado en $\partial V$ para que el denominador sea siempre distinto de cero y el integrando esté siempre bien definido. Por el teorema de la divergencia, $$ I(\vec{r}') = \iiint\limits_V {dV \vec{\nabla} · \frac{\vec{r} - \vec{r}'}{\lVert \vec{r} - \vec{r}' \rVert^3}} . $$ Es un cálculo sencillo demostrar que la divergencia en el integrando desaparece, por lo tanto $I = 0$ .
Si por el contrario $\partial V$ encierra $\vec{r}'$ , dejemos que $S_\delta$ sea una esfera de radio $\delta$ centrado en $\vec{r}'$ y contenida en $\partial V$ . Sea $V - S_\delta$ sea el subconjunto de $V$ en el exterior $S_\delta$ , $\partial(V - S_\delta)$ el límite cerrado de esta región y $\partial S_\delta$ el límite esférico de $S_\delta$ . Podemos escribir $$ I(\vec{r}') = \iint\limits_{\partial(V - S_\delta)} { \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3}} + \iint\limits_{\partial S_\delta} { \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3}}. $$ La primera integral es cero por la misma razón que antes: la divergencia del integrando desaparece en todas partes. La segunda integral puede evaluarse ahora fácilmente utilizando la simetría esférica para obtener
$$ I(\vec{r}') = \frac{\delta\times4\pi\delta^2}{\delta^3} = 4\pi.$$ Aquí es donde el factor de $4\pi$ viene de. Así,
$$ \iint\limits_{\partial V}{\vec{E}(\vec{r})} · d\vec{S} =\frac{1}{4\pi\epsilon_0} \iiint\limits_{\infty}{ dV' \rho(\vec{r}') I(\vec{r}) } = \frac{1}{\epsilon_0} \iiint\limits_{V}{ dV' \rho(\vec{r}') } = \frac{Q}{\epsilon_0} $$ donde la segunda igualdad utiliza el hecho de que el $I(\vec{r}')$ es distinto de cero si y sólo si $V$ contiene $\vec{r}'$ .