# Effective-medium theory for finite-size aggregates

Charles-Antoine Guérin, Pierre Mallet, Anne Sentenac

► To cite this version:

Charles-Antoine Guérin, Pierre Mallet, Anne Sentenac. Effective-medium theory for finite-size aggregates. Journal of the Optical Society of America. A Optics, Image Science, and Vision, 2006, 23 (2), pp.349-358. 10.1364/JOSAA.23.000349. hal-00080460

HAL Id: hal-00080460

https://hal.science/hal-00080460v1

Submitted on 10 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Effective-medium theory for finite-size aggregates

Charles-Antoine Guérin

Institut Fresnel, UMR CNRS 6133, and Université Paul Cézanne, Faculté de Saint-Jérôme,
F-13397 Marseille cedex 20, France

Pierre Mallet

Institut Fresnel, UMR CNRS 6133, Faculté de Saint-Jérôme, F-13397 Marseille cedex 20, France, and Office
National d'Études et de Recherches Aérospatiales Centre de Toulouse, BP 4025-31055 Toulouse cedex 4, France

Anne Sentenac

Institut Fresnel, UMR CNRS 6133, Faculté de Saint-Jérôme, F-13397 Marseille cedex 20, France

We propose an effective-medium theory for random aggregates of small spherical particles that accounts for the finite size of the embedding volume. The technique is based on the identification of the first two orders of the Born series within a finite volume for the coherent field and the effective field. Although the convergence of the Born series requires a finite volume, the effective constants that are derived through this identification are shown to admit of a large-scale limit. With this approach we recover successively, and in a simple manner, some classical homogenization formulas: the Maxwell Garnett mixing rule, the effective-field approximation, and a finite-size correction to the quasi-crystalline approximation (QCA). The last formula is shown to coincide with the usual low-frequency QCA in the limit of large volumes, while bringing substantial improvements when the dimension of the embedding medium is of the order of the probing wavelength. An application to composite spheres is discussed.

# 1. INTRODUCTION

Effective-medium theory (EMT) is a powerful tool for describing the radiative properties of complex heterogeneous media. More precisely, EMT permits one to identify the average field (or coherent field) propagating inside a random medium, obtained by generating several realizations of the random process, to the field propagating inside a homogeneous medium with an effective permittivity. In several models, the effective permittivity possesses an imaginary part even though the random medium is lossless. This effective absorption describes the attenuation of the coherent field and permits the evaluation of the amount of incoherent scattering.

Since the celebrated Maxwell–Garnett (MG) formula established a century ago,¹ many expressions for the effective permittivity of a composite medium have been proposed. The simplest ones are the so-called mixing rules.² They depend solely on the density of each phase (and their corresponding permittivity) and thus are not able to discriminate inhomogeneous media with different phase repartition when their densities are identical.

More advanced theories, such as the effective-field approximation (EFA), the quasi-crystalline approximation (QCA), and the coherent-potential approximation (CP),³⁻⁸ were proposed subsequently to account more precisely for the statistical distribution of the various phases of the random medium. They are essentially based on renormalization techniques⁹ and deliver an effective permittivity as the solution of a dispersion relation in an infinite medium. In the specific problem of field propagation in an

aggregate, i.e., a medium composed of particles randomly placed in a homogeneous matrix, the effective permittivity depends on the particles' scattering operator and the correlations between the particles' positions. These models, which account for multiple scattering, rely on the strong hypothesis that the random medium is statistically homogeneous, thus necessarily infinite.

The question then arises, if one extracts a small volume of this infinite random medium and illuminates it, will its average radiative properties be adequately described by those of the same volume filled by a homogeneous medium with the permittivity given by the EMT? Apart from its theoretical interest, this question is encountered in various practical situations, such as in attempting to find the effective properties of a sphere with multiple spherical inclusions.¹⁰⁻¹⁴ The main application is atmospheric optics, where water droplets and aerosols may include microcontaminants that modify their scattering and absorption properties.¹⁵ In this case the size of the host medium is typically of the order of the wavelength of the probing light while the inhomogeneities are some order of magnitude smaller. Furthermore, all the numerical experiments that have been conducted to check the accuracy of the EMT formulas usually consider finite aggregates whose typical width is about two wavelengths.¹⁶⁻²⁰

The goal of this paper is to propose an EMT that accounts for the finite size of the random medium, a result that seems to have been underresearched in the literature. We compare the coherent and incoherent scattering cross section of a finite volume of a random medium illu-
<figure><img src="image_2.png" /><figcaption>Fig. 1. Random aggregate of particles is obtained by cutting a finite volume inside an infinite distribution.</figcaption></figure>

minated by a plane wave with the scattering and absorption cross section of a homogeneous object of the same shape. Our method relies on the identification of the Born series at second order with the corresponding scattered fields. An effective permittivity is sought as a perturbation of the MG permittivity. We will show that this effective permittivity depends on the size of the volume and is consistent with the QCA in the limit of large volumes. The work is organized as follows. We give a brief description of the random medium under consideration in Section 2 and then detail the procedure of calculation of the coherent field in Section 3. The core of the paper is the homogenization procedure at different levels in Section 4. In Section 5 we extend the method to the case of composites spheres, that is, spherical host media with random spherical inclusions, and derive a finite-size formula for this case, providing some additional approximations. Section 6 is devoted to a numerical illustration of the main results.

## 2. DESCRIPTION OF THE RANDOM MEDIUM

We consider an assembly of N nonoverlapping identical spherical particles of radius $a$ and permittivity $\epsilon_s$ located at random positions $\mathbf{r}_j$ and immersed in vacuum. We assume that this ensemble of particles is obtained by cutting a test volume $V$ from a spatially homogeneous distribution of particles in the whole space (Fig. 1). The density of particles thus obtained is uniform within the test volume, so that the one-point distribution function is given by

$$ \rho_1(\mathbf{r}) = \frac{1}{V} \Pi_V(\mathbf{r}), \qquad (1) $$

where $\Pi_V$ is the characteristic function of the volume [$\Pi_V(\mathbf{r})=1$ if $r \in V$, $\Pi_V(\mathbf{r})=0$ otherwise]. The n-point distribution functions are usually written in the form

$$ \rho_n(\mathbf{r}_1, \dots, \mathbf{r}_n) = \rho_1(\mathbf{r}_1) \dots \rho_1(\mathbf{r}_n) g_n(\mathbf{r}_1, \dots, \mathbf{r}_n), \quad (2) $$

where $g_n$ expresses the correlation between particles. The case $g_n=1$ corresponds to the case of independent particles. These functions are in general unknown for $n \ge 3$

and one contents oneself with the knowledge of the so-called pair-correlation function (PCF), $g_2$. The PCF has two essential asymptotic properties:

$$ g_2(\mathbf{r}_1, \mathbf{r}_2) = 0, \quad |\mathbf{r}_1 - \mathbf{r}_2| < 2a, \qquad (3) $$

which accounts for the absence of overlap between small spherical particles, and

$$ g_2(\mathbf{r}_1, \mathbf{r}_2) \to 1, \quad |\mathbf{r}_1 - \mathbf{r}_2| \to \infty, \qquad (4) $$

which expresses the mutual decorrelation of particles at large separation. For a spatially homogeneous distribution, the PCF depends only on the separation of particles:

$$ g_2(\mathbf{r}_1, \mathbf{r}_2) = g_2(\mathbf{r}_1 - \mathbf{r}_2). \qquad (5) $$

## 3. CALCULATION OF THE SCATTERED FIELD

### A. Extended Born Series

We now assume that the aggregate is illuminated by an incident monochromatic plane wave

$$ \mathbf{E}_{inc}(\mathbf{r}) = e^{i\mathbf{K}_0 \cdot \mathbf{r}} \mathbf{E}^0. \qquad (6) $$

Here and everywhere an implicit time dependence $\exp(-i\omega t)$ is assumed. We denote $\mathbf{K}_0$ the incident wave vector and $\mathbf{E}^0$ the (constant) polarization vector. The electrical properties of the random medium are described by a random permittivity function $\epsilon(\mathbf{r}) = \epsilon_s$ if $\mathbf{r}$ belongs to a particle and $\epsilon(\mathbf{r})=1$ otherwise. The electric field satisfies the well-known integral equation⁸

$$ \mathbf{E}(\mathbf{r}) = \mathbf{E}_{inc}(\mathbf{r}) + K^2 \int d\mathbf{r}' [\epsilon(\mathbf{r}') - 1] G_0(\mathbf{r} - \mathbf{r}') \mathbf{E}(\mathbf{r}'), \quad (7) $$

where $K=2\pi/\lambda$ is the wavenumber in vacuum and $G_0$ is the free-space dyadic Green's function (l is the unit dyad):

$$ G_0(\mathbf{r}) = (1 + K^{-2} \nabla \cdot \mathbf{\nabla}) \frac{e^{iKr}}{4\pi r}. \qquad (8) $$

The expanded expression of $G_0$ away from its singularity $r=0$ is given by

$$ \left\{ \left( 1 + \frac{i}{Kr} - \frac{1}{(Kr)^2} \right) 1 + \left[ \frac{3}{(Kr)^2} - \frac{3i}{Kr} - 1 \right] \hat{\mathbf{r}} \hat{\mathbf{r}}' \right\} \frac{e^{iKr}}{4\pi r}. \quad (9) $$

The Green's tensor can be decomposed into a singular part and a principal value (PV):

$$ G(\mathbf{r} - \mathbf{r}') = \frac{1}{K^2} \mathbf{L} \delta(\mathbf{r} - \mathbf{r}') + \text{PV} G(\mathbf{r} - \mathbf{r}'), \quad (10) $$

where $\mathbf{L}$ is a constant dyad that depends on the shape of the exclusion domain chosen to define the PV.²¹ For a spherical exclusion domain $\mathbf{L}=-1/3\mathbf{I}$, and we may thus rewrite

$$ \mathbf{E}(\mathbf{r}') = \frac{3}{\epsilon(\mathbf{r}') + 2} \mathbf{E}_{inc}(\mathbf{r}') + 3K^2 \oint d\mathbf{r}'' \frac{\epsilon(\mathbf{r}'') - 1}{\epsilon(\mathbf{r}') + 2} G_0(\mathbf{r}') - \mathbf{r}'' \mathbf{E}(\mathbf{r}''), \quad (11) $$

where we have used the notation $\oint$ for the principal value
integral with spherical exclusion domain around the sin-
gularity:

$$
\oint dr' G(\mathbf{r}-\mathbf{r}') := \lim_{a \to 0} \int_{\mathbf{r}-\mathbf{r}' > a} dr' G(\mathbf{r}-\mathbf{r}'). \quad (12)
$$

Plugging the iteration series of Eq. (11) into Eq. (7) we ob-
tain the following expression for the scattered field E<sub>s</sub>
= **E** - **E**<sub>inc</sub> for **r** outside the test volume:

$$
\mathbf{E}_s(\mathbf{r}) = \int d\mathbf{r}_1 G(\mathbf{r}-\mathbf{r}_1) e^{i\mathbf{K}_0 \cdot \mathbf{r}} \varphi[\beta](\mathbf{r}_1), \quad (13)
$$

where we have introduced a rephased Green’s tensor

$$
G(\mathbf{r}) = G_0(\mathbf{r})e^{-i\mathbf{K}_0 \cdot \mathbf{r}}, \quad (14)
$$

and a source term $\varphi[\beta]$ that is written as a functional ex-
pansion in $\beta$,

$$
\varphi[\beta] = \sum_{n=1}^{\infty} \varphi_n[\beta], \qquad (15)
$$

where $\beta$ is the dimensionless function

$$
\beta(\mathbf{r}) = \frac{\epsilon(\mathbf{r}) - 1}{\epsilon(\mathbf{r}) + 2}. \qquad (16)
$$

The successive terms of this series are found to be

$$
\varphi_1[\beta](\mathbf{r}_1) = 3K^2\beta(\mathbf{r}_1)\mathbf{E}^0,
$$

$$
\varphi_n[\beta](\mathbf{r}_1) = 3K^2\beta(\mathbf{r}_1) \oint \mathrm{d}\mathbf{r}_2 3K^2\beta(\mathbf{r}_2) G(\mathbf{r}_1 - \mathbf{r}_2) \\
\times \oint \dots \oint \mathrm{d}\mathbf{r}_n 3K^2\beta(\mathbf{r}_n) G(\mathbf{r}_{n-1} - \mathbf{r}_n) \mathbf{E}^0, \quad n \ge 2. \tag{17}
$$

The series (15) is expected to converge for small values of
$3K^2\beta$, that is, for small contrast, the geometry held fixed.

B. Coherent Scattered Field

Denote $\mathbf{E}_s(\mathbf{K})$ the far-field amplitude in the remote direc-
tion $\hat{\mathbf{K}}$:

$$
\mathbf{E}_s(\mathbf{r}) \approx \frac{e^{iKr}}{r} \tilde{\mathbf{E}}_s(\mathbf{K}), \quad r \to \infty. \qquad (18)
$$

The ensemble average $\langle \tilde{\mathbf{E}}_s(\mathbf{K}) \rangle$ is the so-called coherent field. Using the far-field expression of the Green's tensor, we obtain easily

$$
\langle \tilde{\mathbf{E}}_s(\mathbf{K}) \rangle = \frac{1}{4\pi} \int d\mathbf{r}_1 (1 - \hat{\mathbf{K}}\hat{\mathbf{K}}) e^{i(\mathbf{K}_0 \cdot \mathbf{K}) \cdot \mathbf{r}_1} \langle \varphi[\beta](\mathbf{r}_1) \rangle, \quad (19)
$$

where the dyad $\hat{\mathbf{K}}\hat{\mathbf{K}}$ is the projector along the direction $\hat{\mathbf{K}}$.
Hence the coherent field is simply related to the ensemble
average $\langle\varphi[\beta](\mathbf{r}_1)\rangle$. It is in general impossible to calculate
analytically this last quantity, for it involves an infinite
series and multiple integrals over the n-point distribution
functions of the random process $\beta(\mathbf{r})$. Therefore one usu-
ally restricts it to the first two terms in the series, which
involve only the statistical mean and the correlation func-

tion of the random process $\beta(\mathbf{r})$, respectively:

$$
\langle \varphi_1[\beta](\mathbf{r}_1) \rangle = 3K^2 \beta_1(\mathbf{r}_1) \mathbf{E}^0,
$$

$$
\langle \varphi_2[\beta](\mathbf{r}_1) \rangle = (3K^2)^2 \int d\mathbf{r}_2 B_2(\mathbf{r}_1, \mathbf{r}_2) G(\mathbf{r}_2 - \mathbf{r}_1) \mathbf{E}^0, \quad (20)
$$

where

$$
B_1(\mathbf{r}) := \langle \beta(\mathbf{r}) \rangle, \qquad (21)
$$

$$
B_2(\mathbf{r}, \mathbf{r}') := \langle \beta(\mathbf{r}) \beta(\mathbf{r}') \rangle. \quad (22)
$$

As we will see, the contribution of the second-order term
$\langle\varphi_2[\beta]\rangle$ is necessary to predict an imaginary part of the ef-
fective permittivity that is not accounted for by the first-
order term. The statistical mean $B_1(\mathbf{r})$ is easily calcu-
lated, since it is merely proportional to the probability of
having a particle at position $\mathbf{r}$:

$$
B_1(\mathbf{r}) = \beta_s \text{Proba}[\mathbf{r} \in \text{particle}] = \beta_s \Pi_V(\mathbf{r}), \quad (23)
$$

where $f$ is the volume density of the scatterers,

$$
f = \frac{4\pi a^3 N}{3V}, \qquad (24)
$$

$$
\beta_s = \frac{\epsilon_s - 1}{\epsilon_s + 2}. \qquad (25)
$$

The correlation $B_2(\mathbf{r}, \mathbf{r}')$ is given by the probability that
the positions $\mathbf{r}$ and $\mathbf{r}'$ lie simultaneously in a particle:

$$
B_2(\mathbf{r}, \mathbf{r}') = \beta_s^2 \text{Proba}[\mathbf{r} \in \text{particle}, \mathbf{r}' \in \text{particle}]. \quad (26)
$$

From the calculations in Appendix A we obtain

$$
\begin{equation}
\begin{split}
B_2(\mathbf{r}, \mathbf{r}') = {}& f \beta_s^2 \Phi \left( \frac{|\mathbf{r} - \mathbf{r}'|}{2a} \right) \Pi_V(\mathbf{r}) \Pi_V(\mathbf{r}') \\
& + f^2 \beta_s^2 \Pi_V(\mathbf{r}) \Pi_V(\mathbf{r}') g_2(\mathbf{r} - \mathbf{r}'),
\end{split}
\tag{27}
\end{equation}
$$

where

$$
\Phi(u) = 1 - \frac{3}{2}u + \frac{1}{2}u^3, \quad \text{if } u \le 1, \quad \Phi(u) = 0 \text{ otherwise,} \tag{28}
$$

is the normalized volume intersected by two spheres with
unit diameter at distance $u$.

This entails the following expression for the coherent
field at second order in the Born series:

$$
\begin{equation}
\begin{split}
\langle \varphi_{1+2}[\beta](\mathbf{r}_1) \rangle ={}& 3K^2 f \beta_s \Pi_V(\mathbf{r}_1) \mathbf{E}^0 + (3K^2 \beta_s)^2 f \\
& \times \int d\mathbf{r}_2 \Pi_V(\mathbf{r}_2) \Phi(\mathbf{r}_{12}/a) G(\mathbf{r}_{12}) \mathbf{E}^0 \\
& + (3K^2 f \beta_s)^2 \Pi_V(\mathbf{r}_1) \\
& \times \int d\mathbf{r}_2 \Pi_V(\mathbf{r}_2) g_2(\mathbf{r}_{12}) G(\mathbf{r}_{12}) \mathbf{E}^0,
\end{split}
\tag{29}
\end{equation}
$$

where $\mathbf{r}_{12}=\mathbf{r}_1-\mathbf{r}_2$ and $\varphi_{1+2}=\varphi_1+\varphi_2$. The second term on
the right-hand side can be simplified by using the isotropy
of $\Phi$. Performing a preliminary angular integration we
may rewrite

$$
K^2 \int d\mathbf{r}_2 \Phi(r_{12}/a) G(\mathbf{r}_{12}) \mathbf{E}^0 = \int_0^a d(Kr) \Phi(r/a) \gamma(Kr) \mathbf{E}^0, \quad (30)
$$

where $\gamma$ is related to the average of the Green's tensor
over a sphere ($S_r$ is the sphere of radius $r$):

$$
K \int_{S_r} dS_r G(\mathbf{r}) \mathbf{E}^0 = \gamma(Kr) \mathbf{E}^0, \quad (31)
$$

The expression of $\gamma$ can be easily derived by using the in-
tegration formulas

$$
\int_{S_r} dS_r e^{-i\mathbf{K} \cdot \mathbf{r}} = 4\pi r^2 \frac{\sin(Kr)}{Kr}, \quad (32)
$$

and

$$
\int_{S_r} dS_r e^{-i\mathbf{K} \cdot \mathbf{r}} \hat{\mathbf{r}} \hat{\mathbf{r}} = 4\pi r^2 \left\{ \left[ \frac{\sin(Kr)}{(Kr)^3} - \frac{\cos(Kr)}{(Kr)^2} \right] 1 + \left[ \frac{\sin(Kr)}{Kr} + 3 \frac{\cos(Kr)}{(Kr)^2} - 3 \frac{\sin(Kr)}{(Kr)^3} \right] \hat{\mathbf{K}} \hat{\mathbf{K}} \right\},
$$

and noting that $\hat{\mathbf{K}}_0\hat{\mathbf{K}}_0\mathbf{E}_0=0$. Explicitly, we obtain

$$
\gamma(u) = \frac{e^{iu}}{u} \left\{ (u^2 + iu - 1)\mathrm{sinc}(u) + \frac{3 - 3iu - u^2}{u^2} [\mathrm{sinc}(u) - \cos(u)] \right\}. \quad (33)
$$

Now for small particles ($Ka \ll 1$) we may use $\Phi \approx 1$ in the
integral (30) and retain the Taylor expansion of $\gamma$ about
zero:

$$
\gamma(u) \approx \frac{11}{15}u + \frac{2}{3}iu^2, \quad u \to 0. \qquad (34)
$$

Altogether this leads to the simplified expression, valid
for $Ka \ll 1$:

$$
\begin{equation}
\begin{split}
\langle \varphi_{1+2}[\beta](\mathbf{r}_1) \rangle = {}& 3K^2 f \beta_s \Pi_V(\mathbf{r}_1) \left[ 1 + \frac{11i}{10}(Ka)^2 \beta_s'' + \frac{2i}{3}(Ka)^3 \beta_s' \right. \\
& \left. \qquad + 3K^2 f \beta_s \int d\mathbf{r}_2 \Pi_V(\mathbf{r}_2) g_2(\mathbf{r}_{12}) G(\mathbf{r}_{12}) \right] \mathbf{E}^0,
\end{split}
\tag{35}
\end{equation}
$$

where $\beta_s'$ and $\beta_s''$ are the real and imaginary parts of $\beta_s$,

$$
\beta_s = \beta_s' + i\beta_s'', \quad (36)
$$

Introducing

$$
\beta_0 = f \beta_s \left[ 1 + \frac{11i}{10} (Ka)^2 \beta_s'' + \frac{2i}{3} \beta_s' (Ka)^3 \right], \quad (37)
$$

and discarding the $O(\beta_s^3)$ terms, we may rewrite this in
the simple final expression

$$
\langle \varphi_{1+2}[\beta] \rangle = 3K^2 \beta_0 [1 + 3K^2 \beta_0 \Pi_V * (g_2 G)] \Pi_V \mathbf{E}^0, \quad (38)
$$

where * stands for the convolution operator.

4. HOMOGENIZATION PROCEDURE

**A. Identification**

We will now seek an effective volume with index $\beta_e$ whose
scattered field coincides with the coherent scattered field
of the random medium; that is, $\phi[\beta_e]=\langle\phi[\beta]\rangle$. Note that
we do not restrict the identification to the forward direc-
tion **(K=K₀)** as is frequently done in homogenization
techniques.²²⁻²⁵ The price to pay for an overall identifica-
tion is the occurrence of a nonconstant permittivity. Find-
ing an exact analytical identification seems a formidable
if not impossible task. It is, however, possible as long as
one restricts the analysis to the first two terms in the
Born series. Equating the first-order terms ($\varphi₁[\beta_e]$
= $\langle\varphi₁[\beta]\rangle$) leads to the equality $\beta_e=f\beta_s$, or

$$
\epsilon_e = \epsilon_{MG} := 1 + 3f \frac{\beta_s}{1 - \beta_s f}, \quad (39)
$$

which is the well-known MG effective permittivity. Note
that this effective permittivity remains real for lossless
particles. An identification within the first two orders
($\varphi_{1+2}[\beta_e] = \langle\varphi_{1+2}[\beta]\rangle$) is expressed by a more complex equa-
tion:

$$
\begin{align}
& \beta_e(\mathbf{r}_1)[1 + 3K^2(\Pi_V \beta_e) * G(\mathbf{r}_1)]\mathbf{E}^0 \nonumber \\
& \qquad = \beta_0[1 + 3K^2\beta_0\Pi_V * (g_2G)(\mathbf{r}_1)]\mathbf{E}^0. \tag{40}
\end{align}
$$

Note that this identity is evidently satisfied with $g_2$
=1 and $\beta_e=\beta_0$ or, in terms of effective permittivity,

$$
\epsilon_e = 1 + 3f \frac{\beta_s}{1 - \beta_s f} \left\{ 1 + \left[ \frac{11i}{10} (Ka)^2 \beta_s'' + \frac{2i}{3} (Ka)^3 \beta_s' \right] \frac{1}{1 - \beta_s f} \right\}. \quad (41)
$$

It is interesting to note that this is again the effective con-
stant $\epsilon_{\text{MG}}$ predicted by the classical MG mixing rule,
apart from the imaginary term. This last term can, how-
ever, be reincorporated in the MG formula by introducing
a radiative correction in the polarizability of small
spheres. This was discussed in detail in Ref. 26. There is
also a close connection with the low-frequency formula of
the EFA⁸ for lossless particles ($\beta_s''=0$) which coincides
analytically to within the $O(\beta_s^2)$ term:

$$
\epsilon_{\text{EFA}} = 1 + 3f \beta_s \left[ 1 + \frac{2i}{3} (Ka)^3 \beta_s \right]. \quad (42)
$$

We will now seek the actual index $\beta_e$ as a multiplicative
correction to $\beta_0$ inside the volume V:

$$
\beta_e(\mathbf{r}) = [1 + \eta(\mathbf{r})]\beta_0. \quad (43)
$$

Inserting Eq. (43) into Eq. (40) leads to the following re-
lation for the perturbation $\eta$:

$$
\eta = 3K^2 \beta_0 \{ [\Pi_V * (g_2 G)] - [1 + \eta] [\Pi_V(1 + \eta) * G] \}, \quad (44)
$$

where we have set $G=\mathbf{E}_0 \cdot (\mathbf{G}\mathbf{E}_0)$. Discarding the quadratic
terms in $\eta$ yields
$$ \eta = 3K^2\beta_0\{\Pi_V * [(g_2 - 1)G]\}. \quad (45) $$

Finally the effective permittivity is given by

$$ \epsilon_e(\mathbf{r}_1) = \frac{1 + 2\beta_e(\mathbf{r}_1)}{1 - \beta_e(\mathbf{r}_1)}. \quad (46) $$

inside the test volume V, with

$$ \beta_e(\mathbf{r}_1) = \beta_0(1 + 3K^2\beta_0\{\Pi_V * [(g_2 - 1)G]\}(\mathbf{r}_1)). \quad (47) $$

Equations (46) and (47) constitute the most general result of this paper. Note that they provide a nonconstant effective permittivity inside the medium. We will now study some important specific cases that render this formula more tractable.

## B. Further Approximations and Limiting Cases

### 1. Large Volumes

Although the Born series is not expected to converge at large volumes, we may take the formal limit $V \to \infty$ in Eq. (47), thereby obtaining a constant value

$$ \beta_e = \beta_0 \left[ 1 + 3K^2\beta_0 \int d\mathbf{r} (g_2 - 1)(\mathbf{r}) G(\mathbf{r}) \right]. \quad (48) $$

For isotropic PCF, the last integral can be turned into the radial integral

$$ \beta_e = \beta_0 \left\{ 1 + 3\beta_0 \int_0^\infty d(Kr) \gamma(Kr) [g_2(r) - 1] \right\}. \quad (49) $$

### 2. Short-Range Correlations

Equation (47) can also be made explicit in the case of short-range and isotropic correlations. Denote by $\xi$ the correlation length, that is, the distance above which one can consider $g_2=1$. If the correlation length is of the order of a few particle diameters, $\xi \sim a$, and if furthermore the particles are small ($Ka \ll 1$), then we may retain the Taylor expansion (34) of $\gamma(u)$ about zero and replace the upper bound by the infinity in the integral. After a simple change of variables, we obtain

$$ \beta_e \approx \beta_0 \left[ 1 + \frac{11i}{5} \beta_0 (Ka)^2 M_1 + 2i \beta_0 (Ka)^3 M_2 \right], \quad (50) $$

where

$$ M_n = \int_0^\infty du [\tilde{g}_2(u) - 1] u^n, \quad n = 1, 2, \quad (51) $$

and $\tilde{g}_2$ is the normalized PCF defined for particles with unit radius, that is $g_2(r)=\tilde{g}_2(r/a)$. This leads to the following expression for the effective permittivity, to within $O[(Ka)^3]$ terms:

$$ \epsilon_{QCA} = 1 + 3f \frac{\beta_s}{1 - \beta_s f} \left[ 1 + \frac{2i}{3} (Ka)^3 \frac{\beta_s'}{1 - \beta_s f} (1 + 3fM_2) + \frac{11i}{10} (Ka)^2 \frac{\beta_s''}{1 - \beta_s f} (1 + 2fM_1) \right], \quad (52) $$

which coincides analytically with the low-frequency QCA approximation⁸ for lossless particles ($\beta_s''=0$). In what follows,

we will therefore refer to formula (47) as finite-size QCA (FS-QCA).

### 3. Spatially Averaged Index

In the general case where the correlation is not short-range nor the volume large the analytical calculations cannot be pushed further, and the effective constant depends in a fairly complicated manner on the position. However, a closely related quantity can be considered instead, namely, the spatially averaged index over the volume:

$$ \bar{\beta}_e = \frac{1}{V} \int d\mathbf{r}_1 \beta_e(\mathbf{r}_1) \Pi_V(\mathbf{r}_1). \quad (53) $$

For a spherical test volume with diameter L, this is easily found to be

$$ \bar{\beta}_e = \beta_0 \left\{ 1 + 3K^2\beta_0 \int d\mathbf{r} \Phi(r/L) [g_2(\mathbf{r}) - 1] G(\mathbf{r}) \right\}, \quad (54) $$

where $\Phi$ is the intersection volume of normalized spheres [Eq. (28)]. For isotropic PCF again, the last expression can be simplified to

$$ \bar{\beta}_e = \beta_0 \left\{ 1 + 3\beta_0 \int_0^L d(Kr) \Phi(r/L) [g_2(r) - 1] \gamma(Kr) \right\}, \quad (55) $$

which is consistent with the large-volume formula (49). The corresponding effective permittivity,

$$ \epsilon_e = \frac{1 + 2\bar{\beta}_e}{1 - \bar{\beta}_e}, \quad (56) $$

will be tested in Section 6.

## 5. EXTENSION TO COMPOSITE SPHERES

The method presented so far has addressed random aggregates confined in a test volume in vacuum. It can, to a certain extent, be adapted to the more physical case of a dielectric, spherical host medium with random inclusions and surrounded by vacuum, typically a rain droplet with inside contaminating particles. The procedure essentially mimics the technique employed in the previous sections, and we will only give the main results and emphasize the differences.

We consider the case of a spherical host medium ($\epsilon_h$) of radius L, immersed in vacuum and containing random spherical inclusions ($\epsilon_s$, radius a). The distribution of small particles is identical to that described in Section 2 that is extracted from a spatially homogeneous, hard-sphere distribution. The dielectric properties of the aggregate are now characterized by the random function $\epsilon(\mathbf{r})$ which can take three different values ($\epsilon_s, \epsilon_h, 1$) according to whether $\mathbf{r}$ is inside one particle, outside the particle but in the test volume V, or outside the test volume, respectively. We denote by $\epsilon_{Mie}(\mathbf{r})$ the permittivity function in the absence of inclusions (homogeneous sphere $\epsilon_h$) and by $\mathbf{E}_{Mie}$ the reference solution corresponding to the total field produced by this homogeneous sphere when illuminated by incident field $\mathbf{E}_{inc}$. We denote by $G_{Mie}$ the dyadic
<figure><img src="image_3.png" /><figcaption>Fig. 2. Pair distribution function for two different densities. There is excellent agreement between Monte Carlo simulations (circles and triangles) and Percus–Yevick hard-sphere distribution functions (solid curves).</figcaption></figure>

Green's function of the sphere, whose expression can be found in Tai.²⁷ The electric field integral equation now reads

$$ \mathbf{E}(\mathbf{r}) = \mathbf{E}_{\text{Mie}}(\mathbf{r}) + K^2 \int d\mathbf{r}' [\epsilon(\mathbf{r}') - \epsilon_{\text{Mie}}(\mathbf{r}')] G_{\text{Mie}}(\mathbf{r}, \mathbf{r}') \mathbf{E}(\mathbf{r}'). \quad (57) $$

The calculations presented in Sections 3 and 4 can be reproduced step by step if one makes the following two simplifying assumptions. First, the reference field inside the sphere is given by the first (extended) Born approximation:

$$ \mathbf{E}_{\text{Mie}}(\mathbf{r}) \approx \frac{3}{2 + \epsilon_h} \mathbf{E}_{\text{inc}}(\mathbf{r}). \quad (58) $$

Second, the Green's function of the sphere can be replaced by the free-space Green's function in the host medium close to the diagonal:

$$ G_{\text{Mie}}(\mathbf{r}, \mathbf{r}') \approx G_h(\mathbf{r} - \mathbf{r}'), \quad \text{if} \quad \|\mathbf{r} - \mathbf{r}'\| \ll \lambda, \quad (59) $$

where $G_h$ is given by Eq. (8) with $K$ replaced by the wave-number $K_h$ in the medium ($K_h^2 = \epsilon_h K^2$). This second approximation is justified since $G_{\text{Mie}}$ can actually be decomposed into a homogeneous-space Green's function and an additional nonsingular part that contains the multiple reflections and transmissions at the boundary of the sphere. Hence the behavior of $G_{\text{Mie}}(\mathbf{r}, \mathbf{r}')$ for close points $(\mathbf{r}, \mathbf{r}')$ is imposed by $G_h(\mathbf{r} - \mathbf{r}')$.

Following the same identification procedure as in Section 4, we recover successively the MG mixing rule in the host medium,

$$ \epsilon_c = \epsilon_h + 3\epsilon_h f \frac{\beta_{sh}}{\epsilon_h - \beta_{sf}}, \quad (60) $$

where

$$ \beta_{sh} = \frac{\epsilon_s - \epsilon_h}{\epsilon_s + 2\epsilon_h}; \quad (61) $$

the EFA in the host medium,

$$ \epsilon_{\text{EFA}} = \epsilon_h + 3\epsilon_h f \beta_{sh} \left[ 1 + \frac{2i}{3} (K_h a)^3 \beta_{sh} \right]; \quad (62) $$

and a finite-size QCA

$$ \epsilon_c(\mathbf{r}_1) = \frac{\epsilon_h + 2\beta_c(\mathbf{r}_1)}{\epsilon_h - \beta_c(\mathbf{r}_1)} \quad (63) $$

with

$$ \beta_c(\mathbf{r}_1) = \beta_{0h}(1 + 3K_h^2 \beta_{0h} \{\Pi_V * [(g_2 - 1)G_h]\}(\mathbf{r}_1)), \quad (64) $$

where $G_h = e^{-iK_0 \cdot \mathbf{r}} \mathbf{E}_0 \cdot (G_h \mathbf{E}_0)$ and $\beta_{0h}$ is defined by Eq. (37) with the replacements $\beta_s \to \beta_{sh}$ and $K \to K_h$. A more tractable expression for Eq. (64) can again be obtained by considering the spatially averaged index $\tilde{\beta}_e$ and specializing to isotropic PCF $g_2$. In that case formula (55) is recovered with the replacements $K \to K_h$ and $\beta_0 \to \beta_{0h}$.

# 6. NUMERICAL EXPERIMENTS

We will now give numerical illustrations of the previous theoretical developments. To make rigorous computations possible, we will restrict ourselves to an ensemble of small lossless particles confined in a test volume in vacuum, that is, without an embedding matrix.

## A. Random Medium Generation

For a numerical validation of the effective-medium formulas, it is very important to generate a random medium whose statistical properties, or at least the one- and two-point distribution functions, are perfectly controlled. One such medium is the hard-sphere model that is commonly used to model a system with nonpenetrable particles. Particles are placed at random in a given volume with the only constraint that they not overlap. The simplest generation method is a sequential deposition algorithm with systematic retrial in case of overlap. This technique becomes prohibitive at high density because of the high rejection rate. An efficient algorithm at higher densities is the so-called Metropolis et al.²⁸ shuffling algorithm. We refer to Refs. 17, 18, and 29 for the detailed utilization of this algorithm in the context of hard-sphere medium generation. The shuffling process is performed in a cubical region with periodic boundary conditions. Only the scatterers lying inside the inscribed spherical volume are finally retained. The PCF of the resulting medium is accurately described by the Percus–Yevick (PY) integral equation,³⁰ for which an exact solution is known.³¹,³² A comparison between the numerical PCF based on Monte Carlo simulations and the PY PCF is shown in Fig. 2.

## B. Field Computations

A rigorous computation of the field scattered by an assembly of spheres requires the involved formalism of multiple Mie scattering³³ or the method of moments¹⁸ and is numerically very demanding, especially when it comes to Monte Carlo simulations. For electrically small particles, however, the computation can be greatly simplified by making a dipolar approximation for the scatterers. This is possible for small size parameters $\sqrt{\epsilon}Ka$, where $\epsilon$ and $a$
are the permittivity and radius of the particle, respectively. In this case the field scattered by the aggregate is given merely by

$$ \mathbf{E}_s(\mathbf{r}) = \alpha_s K^2 \sum_{j=1}^{N} G_0(\mathbf{r} - \mathbf{r}_j) \mathbf{E}_j, \quad (65) $$

where $\mathbf{E}_j$ is the exciting field on the jth dipole and $\alpha_s$ its polarizability. The exciting fields are mutually related by the system of Foldy-Lax equations:

$$ \mathbf{E}_i = \mathbf{E}_{\text{inc}}(\mathbf{r}_i) + \alpha_s K^2 \sum_{j \neq i} G_0^{ij} \mathbf{E}_j. \quad (66) $$

The polarizability of a spherical particle is given by the Clausius-Mossotti formula, together with a radiative correction that accounts for the finite size of the scatterer³⁴:

$$ \alpha_s = 4\pi a^3 \beta_s \left[ 1 + \frac{2}{3} i \beta_s (Ka)^3 \right]. \quad (67) $$

Some other dynamical corrections have been proposed,³⁵,³⁶ but the difference between the corresponding imaginary parts is negligible for small size parameters. Taking the far-field expression of the Green's tensor leads to the following expression for the scattered amplitude:

$$ \tilde{\mathbf{E}}_s(\mathbf{K}) = \frac{\alpha_s K^2}{4\pi} \sum_{j=1}^{N} (1 - \hat{\mathbf{K}}\hat{\mathbf{K}}) \mathbf{E}_j e^{-i\mathbf{K}\mathbf{r}_j}. \quad (68) $$

### C. Validation Procedure

The validation of a homogenization procedure relies in general on the comparison of the different cross sections of the heterogeneous medium and a homogeneous one with the same shape. Usually, the comparison is made between the extinction of the random medium and the equivalent homogeneous sphere,¹⁰,¹⁴,¹⁹,³⁷,³⁸ as depicted in Fig. 3. We will proceed in the same way but in addition we will compute the absorption cross section. The extinction, scattering, and absorption cross section³⁹ of a scatterer at incident wave vector $\mathbf{K}_0$ and incident polarization $\mathbf{E}_0$ are respectively defined by

$$ \sigma_e = \frac{4\pi}{K} \Im\{[\tilde{\mathbf{E}}_s(\mathbf{K}) \cdot \mathbf{E}_0]\}_{\mathbf{K}=\mathbf{K}_0}, \quad (69) $$

$$ \sigma_s = \int_{4\pi} |\tilde{\mathbf{E}}_s|^2 d\Omega, \quad (70) $$

$$ \sigma_a = \sigma_e - \sigma_s. \quad (71) $$

Now the field scattered by a random aggregate can be decomposed into an average (coherent part) and a fluctuating value (incoherent part):

$$ \tilde{\mathbf{E}}_s = \langle \tilde{\mathbf{E}}_s \rangle + \delta \tilde{\mathbf{E}}_s. \quad (72) $$

For nonabsorbing inclusions we can equate the extinction and scattering cross section:

$$ \int_{4\pi} |\langle \tilde{\mathbf{E}}_s \rangle + \delta \tilde{\mathbf{E}}_s|^2 d\Omega = \frac{4\pi}{K} \Im[(\langle \tilde{\mathbf{E}}_s \rangle + \delta \tilde{\mathbf{E}}_s) \cdot \mathbf{E}_0]_{\mathbf{K}=\mathbf{K}_0}. \quad (73) $$

Taking the ensemble average of the last equation leads to the identity

$$ \int_{4\pi} |\langle \tilde{\mathbf{E}}_s \rangle|^2 d\Omega + \int_{4\pi} |\langle \delta \tilde{\mathbf{E}}_s \rangle|^2 d\Omega = \frac{4\pi}{K} \Im(\langle \tilde{\mathbf{E}}_s \rangle \cdot \mathbf{E}_0)_{\mathbf{K}=\mathbf{K}_0}. \quad (74) $$

Since the coherent scattered field of the random medium is identified with the field scattered from the homogeneous equivalent volume, an absorption term must be introduced to account for the incoherent part (which implies a complex effective permittivity). Hence the validation of the homogenization procedure relies on the identifications (with obvious notations):

$$ \sigma_s[V, \epsilon_e] = \sigma_s^{\text{coh}}[\text{Aggregate}], \quad (75) $$

$$ \sigma_a[V, \epsilon_e] = \sigma_s^{\text{incoh}}[\text{Aggregate}], \quad (76) $$

$$ \sigma_e[V, \epsilon_e] = (\sigma_s^{\text{coh}} + \sigma_s^{\text{incoh}})[\text{Aggregate}]. \quad (77) $$

To allow a rigorous computation of the field scattered by the homogenized medium, the embedding volume $V$ will be chosen as spherical. In that case, the different cross sections are given by the exact Mie formalism.

### D. Numerical Results

We will now proceed to a numerical validation of the identifications of Eqs. (75)–(77) for the homogenized medium based on formulas (55) and (56). A comparison will also be made with the homogenization according to the classical theories QCA and QCA-CP. For the hard-sphere homogeneous distribution, the second moment of the PCF can be evaluated analytically in the PY approximation³¹,³² by

$$ M_2 = \frac{1}{3f} \left[ \frac{(1-f)^4}{(1+2f)^2} - 1 \right], \quad (78) $$

which leads to a completely explicit formula for the effective permittivity of Eq. (52) predicted by QCA for lossless particles ($\beta_s''=0$). In the QCA-CP approximation, however, an implicit equation must be solved numerically:

<figure><img src="image_4.png" /><figcaption>Fig. 3. Illustration of the validation procedure. An identification is made between the different cross sections of the random and the homogenized medium.</figcaption></figure>
$$ \epsilon_e = 1 + 3f \frac{\epsilon_e(\epsilon_s - 1)}{3\epsilon_e + (\epsilon_s - 1)(1-f)} \left[ 1 + \frac{2i}{3}(Ka)^3\beta_s \frac{\epsilon_e^{3/2}}{3\epsilon_e + (\epsilon_s - 1)(1-f)}(1 + 3fM_2) \right]. \quad (79) $$

The performance of QCA and QCA-CP has been investigated in detail through numerical<sup>17,19,20,29</sup> as well as experimental<sup>40</sup> trials. Emphasis has been put on the extinction rate, which describes the attenuation of waves propagating in the medium. For the heterogeneous medium, it is given by the incoherent scattering cross section normalized by the volume of the bounding medium and is computed through Monte Carlo averages. For an infinite homogeneous medium, it is given merely by two times the imaginary part of the effective wavenumber normalized by the incident wavenumber.

A comparison has been made between the extinction rate of both random and effective media with EFA, QCA, and QCA-CP for fixed-size parameter of the particles but varying density. While the homogenization formulas are equivalent at low density, only QCA and QCA-CP show satisfactory comparison with Monte Carlo results at higher densities. We will reproduce closely these numerical experiments for the finite-size formula, with the following differences, however. Since the homogenized medium is finite, we compute its true absorption and extinction cross section instead of extinction rate. To this end, we use embedding volumes with spherical shape, instead of cubical as in Refs. 17 and 29. This makes the different cross sections available through the Mie solution and allows for a simplified homogenization formula [Eq. (55)].

The nonoverlapping spherical particles are randomly distributed in the spherical test volume as described in Section 2. The number N of scatterers is fixed at 2000 for each realization and the radius L of the embedding volume is taken to match a given density, $L=a(N/f)^{1/3}$, so that larger densities correspond to smaller volumes. The Monte Carlo simulations have been done for three different particle permittivities ($\epsilon_s=2.25, 3.2, 16$) and the size parameter of the small particles has been set to $Ka=0.1$. We have computed the total (i.e., coherent+incoherent) and incoherent scattering cross section for the corresponding random media at increasing filling ratio (from 10% to 40%). To this end, an average has been performed over 250 realizations. As explained earlier, this is to be compared with the extinction and absorption cross sections of a homogeneous sphere with permittivity given by formulas (52) (QCA), (79) (QCA-CP), and (55) (FS-QCA). Results are displayed in Figs. 4 and 5 for the extinction and Figs. 6 and 7 for the absorption. We have not shown the results for $\epsilon=2.25$, which are very similar to those of $\epsilon=3.2$. For the field computations, the incident wave vector and polarization have been taken along principal directions ($\hat{K}_0=\hat{z}$ and $\hat{E}_0=\hat{y}$). However, because of the perfect isotropy of the aggregate, there is no polarization dependence.

The extinction is predicted correctly by both QCA and FS-QCA, which are indistinguishable. QCA-CP, however, overestimates by far the extinction at higher densities, an

effect that increases dramatically with the dielectric constant. For the absorption, QCA and FS-QCA are both satisfactory at low densities and moderate permittivity ($\epsilon \le 3.2$), and very close to one another. This is to be expected since low densities here correspond to large vol-

<figure><img src="image_5.png" /><figcaption>Fig. 4. Extinction cross section of the homogeneous medium and average total scattered field of the aggregate for $\epsilon_s=3.2$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Same as Fig. 4 with $\epsilon_s=16$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Absorption cross section of the homogeneous medium and incoherent scattering cross section of the aggregate for $\epsilon_s=3.2$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Same as Fig. 6 with $\epsilon_s=16$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Imaginary part of the finite-size effective permittivity (FS-QCA) as a function of the radius of the embedding sphere, for otherwise fixed geometry ($K\alpha=0.1$, $\epsilon_s=3.2$, $f=0.3$). QCA is approached in the limit of large volumes.</figcaption></figure>

umes, in the limit of which both methods coincide. For higher densities, however, a clear separation between these last two methods appears, to the advantage of FS-QCA, which follows very closely the Monte Carlo values. QCA-CP is surprisingly disappointing at low densities, while it approaches FS-QCA at higher densities (hence smaller volumes). For higher permittivity, none of the methods is able to give a satisfactory absorption, and QCA-CP even shows a “pathological” behavior. The deterioration of FS-QCA and QCA is due to a more significant contribution of orders higher than two in the Born series, which is not taken into account by these methods. We have indeed checked that a Monte Carlo computation based on the first iteration of the Foldy–Lax equation instead of a rigorous inversion makes a perfect agreement with FS-QCA.

It is interesting to see how the predicted value of FS-QCA evolves with the size of the test volume. Figure 8 shows the imaginary part of $\epsilon_e$ [after Eq. (55)] as a function of the medium size, for fixed parameters $f=0.3$, $\epsilon_s=3.2$, and $K\alpha=0.1$, together with the imaginary part predicted by both QCA [Eq. (52)] and the complex MG [Eq. (41)]. FS-QCA makes a continuous transition from MG to

QCA as the size $L$ of the medium increases. The finite-size effect is marked until roughly 30 particle diameters, with about 15% difference at 20 diameters.

## 7. CONCLUSION

We have derived a homogenization formula for random aggregates that takes into account the finite size of the embedding volume, and approaches the classical QCA formula in the limit of large volumes. When the volume is smaller or comparable to the wavelength, numerical experiments have shown significant improvements over the classical theories (QCA and QCA-CP) for the estimation of the absorption. A finite-size formula for compound spheres has also been derived using the same technique and remains to be validated.

## APPENDIX A: CALCULATION OF JOINT PROBABILITIES

The event of having two positions $\mathbf{r}, \mathbf{r}'$ lying simultaneously in a particle can be decomposed into the two disjoint events according to whether these positions belong to the same particle or to distinct particles. Labeling the different particles from 1 to $N$, we may therefore write

$$ \begin{aligned} & \text{Proba}[\mathbf{r} \in \text{particle}, \mathbf{r}' \in \text{particle}] \\ &= N \text{Proba}[\mathbf{r}, \mathbf{r}' \in 1] + N(N-1)\text{Proba}[\mathbf{r} \in 1, \mathbf{r}' \in 2]. \end{aligned} \tag{A1} $$

Now, introducing the random variable $\mathbf{r}_j$, the center of particle $j$, we have for the first term

$$ \begin{aligned} & N \text{Proba}[\mathbf{r}, \mathbf{r}' \in 1] \\ &= \Pi_V(\mathbf{r})\Pi_V(\mathbf{r}')N \text{Proba}[|\mathbf{r} - \mathbf{r}'| < a, |\mathbf{r}' - \mathbf{r}'| < a] \\ &= \Pi_V(\mathbf{r})\Pi_V(\mathbf{r}')N \int d\mathbf{u}_1\rho_1(\mathbf{u}_1)\Pi_a(\mathbf{r} - \mathbf{u}_1)\Pi_a(\mathbf{r}' - \mathbf{u}_1) \\ &= f\Pi_V(\mathbf{r})\Pi_V(\mathbf{r}')\Phi\left(\frac{|\mathbf{r} - \mathbf{r}'|}{2a}\right), \end{aligned} $$

where we have used $f=NV_a/V$ and denoted $V_a=4/3\pi a^3$ the volume of one particle and $\Pi_a$ the characteristic function of such a centered particle [$\Pi_a(\mathbf{r})=1$ if $r<a$, $\Pi_a(\mathbf{r})=0$ otherwise]. For the second term we have

$$ \begin{aligned} \text{Proba}[\mathbf{r} \in 1, \mathbf{r}' \in 2] &= \text{Proba}[|\mathbf{r} - \mathbf{r}'| < a, |\mathbf{r}' - \mathbf{r}'| < a] \\ &= \Pi_V(\mathbf{r})\Pi_V(\mathbf{r}') \int d\mathbf{u}_1 d\mathbf{u}_2 \rho_2(\mathbf{u}_1, \mathbf{u}_2) \Pi_a(\mathbf{r} - \mathbf{u}_1) \Pi_a(\mathbf{r}') \\ &\quad - \mathbf{u}_2) = \Pi_V(\mathbf{r})\Pi_V(\mathbf{r}') \int d\mathbf{u}_1 d\mathbf{u}_2 \rho_1(\mathbf{u}_1) \rho_1(\mathbf{u}_2) g_2(\mathbf{u}_1) \\ &\quad - \mathbf{u}_2) \Pi_a(\mathbf{r} - \mathbf{u}_1) \Pi_a(\mathbf{r}' - \mathbf{u}_2). \end{aligned} $$
To proceed further we make the approximation $\Pi_a(\mathbf{r}-\mathbf{u}) = V_a \delta(\mathbf{r}-\mathbf{u})$ which, together with the approximation $N(N-1) \approx N^2$, leads to the expression

$$N(N-1)\text{Proba}[\mathbf{r} \in 1, \mathbf{r}' \in 2] = f^2\Pi_V(\mathbf{r})\Pi_V(\mathbf{r}')g_2(\mathbf{r}-\mathbf{r}'), \quad (A2)$$

and terminates the calculation of the joint probability.

## ACKNOWLEDGMENTS

The present work was supported by a grant of the ONERA, Toulouse. Many thanks also go to Matthias Holschneider for helpful discussions. Corresponding author C.-A. Guerin's e-mail address is charles-antoine.guerin@fresnel.fr.

## REFERENCES

1. J. Maxwell Garnett, “Colours in metal glasses and in metallic films,” Philos. Trans. R. Soc. London **203**, 385 (1904).

2. A. Sihvola, *Electromagnetic Mixing Formulas and Applications*, IEE Electromagnetic Waves Series (The Institute of Electrical Engineers, 1999).

3. L. Foldy, “The multiple scattering of waves,” Phys. Rev. **107**, 107–119 (1945).

4. M. Lax, “Multiple scattering of waves,” Rev. Mod. Phys. **23**, 287–310 (1951).

5. M. Lax, “Multiple scattering of waves. II. The effective field in dense systems,” Phys. Rev. **85**, 621–629 (1952).

6. V. Varadan, V. Bringi, and V. Varadan, “Coherent electromagnetic wave propagation though randomly distributed dielectric scatterers,” Phys. Rev. D **19**, 2480–2489 (1979).

7. L. Tsang and J. A. Kong, “Multiple scattering of electromagnetic waves by random distributions of discrete scatterers with coherent potential and quantum mechanical formalism,” J. Appl. Phys. **51**, 3465–3485 (1980).

8. L. Tsang and J. A. Kong, *Scattering of Electromagnetic Waves: Advanced Topics* (Wiley, 2001).

9. U. Frisch, “Wave propagation in random media,” in *Probabilistic Methods in Applied Mathematics*, A. T. Barucha-Reid, ed. (Academic, 1953), Vol. 1, pp. 75–198.

10. A. Doicu and T. Wriedt, “Equivalent refractive index of a sphere with multiple spherical inclusions,” J. Opt. A, Pure Appl. Opt. **3**, 204–209 (2001).

11. K. Fuller, “Scattering and absorption cross sections of compounded spheres III. Spheres containing arbitrarily located spherical inhomogeneities,” J. Opt. Soc. Am. A **12**, 893–904 (1995).

12. G. Videen, D. Ngo, P. Chýlek, and R. Pinnick, “Light scattering from a sphere with an irregular inclusion,” J. Opt. Soc. Am. A **12**, 922–928 (1995).

13. G. Videen, P. Pellegrino, D. Ngo, J. Videen, and R. Pinnick, “Light-scattering intensity fluctuations in microdroplets containing inclusions,” Atmos.-Ocean. **36**, 6115–6118 (1997).

14. P. Chýlek and G. Videen, “Scattering by a composite sphere and effective medium approximations,” Opt. Commun. **146**, 15–20 (1998).

15. P. Chýlek, G. Videen, and D. Ngo, “Effect of air bubbles on absorption of solar radiation by water droplets,” J. Atmos. Sci. **55**, 340–343 (1998).

16. K. H. Ding, L. M. Zurk, and L. Tsang, “Pair distribution functions and attenuation rates for sticky particles in dense media,” J. Electromagn. Waves Appl. **8**, 1585–1604 (1994).

17. L. M. Zurk, L. Tsang, K. H. Ding, and D. P. Winebrenner,

“Monte Carlo simulations of the extinction rate of densely packed spheres with clustered and nonclustered geometries,” J. Opt. Soc. Am. A **12**, 1772–1781 (1995).

18. L. Tsang, K. H. Ding, S. Shih, and J. Kong, “Scattering of electromagnetic wave from dense distributions of spheroidal particles based on Monte Carlo simulations,” J. Opt. Soc. Am. A **15**, 2660–2669 (1998).

19. B. E. Barrowes, C. O. Ao, F. L. Teixeira, J. A. Kong, and L. Tsang, “Monte Carlo simulation of electromagnetic waves propagation in dense random media with dielectric spheroids,” IEICE Trans. Electron. **E83-C**, 1797–1802 (2000).

20. C. Ao and J. Kong, “Analytical approximations in multiple scattering of electromagnetic waves by aligned dielectric spheroids,” J. Opt. Soc. Am. A **19**, 1145–1156 (2002).

21. A. Yaghjian, “Electric dyadic Green’s function in the source region,” Proc. IEEE **68**, 248–263 (1980).

22. G. Smith, “Dielectric constant for mixed media,” J. Phys. D **10**, L39–42 (1977).

23. D. Stroud and F. Pan, “Self-consistent approach to electromagnetic wave propagation in composite media: application to model granular theories,” Phys. Rev. B **17**, 1602 (1978).

24. G. Niklasson, C. Granqvist, and O. Hunderi, “Effective medium models for the optical properties of inhomogeneous materials,” Appl. Opt. **20**, 26–30 (1981).

25. P. Chýlek and V. Srivastava, “Dielectric constant of a composite inhomogeneous medium,” Phys. Rev. B **27**, 5098 (1983).

26. P. Mallet, C. A. Guérin, and A. Sentenac, “Maxwell–Garnett mixing rule in the presence of multiple scattering: derivation and accuracy,” Phys. Rev. B **72**, 014205 (2005).

27. C. Tai, *Dyadic Green's Functions in Electromagnetic theory* (Intext, 1971).

28. N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, “Equation of state calculation by fast computing machines,” J. Chem. Phys. **21**, 1087–1092 (1953).

29. L. Tsang, J. A. Kong, K. Ding, and C. O. Ao, *Scattering of Electromagnetic Waves: Numerical Simulations* (Wiley, 2001).

30. J. Percus and G. Yevick, “Analysis of classical statistical mechanics by means of collective coordinates,” Phys. Rev. **110**, 1–13 (1958).

31. M. Wertheim, “Exact solution of the Percus-Yevick integral equation for hard spheres,” Phys. Rev. Lett. **10**, 321–323 (1963).

32. R. Baxter, “Method of solution of the Percus-Yevick, hypernetted-chain or similar equation,” Phys. Rev. **154**, 170–175 (1967).

33. Y. Xu, “Electromagnetic scattering by an aggregate of spheres,” Appl. Opt. **34**, 4573–4588 (1995).

34. B. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. **333**, 848–872 (1988).

35. W. T. Doyle, “Optical properties of a suspension of metal spheres,” Phys. Rev. B **39**, 9852–9858 (1989).

36. A. Lakhtakia, “Size-dependent Maxwell-Garnett formula from an integral equation formalism,” Optik (Stuttgart) **91**, 134–137 (1992).

37. G. Videen and P. Chýlek, “Scattering by a composite sphere with an absorbing inclusion and effective medium approximations,” Opt. Commun. **158**, 1–6 (1998).

38. L. Kolokolova and B. A. Gustafson, “Scattering by inhomogeneous particles: microwave analog experiments and comparison to effective medium theories,” J. Quant. Spectrosc. Radiat. Transf. **70**, 611–625 (2001).

39. C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles* (Wiley, 1983).

40. R. West, D. Gibbs, L. Tsang, and A. Fung, “Comparison of optical scattering experiments and the quasi-crystalline approximation for dense media,” J. Opt. Soc. Am. A **11**, 1854–1857 (1994).