Relation of bound charges and polarization In this Wikipedia article the polarization
$$
  \vec P = \rho_\mathrm{b} \vec d
$$
is determined from the density of bound charges $\rho_\mathrm{b}$ and the displacement vector $\vec d$, which denotes the displacement of positive and negative charges in a small volume. If I recall correctly, the volume as a whole has neutral charge.
Now if I am not mistaken, the relation $\rho_\mathrm{b} = - \nabla \vec P$ should hold as well. But how can one verify this?
 A: Imagine some volume filled with point-like charges at the location $\mathbf{R}$, the charge density is then
$$\rho = \sum_{n=1}^N q_n \delta(\mathbf{r} - \mathbf{r}_n - \mathbf{R}) $$
where $\mathbf{r}_n$ is the vector that shows from the centre of the volume to the nth charge. The macroscopic charge density however is an average over this microscopic charge density, so $\rho_m = \langle \rho \rangle$. To average the microscopic charge density use some smearing function $g(\mathbf{r})$ and take the convolution with $\rho$ to get the averaged charge density at the location $\mathbf{r}$
\begin{align*}
\langle \rho\rangle (\mathbf{r}) &= \int \mathrm{d}^3r^{\prime}\, g(\mathbf{r}^{\prime}-\mathbf{r}) \sum_{n=1}^{N}q_n \delta(\mathbf{r}^{\prime}-\mathbf{r}_n - \mathbf{R}) \\
&= \sum_{n=1}^{N} q_n \int \mathrm{d}^3r^{\prime}\, g(\mathbf{r}^{\prime}-\mathbf{r}) \delta(\mathbf{r}^{\prime}-\mathbf{r}_n -\mathbf{R})
\end{align*}
If you now taylor expand the integrand, assuming that $\mathbf{r}<<\mathbf{R}$ you get the following terms
\begin{align*}
\langle \rho \rangle (\mathbf{r}) &= \sum_{n=1}^{N} q_n \int \mathrm{d}^3r^{\prime}\, g(\mathbf{r}^{\prime}-\mathbf{r}) \delta(\mathbf{r}^{\prime}-\mathbf{R}) \\
& + \sum_{n=1}^{N} q_n \int \mathrm{d}^3r^{\prime} g(\mathbf{r}^{\prime}-\mathbf{r}) \Big(-\mathbf{r}_n \cdot \nabla_{\mathbf{r}^{\prime}} \delta(\mathbf{r}^{\prime} - \mathbf{R}) \Big)
\end{align*}
Now the first term gives you a sum over the smeared out charges, which is simply 0 if the Volume has no net charge. The second term however can be further simplified with integration by parts
\begin{align*}
&\quad\ \sum_{n=1}^{N}q_n \mathbf{r}_n \cdot \int \mathrm{d}^3r^{\prime} \delta(\mathbf{r}^{\prime}-\mathbf{R}) \nabla_{\mathbf{r}^{\prime}} g(\mathbf{r}^{\prime}-\mathbf{r}) \\
&= \sum_{n=1}^{N}q_n \mathbf{r}_n \cdot \nabla_{\mathbf{R}}\,g(\mathbf{R}-\mathbf{r}) \\
&= \mathbf{d} \cdot \nabla_{\mathbf{R}}\, g(\mathbf{R}-\mathbf{r})
\end{align*}
Where $\mathbf{d}$ is the microscopic dipole moment. We can now change the variables and get that the averaged bound charge density can be expressed as
\begin{align*}
\langle\rho\rangle(\mathbf{r}) &= -\mathbf{d} \cdot \nabla_{r}\, g(\mathbf{R}-\mathbf{r}) \\
&= - \nabla \cdot \mathbf{d}\, g(\mathbf{R}-\mathbf{r})
\end{align*}
This smeared out dipole vector is the polarisation $\mathbf{P}$, so we get that $\langle \rho \rangle = - \nabla \cdot \mathbf{P}$. And recalling that $\mathbf{d} = \sum_{n=1}^{N} q_n \mathbf{r}_n$, so the polarisation can be written as
\begin{align*}
\mathbf{P}(\mathbf{r}) &= \sum_{n=1}^{N} q_n \mathbf{r}_n g(\mathbf{R}-\mathbf{r}) 
\end{align*}
For simplicity let's look at two charges with opposite sign now so then
\begin{align*}
\mathbf{P}(\mathbf{r}) &= \Big( q \mathbf{r}_1 - q \mathbf{r}_2 \Big) g(\mathbf{R}-\mathbf{r}) \\
&= q \mathbf{D} g(\mathbf{R}-\mathbf{r}) 
\end{align*}
where $\mathbf{D}$ is the displacement vector. Now for the units to work out $g$ has to have units of $1 / V$ so $q g$ is charge per Volume which is the charge density so $\mathbf{P} = \rho \mathbf{D}$
