Physical intuition for quadrupole source In his Theory of Vortex Sound M. S. Howe defines sources "mathematically" (i.e. dipole is a source that could be described as a vector and than there is proved that it's equivalent to a two point sources with usual features). The problem arises when it comes to quadrupoles:

A source distribution involving two space derivatives is equivalent to a combination of four monopole sources (whose net volume source strength is zero), and is called a
  quadrupole. A general quadrupole is a source of the form: $$ F = \frac{\partial \ T_{i,j}}{\partial x_i  \partial x_j}$$ in the equation: $$\Box p = F$$

I don't see a link between that and the multipole-expansion and to be honest I am missing the physical intuition on this topic.
Thanks in advance for any hints or clarification!
Note: This question is not really about electrostatics. I have added this tag just because the math is basically the same.
 A: One of the ways of motivating the multipole expansion is as follows: consider a system of charges (more generally, sources, but let's consider the electrostatic case in particular) $q_i$ with position vectors $\mathbf{r}_i$. We want to calculate the electrostatic potential (more generally, whatever field you are interested in) $\phi$ at a point $\mathbf{R}$, such that $r_i \ll R$.
The most convenient case is for all these charges to be at the same location $\mathbf{r}_i = 0$. However, here there is a small deviation from this due to each charge. We can still get a fairly simple approximation to the potential.
For an ordinary function, this is done by means of a Taylor expansion:
$$f(a+x) = f(a) + xf'(a) + \frac{1}{2}x^2f''(a) + ...$$
i.e. essentially a power series in $x$, where $x \ll a$.
We try something similar with the potential, but with the added complication that the variables are vectors. This can still be taken care of however, with an expansion of the schematic type (suppressing the index due to the number of charges):
$$\phi = \sum_{\text{charges}}(\phi_0 + \phi_1\mathbf{r} + \frac{1}{2}\phi_2\mathbf{r}^2+...)$$
We interpret the terms as follows: the first one as if the charges were all at $\mathbf{r} = 0$, the second being the first-order correction for the fact that they are not, and so on. $\phi_n$ contains $n$ derivatives of $\phi$ with respect to the $\mathbf{r}_i$.
Remember that the above expansion is only heuristic. A better way to write it would be as follows:
$$\phi = \sum_{\text{charges}}(\phi_0 + \sum_a (\phi_1)_ar_a + \sum_{a,b}(\phi_2)_{ab}\ r_ar_b +...)$$
(where $a$, $b$ etc. label the components of the vectors).
Each coefficient is more precisely (see, for example, here - the 3D example in this link directly corresponds to the multipole expansion of the electrostatic potential):
$$(\phi_0) = \phi(\mathbf{r}_i = 0)$$
$$(\phi_1)_a = \frac{\partial \phi(\mathbf{r}_i = 0)}{\partial x_a}$$
$$(\phi_2)_{ab} = \frac{\partial^2 \phi(\mathbf{r}_i = 0)}{\partial x_a \partial x_b}$$
Not surprisingly, these are successively called the monopole, dipole, quadrupole... contributions. In particular, let us look at the quadrupole contribution to the potential, which has the form:
$$\sum_{a,b}\frac{\partial^2 f}{\partial x_a \partial x_b}\ r_ar_b$$
Absorbing $r_a$, $r_b$ as constants into $f$, and changing notation, we see that this is essentially of the form:
$$\frac{\partial^2 A_{i,j}}{\partial x_i \partial x_j}$$
This is the quadrupole contribution to the field, given by a source of the type (where we note that the Laplacian acts on the coordinates of the field, $\mathbf{R}$, and not on that of the source distribution):
$$\frac{\partial^2 T_{i,j}}{\partial x_i \partial x_j}$$
The reason behind requiring successive terms in the multipole expansion to resemble some kind of corresponding higher-order derivatives is therefore the Taylor expansion for "nearby sources".
Finally, to (roughly) see how this term corresponds to "four poles", we consider the following case:
$$T_{i,j} = t_{i,j}\delta(\mathbf{r})$$
Now, $\delta(\mathbf{r})$ is a function that is spiked at the origin. We now consider its derivatives. For this, it is easier to consider the one-dimensional case $\delta(x)$. Here, the function is rapidly rising when $x \lt 0$ and rapidly decreasing when $x \gt 0$. The first derivative is then a positive spike just 'left' of the origin, and a negative spike just 'right' of the origin. This can be substantiated by consider the delta function as a limiting case of finite distributions (this article has a relevant image).
Therefore, differentiating a 'delta' spike once splits it into two 'poles' of opposite polarity (a dipole, as in fact used in the reference you quoted). Therefore, a double derivative gives us four 'poles': two of one polarity, and two of another i.e. a quadrupole. For example, with $t_{x, y}$ being the only non-vanishing component, we get a quadrupole in the $x$-$y$ plane.
