Reference: "Theoretische Physik" (2015) by Bartelsmann and others, page 391, equation (11.23).
While deriving the first Maxwell equation based on Coulomb's law, the authors are using the following integral identity:
$$ \begin{align} \oint_{\partial V} \mathbf{d}f \cdot \mathbf{E}(\mathbf{r}) & = k\oint_{\partial V} \mathbf{d}f \cdot \int dV' \rho(\mathbf{r}') \frac{\mathbf{r} - \mathbf{r}'}{|\mathbf{r} - \mathbf{r}'|^3} \\ & = k \int dV' \rho(\mathbf{r}') \oint_{\partial V} \mathbf{d}f \cdot \frac{\mathbf{r} - \mathbf{r}'}{|\mathbf{r} - \mathbf{r}'|^3}, \end{align} $$
where $\mathbf{d}f$ is a surface element, $\mathbf{E}$ is electric field, $k$ is some constant, $\rho$ is volume specific charge density, $\mathbf{r}$ is location of reference charge, and $\mathbf{r}'$ is the integration variable.
Is it true that they first swapped the boundary and volume integral, and then pulled $\rho$ out of the boundary integral? Is it true that $\rho$ is a function of $\mathbf{r}'$, and can't be pulled out of the integral?
Since I could not follow the argument of the answer given by naturallyInconsistent, I made up an example and evaluated it once with integrals swapped, once left in the original order. Here is the result for $\rho = r$, evaluated at $\mathbf{r} = \mathbf{0}$, and integrated over a unit sphere:
Original order of integration:
$$ \Phi(\mathbf{r}) = \int_{\partial V} \left( \int_V \rho \frac{\mathbf{r} - \mathbf{r}'}{||\mathbf{r} - \mathbf{r}'||^3} \text{d}V \right) \cdot \mathbf{n} \text{d}A. $$
How not to do it
In what follows, the integration is done wrong. The misconception here is that the electric field $\mathbf{E}$ can be evaluated at a certain point $\mathbf{r}=\mathbf{0}$, and then integrated over the surface. But this would require the electric field to be constant everywhere. Instead, it is a function of position $\mathbf{E}(\mathbf{r})$.
First, compute $\mathbf{E}(\mathbf{0}):$ $$ \mathbf{E}(\mathbf{0}) = 8 \int_{\varphi=0}^{\pi/2} \int_{\theta=0}^{\pi/2} \int_{r=0}^1 \frac{r}{r^3} \left( \begin{array}{c} r \sin(\theta)\cos(\varphi) \\ r \sin(\theta)\sin(\varphi) \\ r \cos(\theta) \end{array} \right) r^2 \sin(\theta) \text{d}r \text{d}\theta \text{d}\varphi = \left( \begin{array}{c} 2 \pi/3 \\ 4 \pi/3 \\ \pi \end{array} \right) $$
Then, compute the flux:
$$ \Phi(\mathbf{0}) = 8 \int_{\varphi=0}^{\pi/2} \int_{\theta=0}^{\pi/2} \mathbf{E}(\mathbf{0}) \cdot \left( \begin{array}{c} \sin(\theta) \cos(\varphi) \\ \sin(\theta) sin(\varphi) \\ \cos(\theta) \end{array} \right) \sin(\theta) \text{d}\theta \text{d}\varphi = 2 \pi^2 + 16 \pi $$
Now if I swap the order of integration, I'm evaluating the expression
$$ \int_V \int_{\partial V} \left( \rho \frac{\mathbf{r}'}{||\mathbf{r}'||^3} \right) \cdot \mathbf{n} \text{d}A \text{d}V = \int_V \int_{\partial V} \left( \frac{\rho}{||\mathbf{r}'||^2} \right) \mathbf{n} \cdot \mathbf{n} \text{d}A \text{d}V = 8 \pi^2. $$
So there is at least one case where the order of integration may not be changed?
