Gauss law from physics I would like to apologise for the question that I am going to pose; and for my curiosity.
Q: How to evaluate the integral below of the function $E(x, y, z) $over $\Omega$ which is the surface of a ball or sphere. Suppose that the radius of sphere is 1; and that it is set at the origin. Define:
$$\vec{E} = \frac{1}{4\pi \varepsilon_0 \varepsilon} \frac{Q} {r^3} \frac{\vec{r}}{r} $$
Gauss law in physics is based on the divergence theorem, stating that
$$\iint_\Omega \vec{E} d\vec{S} =  \frac{Q}{\varepsilon_0}$$
where $\vec{E} $ is the intensity of the electrical field defined by the Coulomb law as $\vec{F} = q\vec{E} $. This question is all about math. From this, we get the definition of DIV by the limit operation:
$$\operatorname{div} \vec{F} = \lim_{\varepsilon \to 0} \frac{1}{V_\varepsilon}\iint_{\Omega_\varepsilon} \vec{F} d\vec{S} $$
 A: As per divergence theorem,
$\displaystyle \iint_\Omega \vec{E} \cdot d\vec{S} =  \iiint_V div(\vec{E}) \ dV$
If we take a single point charge at the origin, $\displaystyle \vec{E} = \frac{1}{4\pi \varepsilon_0} \frac{q} {r^2} \hat{r} \ $ where $\hat{r} = \frac{\vec{r}}{r}$ is a unit normal vector.
Now the outward unit normal vector for a sphere of radius $1$ centered at the origin is $\hat{r}$.
So the surface integral $\iint_\Omega \vec{E} \cdot d\vec{S} = \displaystyle \int_0^{2\pi} \int_0^{\pi} \frac{q}{4\pi \varepsilon_0} \hat{r} \cdot \hat{r} \sin \phi \ d\phi \ d\theta$
$ = \displaystyle \int_0^{2\pi} \int_0^{\pi} \frac{q}{4\pi \varepsilon_0} \sin \phi \ d\phi \ d\theta = \frac{q}{\varepsilon_0}$
If you are applying divergence theorem, in spherical coordinates,
$\displaystyle \iint_\Omega \vec{E} \cdot d\vec{S} =  \frac{q}{4 \pi \varepsilon_0} \iiint_V \nabla \cdot \big(\frac{\vec{r}}{r^2}\big) \ dV = \frac{q}{4 \pi \varepsilon_0} \int_0^{2 \pi} \int_0^{\pi} \int_0^1 \nabla \cdot \big(\frac{\vec{r}}{r^2}\big) \ r^2 \ \sin \phi \ dr \ d\phi \ d\theta$
$\displaystyle  = \frac{q}{\varepsilon_0}$
EDIT: one of the things I forgot to add in my original answer as to how the volume integral comes to that when the vector field has singularity at the origin. In fact the whole contribution to the flux is from the origin and it is zero everywhere else. For that you have to refer to Dirac Delta function. Let me know if you need more details on that.
A: Your formula for $E$ has one too many factors of $1/r$: you can either write $\frac{Q}{4\pi\varepsilon_0r^2}\frac{\vec{r}}{r}$ or $\frac{Q}{4\pi\varepsilon_0r^3}\vec{r}$. In terms of infinitesimal solid angle $d\Omega$, for a radius-$r$ circle centred on the origin$$d\vec{S}=r\vec{r}d\Omega\implies\vec{E}\cdot d\vec{S}=\frac{Q}{4\pi\varepsilon_0}d\Omega\implies\int\vec{E}\cdot d\vec{S}=\frac{Q}{\varepsilon_0}.$$
