How to prove $\oint_{A}\mathbf{E}\cdot d\mathbf{A}=\frac{q}{\epsilon_{0}}$ mathematically and rigorously?
This equation is called “Gauss's Law” in physics. I seek for a rigorous mathematical proof for it (my original physics textbook doesn't provide a rigorous proof).
Now, I want to describe my question in detail.
Let $\mathbf{r}=x\hat{\mathbf{x}}+y\hat{\mathbf{y}}+z\hat{\mathbf{z}}$ be a vector in $\mathbf{R}^{3}$. $\mathbf{E}$ is function from $\mathbf{R}^{3}\backslash{\mathbf{0}}$ to $\mathbf{R}^{3}$ and $\mathbf{E}(\mathbf{r})=\frac{q}{4\pi\epsilon_{0}|\mathbf{r}|^{2}}\hat{\mathbf{r}}$. $A$ is an irregular closed surface enclosing $\mathbf{0}$. Prove $\oint_{A}\mathbf{E}\cdot d\mathbf{A}=\frac{q}{\epsilon_{0}}$.
You can obtain a rigorous proof by first obtaining the result for $\ A=RS^2 ,$ the surface of a sphere of radius $ R\ $ centred on the origin, and then using Gauss's divergence theorem to extend it to an arbitrary surface surrounding the origin. First \begin{align} \oint_{RS^2}\mathbf{E}\cdot d\mathbf{A}&=\frac{q}{4\pi\epsilon_0}\oint_{RS^2}\frac{\hat{\mathbf{r}}\cdot\hat{\mathbf{r}}}{|\mathbf{r}|^2}dA\\ &=\frac{q}{4\pi\epsilon_0R^2}\oint_{RS^2}dA\\ &=\frac{q}{4\pi\epsilon_0}\ , \end{align} because $\ |\mathbf{r}|=R\ $ is constant on $\ RS^2\ ,$ and the surface area of a sphere of radius $\ R\ $ is $\ 4\pi R^2\ .$ Now suppose $\ A\ $ is any compact rectifiable surface surrounding the origin, and therefore at a positive minimum distance from it, and let $\ \delta\stackrel{\text{def}}{=}\min_\limits{\mathbf{r}\in A}\frac{|\mathbf{r}|}{2}\ .$ Let $\ V\ $ be the set enclosed between the surfaces $\ A\ $ and $\ \delta S^2\ .$ Then, since $\ \mathbf{\nabla}\cdot\mathbf{E}=0\ $ throughout $\ V\ ,$ Gauss's theorem gives us \begin{align} 0&=\iiint_V\mathbf{\nabla}\cdot\mathbf{E}\,dV\\ &=\oint_{\partial V}\mathbf{E}\cdot d\mathbf{A}\\ &=\oint_{A\,\cup\,\delta S^2}\mathbf{E}\cdot d\mathbf{A}\\ &=\oint_A\mathbf{E}\cdot d\mathbf{A}+\oint_{\delta S^2}\mathbf{E}\cdot d\mathbf{A}\\ &=\oint_A\mathbf{E}\cdot d\mathbf{A}-\frac{q}{4\pi\epsilon_0}\ , \end{align} the negative sign here arising from the fact that the outward normal to $\ \delta S^2\ $ as the boundary of $\ V\ $ is in the opposite direction to its outward normal as the boundary of a ball centred on the origin.