The "simplest" classical explanation I know is the van der Waals interaction described by Keesom between two permanent dipoles. Let us consider two permanent dipoles $\vec{p}_1$ (located at $O_1$) and $\vec{p}_2$ located at $O_2$. Their potential energy of interaction is:
\begin{equation}
U(\vec{p}_1,\vec{p}_2,\vec{O_1 O_2}) = -\vec{p}_1\cdot \vec{E}_2 = -\vec{p}_2 \cdot \vec{E}_1
= -\frac{[3(\vec{p}_1\cdot\vec{u})(\vec{p}_2\cdot\vec{u})-\vec{p}_1\cdot \vec{p}_2]}{4\pi \epsilon_0 ||\vec{O_1 O_2}||^3}
\end{equation}
where $\vec{u}= \vec{O_1 O_2}/||\vec{O_1O_2}||$ and $\vec{O_1 O_2}$ is a vector going from $O_1$ to $O_2$.
Now, if the system is subject to thermal fluctuations then the dipoles can change their orientations according to a Boltzmann weight. If the time scale associated to the motion of the particles is much longer than the one associated to the dipole orientations then for each distance $r = ||\vec{O_1O_2}||$ on can trace over all the possible orientations.
The effective interaction between the dipoles is then characterized by the free energy of the system i.e. by:
\begin{equation}
\mathcal{F}(r|p_1,p_2) \equiv -k_B T \ln \int \frac{d\Omega_1 d\Omega_2}{(4\pi)^2} \:e^{-\beta U(\vec{p}_1,\vec{p}_2,\vec{O_1 O_2})}
\end{equation}
where $d\Omega_i$ is the solid angle element associated with the direction of dipole $\vec{p_i}$ and $\beta=1/k_B T$.
This quantity is super difficult to compute. However, if the dipoles are sufficiently far apart, their interaction energy is small and one may expand the exponential in powers of $\beta U$:
\begin{eqnarray}
\mathcal{F}(r|p_1,p_2) \approx -k_B T \ln \int \frac{d\Omega_1 d\Omega_2}{(4\pi)^2} && \times\\
&& \:\left(1 \right.\\
&& -\beta U(\vec{p}_1,\vec{p}_2,\vec{O_1 O_2}) \\
&& \left.+\frac{1}{2}\beta^2 U(\vec{p}_1,\vec{p}_2,\vec{O_1 O_2})^2\right)
\end{eqnarray}
The first integral $\int \frac{d\Omega_1 d\Omega_2}{(4\pi)^2} 1 $ is trivial and equal to $1$. The second integral of the potential is zero by symmetry (each dipole is as likely to be at either a positive or negative angle from $\vec{u}$ and the integral of a scalar product of two independent vectors is zero). The only remaining term is of order 2 in $\beta U$ which decays with $1/r^6$. At the end of the day we end up with something like:
\begin{equation}
\mathcal{F}(r|p_1,p_2) \approx -k_B T \ln \left(1+\beta^2\frac{C}{r^6} \right) \approx -\frac{C}{k_B T}\frac{1}{r^6}
\end{equation}
where
\begin{equation}
C = \int \frac{d\Omega_1 d\Omega_2}{(4\pi)^2}\left( \frac{[3(\vec{p}_1\cdot\vec{u})(\vec{p}_2\cdot\vec{u})-\vec{p}_1\cdot \vec{p}_2]}{4\pi \epsilon_0}\right)^2
\end{equation}
is a number not important in terms of physical insight.
What matters here is that the average value of the dipole-dipole interaction is zero while the first non zero term not to be zero is owing to fluctuations. This is roughly the rule with van der Waals like interactions which is that they are fluctuation driven (the same is true in the quantum case).
Remark: the previous derivation assumes that the interactions travel at infinite speed which is not true. When the limitation of the speed of light is introduced then at large distances, the van der Waals interaction decays as $1/r^7$.
