It is known that Gauss's law for the electrostatic field $\mathbf{E}$, in the SI, is given by the equation
$$ \int_S \mathbf{E}\cdot \mathrm{d}\mathbf{a}=4\pi k_e Q_{\text{ encl}} \tag{1} $$
where $k_e$ it is the electric constant, $S$ it is the gausssian surface and $Q_{\text{ encl}}$ is the quantity of charge contained inside $S$. What is the general rigorous proof-explanation of the $(1)$?
Being that $\int_S \mathbf{E}\cdot \mathrm{d}\mathbf{a}=\int_S E_x(da)_x+E_y(da)_y+E_z(da)_z$ is there some correlations, for the proof, with the differential quadratic forms?
Gauss's law is equivalent to divergence theorem and since since $\nabla \mathbf{E}=\frac{\rho(\vec P)}{\epsilon_0}$ we have
$$\int_S \mathbf{E}\cdot \mathrm{d}\mathbf{a}=\int_V \nabla \mathbf{E}\, dV=\int_V \frac{\rho(\vec P)}{\epsilon_0} \, dV=\frac{Q_{ \text{ encl} }}{\epsilon_0}=4\pi k_e Q_{ \text{ encl} }$$