  We deal with the electromagnetic waves propagation in the harmonic regime. We
derive the Foldy-Lax approximation of the scattered fields generated by a
cluster of small conductive inhomogeneities arbitrarily distributed in a
bounded domain $\Omega$ of $\mathbb{R}^3$.
  This approximation is valid under a sufficient but general condition on the
number of such inhomogeneities $m$, their maximum radii $\epsilon$ and the
minimum distances between them $\delta$, the form $$(\ln
m)^{\frac{1}{3}}\frac{\epsilon}{\delta} \leq C,$$ where $C$ is a constant
depending only on the Lipschitz characters of the scaled inhomogeneities. In
addition, we provide explicit error estimates of this approximation in terms of
aforementioned parameters, $m, \epsilon, \delta$ but also the used frequencies
$k$ under the Rayleigh regime. Both the far-fields and the near-fields (stated
at a distance $\delta$ to the cluster) are estimated. In particular, for a
moderate number of small inhomogeneities $m$, the derived expansions are valid
in the mesoscale regime where $\delta \sim \epsilon$.
  At the mathematical analysis level and based on integral equation methods, we
prove a priori estimates of the densities in the $L^{2,Div}_t$ spaces instead
of the usual $L^2$ spaces (which are not enough). A key point in such a proof
is a derivation of a particular Helmholtz type decomposition of the densities.
Those estimates allow to obtain the needed qualitative as well as quantitative
estimates while refining the approximation. Finally, to prove the invertibility
of the Foldy-Lax linear algebraic system, we reduce the coercivity inequality
to the one related to the scalar Helmholtz model. As this linear algebraic
system comes from the boundary conditions, such a reduction is not
straightforward.
