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*} $$