The basis is the assumption that the flux density \mathbf{\rho} obeys the Gauss-Green theorem \begin{equation*} \int_V \mathrm{div}\mathbf{\rho} dx = \int_{\partial V} \mathbf{\rho}\mathbf{n}dS, \end{equation*} where n is the unit outer normal field and V is a smooth region with the space where we seek a solution. The Gauss-Green theorem leads to \begin{equation*} \mathrm{div} \mathbf{\rho} = 0. \end{equation*}