Showing that $\mathbf{P}(\textbf{r})=\sum_n q_n \mathbf{r}_n \int^1_0 \delta(\mathbf{r}-s\mathbf{r}_n)\mathrm ds$ for a collection of point charges There is one mathematical step that I do not follow in the book Principles of Nano-optics,  page 225. It says that if the charge density is
$\rho=\sum_n q_n \delta(\mathbf{r}-\mathbf{r}_n)$ and the polarisation is $\rho=-\nabla\cdot\mathbf{P}$, then 
$$\mathbf{P}(\textbf{r})=\sum_n q_n \mathbf{r}_n \int^1_0 \delta(\mathbf{r}-s\mathbf{r}_n)\mathrm ds.$$
Could anyone mathematically show how the author arrives at the last expression?
 A: In essence, all you need to do is show that $\mathbf P$ has the correct divergence. Once you do that, it is an appropriate model to underlie your charge density. This is an important distinction: the polarization there is a model, with the operative word being "a": there are multiple equivalent such models, since if you change the origin you will get a different polarization. At heart, you're thinking of your charge density as lots of little head-to-toe dipoles radiating off from the origin and terminating at each of the charges.
The proof can be done in all its vector glory, but I suspect that it will be more illuminating in an explicit coordinate frame. Thus, consider the effect of a single charge $q_0$ placed at $\mathbf r_0=(0,0,z_0)$, for which the polarization reads
$$
\mathbf P_0(\mathbf r) = q_0z_0\hat{\mathbf e}_z \int_0^1\delta(x)\delta(y)\delta(z-sz_0)\mathrm ds.
$$
As I said, the trick is to just jump in and calculate the divergence, which in this instance is just the $z$ derivative of the polarization:
\begin{align}
\nabla\cdot\mathbf P_0(\mathbf r)
& = q_0z_0 \frac{\partial}{\partial z} \int_0^1\delta(x)\delta(y)\delta(z-sz_0)\mathrm ds
\\ & = q_0 \delta(x)\delta(y)\int_0^1z_0\delta'(z-sz_0)\mathrm ds
\\ & = q_0 \delta(x)\delta(y)\int_0^1-\frac{\mathrm d}{\mathrm ds}\left[\delta(z-sz_0)\right]\mathrm ds
\\ & = -q_0 \delta(x)\delta(y)\left[\delta(z-sz_0)\right]_0^1
\\ & = -q_0 \delta(x)\delta(y)\left[\delta(z-z_0)-\delta(z)\right]
\\ & = -q_0\left( \delta(\mathbf r-\mathbf r_0)-\delta(\mathbf r)\right).
\end{align}
As you can see, the calculation relies on the existence of the derivative of the delta function (which is a perfectly legitimate distribution), but it only really pivots off it to pass the derivative from $z$ to $s$, and from there to just kill off the integral. It should also be clear why I went for an explicit coordinate calculation, which is no longer needed after the result is in full vector form, but that this can also be repeated in its full gory details for arbitrary $\mathbf r_0$.
Anyways, once you've done this, you get that the full polarization density, once you've added all the charges, gives you
$$
-\nabla\cdot\mathbf P(\mathbf r) = \sum_n q_n \delta(\mathbf r-\mathbf r_n) - \sum_n q_n \delta(\mathbf r).
$$
The first term is obviously the charge density you were after, but the other term also sticks out like a sore thumb unless the total charge is zero; if that's the case then the model does not work. As I said earlier, the integral over $s$ corresponds to taking lots of little identical dipoles and stacking them head-to-toe, so that the positive charge of one cancels out with the negative charge from its neighbour. This leaves only two charges, at the two ends: one at $\mathbf r_n$, as you wanted, and one at the origin for each and every one of your original charges; if those don't add to zero then you're obviously in trouble.
On the other hand, the whole name of the $\rho=-\nabla\cdot\mathbf P$ game is representing your charge density as resulting from lots of little dipoles, and if your system is not neutral then that's obviously a doomed undertaking anyway.
