What is the charge distribution of an Electric Quadrupole?

Question

I'm trying to compute the charge distribution corresponding to the (point) quadrupole moment to gain some intuition (I wanted to know if there can be a quadrupole surface layer like there can be charge e dipole surface layers):

$$ {\newcommand{\r}{\vec{r}}} \phi_2(\r) = \frac{1}{4\pi\varepsilon_0} \frac{3r_hr_k - r^2\delta_{hk}}{r^5}Q_{hk}$$

By Poisson's equation (derived bia Gauss' equation) the corresponding distribution is:

$$ \rho_2 = -\varepsilon_0\nabla^2\phi_2$$

What I end up with is (using vector/matrix notation instead of the indices):

$$ \rho_2 = \frac{3\r\cdot Q\cdot\r}{r^7} + 3r \left((-4\pi\nabla\delta(\r))\cdot Q\cdot \frac{\r}{r^3}\right) + 3r \left(\frac{\r}{r^3}\cdot Q\cdot(-4\pi\nabla\delta(\r))\right) - Tr(Q)\frac{6}{r^5} $$

I can't make much sense of this. Are my calculations wrong? If yes, can anyone provide the correct result? If not, can anyone help understand what the various terms correspond to?

For instance, the charge distribution of a point charge and point dipole are pretty straightforward to understand:

$$ \rho_0(\r) = q\delta(\r) \quad \rho_1(\r) = -\vec{p}\cdot\nabla\delta(\r) $$

The first term represents a charge concentrated in $\r = 0$ of value of $q$. The second term represents a double peak along the direction of $\vec{p}$, with a negative peak in $\r\to 0^-$ and a positive peak in $\r\to 0^+$. This is analogous to the finite dipole, with two point charges of negative value, but with a distance going to zero. Is there a similar interpretation for the electric quadrupole?

EDIT: By looking through my book I've found the expression of the charge distribution:

$ \rho_2(\r) = Q_{hk}\partial_h\partial_k\delta(\r) $

This is also has a nice explanation. For $h\neq k$ there corresponds quadrupoles given by two dipoles in the directions $h$ and $k$, which lie in the plane $r_i = 0$ (with $i\neq h,k$). For instance the term $Q_{12}\delta(z)\partial_1\delta(x)\partial_2\delta(y)$ represents two dipoles along directions $x$ and $y$, which lie in $z=0$ plane. For $h=k$ instead there is a "linear quadrupole", with two dipoles overlapping along one direction. For instante the term $Q_{11}\delta(y)\delta(z)\partial_1\partial_1\delta(x)$ is a "triple peak" along the $x$ direction, with the peak in $\r = 0$ being "double" the ones in $\r\to 0^{\pm}$. For a clearer explanation one would have to use the smooth approximations to the delta function and its derivatives.

However the book did not provide a proof (more of a long calculation) for the result. While I did find various errors in the expression above provided (which assumes that the chain rule still applies to the laplacian) and others I've obtained, I haven't yet gotten to the cited result. Any help would be appreciated.

I'd also add that the derivation be done using the result $\nabla^2(1/r) = \partial_i\partial_i(1/r) = -4\pi\delta(\r)$. I know one may circumvent part of the problem by using Fourier transforms, for which the distributional derivative becomes multiplication, but I'd rather not use this approach.

Answer 1

Starting from the given potential:

$$ {\newcommand{\r}{\vec{r}}} \phi_2(\r) = \frac{1}{4\pi\varepsilon_0} \sum_{ij} \frac{3r_ir_j - r^2\delta_{ij}}{r^5}Q_{ij} \tag{1}$$

First examine \begin{align} \frac{\partial^2}{\partial x_i\partial x_j} \frac{1}{r} =& \frac{\partial}{\partial x_i} \left\{ -\frac{x_j}{r^3}\right\}\\ =& -\frac{\delta_{ij}}{r^3} + 3\frac{x_i x_j}{r^5}\\ =& \frac{3x_i x_j - r^2 \delta_{ij}}{r^5}. \tag{2} \end{align} Note that Eq.(2) resemble the potential part in Eq.(1). Therefore, we can rewrite Eq.(1) in term of Eq.(2): \begin{align} \phi_2(\r) =& \frac{1}{4\pi\varepsilon_0} \sum_{ij}\frac{3r_ir_j - r^2\delta_{ij}}{r^5}Q_{ij}\\ =& \frac{1}{4\pi\varepsilon_0} \sum_{ij} Q_{ij} \frac{\partial^2}{\partial x_i\partial x_j} \frac{1}{r} \end{align}

Now, we may apply the Laplace onto potential $\phi_2$:

\begin{align} \nabla^2 \phi_2(\r) =& \nabla^2 \frac{1}{4\pi\varepsilon_0} \sum_{ij} Q_{ij} \frac{\partial^2}{\partial x_i\partial x_j} \frac{1}{r}\\ =& \frac{1}{4\pi\varepsilon_0} \sum_{ij} Q_{ij} \frac{\partial^2}{\partial x_i\partial x_j} \nabla^2 \left(\frac{1}{r}\right)\\ =& \frac{1}{4\pi\varepsilon_0} \sum_{ij} Q_{ij} \frac{\partial^2}{\partial x_i\partial x_j} \left\{-4\pi\delta^3(\r) \right\}\\ =& -\frac{1}{\varepsilon_0} \sum_{ij} Q_{ij} \frac{\partial^2}{\partial x_i\partial x_j} \delta^3(\r). \end{align}

From Gauss's law, we have the charge density for quadrupole: $$ \rho(\r) = \sum_{ij} Q_{ij} \frac{\partial^2}{\partial x_i\partial x_j} \delta^3(\r). $$