Divergence in interaction between induced dipoles Suppose I have two polarizable point-particles 1 and 2 arranged on the z-axis and separated by $r$. Each particle has a dipole moment $\mathbf{p}_i=\alpha \mathbf{E}_i$. If we apply an external field $E_\mathrm{ext}$ also along the z-axis, then by symmetry the moments must point along $z$. A dipole's field along its axis is $\mathbf{E}_p = \frac{2\mathbf{p}}{r^3}$. By symmetry we know the x and y components of each $\mathbf{p}$ are zero, so the problem reduces to two scalar equations:
$$  p_0 = \alpha \left ( \frac{2p_1}{r^3} + E_\mathrm{ext} \right ) $$
$$  p_1 = \alpha \left ( \frac{2p_0}{r^3} + E_\mathrm{ext} \right ) $$
which are solved by
$$ p_0 = p_1 = \alpha E_\mathrm{ext} \left (1 - \frac{2\alpha}{r^3} \right )^{-1} $$.
Not only is the solution singular at a special (finite) separation, but the sign flips so the induced dipoles actually oppose the external field at small separations. You can get out of the singularity by including nonlinear terms in the polarization, but I'm not sure you can get out of the sign flipping. And anyway, there would still be a special point in the solution where the response becomes nonlinear regardless of the strength of the external field. 
It's also weird because we usually assume that we can model polarizable media with the discrete dipole approximation - treating chunks of volume as point dipoles. With strong polarizability, the separation between points in the approximation might be less than the critical value $(2\alpha)^{\frac{1}{3}}$, and then weird things can happen.
So is this physical? Did I make an obvious mistake? I can't see anything wrong with it mathematically. 
 A: Polarizability of a small particle is proportional to volume, i.e. to $a^3$, where $a$ is some characteristic radius. So the weird behavior happens only for distances $r$ comparable to $a$. But such distances necessarily implies that higher multipoles will become important rendering the underlying dipole approximation inaccurate. Using strongly polarizable materials we can somewhat increase the critical $r/a$ but it will also increase the strength of multipoles, so they will remain important for larger $r$ as well.
To conclude this part, the point dipole is a convenient model, but it has limitations with respect to distance to observer ($r\gg a$), while some internal size $a$ always exists even if it is not explicitly specified. If a realistic model is considered, e.g. two slightly separated metallic spheres (with sufficient number of multipoles in numerical calculation) - it may have non-trivial dependence on the separation distance. Some of the observed features may even be qualitatively similar to the weirdness you see in the dipole model, but I have not checked that.
Concerning the discrete dipole approximation (DDA), the problem is alleviated by the fact that many dipoles are used. And to have good numerical accuracy, one usually tries to make dipole size (or, equivalently, their separation) small compared to any characteristic size of the particle. That means that most dipoles will be inside the bulk (locally homogeneous medium), then they are surrounded by other dipoles. And the electrostatic interaction of these dipoles (with similar polarizabilities) largely cancel each other. Numerically, this is exemplified by the convergence of the spectrum of the DDA interaction matrix to that of the integral operator with increasing number of dipoles (for fixed particle). However, using the DDA with only a few dipoles may indeed lead to weird results (non-physical eigenmodes/frequencies). By contrast, the DDA with a single dipole is almost perfect representation for a small sphere.
