# Maxwell-Garnett mixing rule in the presence of multiple scattering: Derivation and accuracy

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

► To cite this version:

Pierre Mallet, Charles-Antoine Guérin, Anne Sentenac. Maxwell-Garnett mixing rule in the presence of multiple scattering: Derivation and accuracy. Physical Review B: Condensed Matter and Materials Physics (1998-2015), 2005, 72 (1), pp.014205. 10.1103/PhysRevB.72.014205. hal-00080462

HAL Id: hal-00080462

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

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
# Maxwell-Garnett mixing rule in the presence of multiple scattering: Derivation and accuracy

Pierre Mallet,<sup>1,2</sup> C. A. Guérin,<sup>1</sup> and Anne Sentenac<sup>1</sup>

<sup>1</sup>Institut Fresnel and Université Paul Cézanne, Faculté de Saint-Jérôme, F-13397 Marseille cedex 20, France

<sup>2</sup>ONERA Centre de Toulouse, BP 4025, 31055 Toulouse Cedex 4, France

We give a rigorous and original derivation of the Maxwell-Garnett mixing rule in the dynamical regime for a composite dielectric random medium with small spherical inclusions. For certain configurations of scatterers, we show that contrarily to the common belief, the Maxwell-Garnett formula can remain very accurate at a high concentration of scatterers and incorporate multiple-scattering effects as well as attenuation of the mean field. We provide a realistic numerical example for which this is the case.

## I. INTRODUCTION

Effective medium theory (EMT) is a powerful way to handle the radiative properties of composite materials such as, e.g., cosmic dusts, aerosols, porous media, and recently metamaterials. One of the oldest and most popular EMT is the Maxwell-Garnett (MG) mixing rule¹ (see Ref. 2 for a classical derivation and Refs. 3 and 4 for a review). MG theory was originally derived by invoking the dipolar character of scatterers and neglecting their density fluctuations about a mean value. A modern heuristic approach consists in requiring that the forward-scattering amplitude of a cell of the composite medium, immersed inside the effective homogeneous background, be zero.⁵⁻⁸ These approaches can incorporate finite-size effects of scatterers. Their main limitation is that they do not discriminate between two random media with same density of scatterers but different statistical distributions. Neglecting the *n*-particle distribution functions is justified in the single-scattering regime, where the mean field inside the material actually only depends on the particle density. This is perhaps the reason why it is often claimed that a weak particle interaction is a condition for the validity of MG theory.

More sophisticated models based on multiple-scattering theory have been proposed (see Refs. 9–12 for a few milestones, Ref. 13 for an exhaustive review) to take into account the higher-order statistics as well as the finite-size effects of the scatterers. Tractable expressions for the effective permittivity, however, can only be obtained by retaining the sole pair-correlation function and taking the long-wavelength regime. For most of these expressions, the MG formula is found to be the limiting case as the size of inclusions goes to zero (in the quasistatic framework, it has even been shown analytically¹⁴,¹⁵ that the MG formula can incorporate multipolar effects within a mean-field approximation). However, we are not aware of any direct derivation of the MG formula in the dynamical regime that takes properly into account the *n*-particle distribution and sets explicitly the limits of this approximation. A clarification is needed all the more as there is a persistent belief in the literature that the MG formula only remains valid for dilute systems of Rayleigh scatterers.

Yet the accuracy of the MG mixing rule has been largely investigated and tested against other EMT's, from a numerical¹⁶⁻²² as well as experimental²³⁻²⁶ point of view.

While the MG model is in general observed to deteriorate at high filling ratio of inclusions, it has been shown for certain configurations<sup>22,23,27</sup> to remain surprisingly accurate at high particle density, in the presence of strong mutual interactions. The purpose of this paper is to remove this apparent contradiction by providing a rigorous and at the same time simple justification of the MG model in the dynamical case for Rayleigh scatterers in the so-called cermet topology—that is, for separated inclusions in the host medium. We show that in certain statistical configurations, the radiative properties of a random medium can be accurately described by the MG effective permittivity, even for a high concentration of scatterers, where strong multiple scattering is present. Furthermore, we show that this effective permittivity can account correctly for the attenuation of the mean field in the material, without resorting to more advanced effective medium theories.<sup>13</sup>

Our approach consists in comparing the coherent field radiated by a random aggregate of particles confined in a given test volume *V* with the field scattered by an homogeneous object, with same shape and same volume *V*, and effective permittivity *ε*<sub>e</sub>. The method is based on the matching of the iteration series for the coherent field of the aggregate and the Born series of the field scattered by an equivalent homogeneous medium, a procedure which is to our best knowledge original. We present Monte Carlo simulations for configurations that correspond to random defects in periodic composite media and investigate the role of multiple scattering and the influence of the statistical distribution. In that case we show that the MG mixing rule remains accurate at very high density, contrarily to the alternative classical mixing rule: namely, the Bruggeman law.

## II. RIGOROUS APPROACH TO THE MAXWELL-GARNETT FORMULA

### A. Description of the random medium

We consider an aggregate of *N* nonoverlapping dielectric identical particles, located at random positions *r*<sub>j</sub> within a volume *V* and illuminated by a monochromatic radiation of wave number *K*=2π/λ. The particles have relative permittivity *ε*<sub>s</sub> and are small with respect to the wavelength λ, so that they can be assumed to radiate like dipoles. For the sake of clarity we will make two simplifying assumptions, which
can be relaxed in the demonstration (see the discussion at the end of the section). The host medium will be a vacuum and the particles will be chosen to be spherical. The polarizability of a sphere of radius a is given by the Clausius-Mossotti (CM) formula, with radiative correction,²⁸

$$ \alpha_s = \alpha_0 \left( 1 + \frac{2i}{3} \beta_s (Ka)^3 \right), \qquad (1) $$

where

$$ \alpha_0 = 4\pi a^3 \frac{\epsilon_s - 1}{\epsilon_s + 2} = 4\pi a^3 \beta_s. \qquad (2) $$

Some other dynamical corrections to the CM formula have been proposed.²⁹,³⁰ For small-size parameter $Ka$, the difference between the corresponding imaginary parts is negligible. The influence of this finite-size correction will, however, be discussed subsequently.

## B. Mean field for the aggregate

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

$$ E^0(\mathbf{r}) = e^{iK_0 R} E^0. \qquad (3) $$

Here $K_0$ is the incident wave vector and $E^0$ is a (constant) polarization vector. We assume that the particles are actually dipoles, with a polarizability (1) reminiscent of their finite size. The impact of this simplification will be discussed in

Sec. III D. Denote $E_j$ the exciting field on the jth dipole and $E_j^0=E^0(r_j)$ the incident field. The exciting fields for point scatterers are mutually related by the Foldy-Lax equations⁹,³¹

$$ E_j = E_j^0 + \alpha_s K^2 \sum_{i=1, i \neq j}^{N} G_{ij} E_j, \quad j = 1, \dots, N, \qquad (4) $$

where $G_{ij}=G_0(r_{ij})$ is the free-space dyadic Green's function evaluated at the separation $r_{ij}=r_i-r_j$. For $r \neq 0$, $G_0(r)$ is given by

$$ \left[ \left( 1 + \frac{i}{Kr} - \frac{1}{(Kr)^2} \right) \hat{1} + \left( \frac{3}{(Kr)^2} - \frac{3i}{Kr} - 1 \right) \hat{r}\hat{r} \right] \frac{e^{iKr}}{4\pi r}, \qquad (5) $$

where $\hat{1}$ is the unit dyad and $\hat{r}\hat{r}$ is the projector on the direction $\hat{r}$. We recall that $G_0(r)p$ is the electric field produced at distance $r$ by a dipole of moment $p$ (see, e.g., Ref. 32, p. 411) placed at the origin.

The scattered field $E_s$ is given by

$$ E_s(\mathbf{r}) = \alpha K^2 \sum_{j=1}^{N} G(\mathbf{r} - \mathbf{r}_j) E_j. \qquad (6) $$

The system (4) can be solved by iteration, provided the corresponding series converges. The scattered field outside the spherical volume can then be rewritten:

$$ E_s(\mathbf{r}) = E_s^{(1)}(\mathbf{r}) + E_s^{(2)}(\mathbf{r}) + E_s^{(3)}(\mathbf{r}) + \dots, \qquad (7) $$

where

$$ E_s^{(n)}(\mathbf{r}) = (\alpha K^2)^n \sum_{i_1} \sum_{i_2 \neq i_1} \cdots \sum_{i_n \neq i_{n-1}} G(\mathbf{r} - \mathbf{r}_{i_1}) \cdots G(\mathbf{r}_{i_{n-1}} - \mathbf{r}_{i_n}) E^0(\mathbf{r}_{i_n}) \qquad (8) $$

is the contribution of the nth iteration to the scattered field. The first term ($n=1$) is the single-scattering approximation. The ensemble average $\langle E_s(\mathbf{r}) \rangle$ is the so-called coherent scattered field. Averaging Eq. (7) we obtain the corresponding iteration series, with generic term:

$$ \langle E_s^n(\mathbf{r}) \rangle = (\alpha N K^2)^n \int d\mathbf{r}_1 G(\mathbf{r} - \mathbf{r}_1) \cdots \int_{|\mathbf{r}_{n-1} - \mathbf{r}_n| > 2a} d\mathbf{r}_n \rho_n(\mathbf{r}_1, \dots, \mathbf{r}_n) G(\mathbf{r}_{n-1}, \dots, \mathbf{r}_n) E^0(\mathbf{r}_n), \qquad (9) $$

where $\rho_n$ denotes the n-point probability distribution function (PDF) of particle positions. In this averaging procedure we have used $N-1 \approx N$ and have neglected the occurrence of repeated entries in the n-uple $(\mathbf{r}_{i_1}, \dots, \mathbf{r}_{i_n})$, so that the configuration space is actually described by the n-point PDF. This is justified when $N \gg n$—that is, for a large number of particles.

## C. Effective homogeneous medium

We now consider an homogeneous medium of complex relative permittivity $\epsilon_e$ with the same volume $V$, surrounded by a vacuum and illuminated by the same incident field. From the combined Maxwell harmonic curl equations, the electric field is seen to satisfy the well-known curl-curl equation

$$ \mathrm{curlcurl} E(\mathbf{r}) - \epsilon(\mathbf{r}) K^2 E(\mathbf{r}) = i\omega\mu_0\mathbf{J}, \qquad (10) $$

where $\omega$ is the harmonic time pulsation, $J$ is the source of current that produces the incident field, and $\epsilon(\mathbf{r})=\epsilon_e$ if $\mathbf{r} \in V$ and $\epsilon(\mathbf{r})=1$ elsewhere. The previous equation may be conveniently rewritten as a free space-equation with additional source terms:

$$ \mathrm{curlcurl} E(\mathbf{r}) - K^2 E(\mathbf{r}) = i\omega\mu_0\mathbf{J} + K^2(\epsilon_e - 1)\Pi(\mathbf{r})E(\mathbf{r}), \qquad (11) $$

where $\Pi$ is the characteristic function of the volume $[\Pi(\mathbf{r})=1 \text{ if } \mathbf{r} \in V, \Pi(\mathbf{r})=0 \text{ elsewhere}]$. The unique solution that
satisfies the outgoing wave condition is then given by an integral equation involving the free-space dyadic Green's function:

$$E(\mathbf{r}) = E^0(\mathbf{r}) + K^2(\varepsilon_e - 1) \int_V d\mathbf{r}' \Pi(\mathbf{r}') G(\mathbf{r} - \mathbf{r}') E(\mathbf{r}'). \quad (12)$$

Due to its singularity on the diagonal, the Green's kernel is decomposed into a singular part and a principal value $\mathcal{P}$ (Ref. 33):

$$G(\mathbf{r}-\mathbf{r}') = \frac{1}{K^2} \mathcal{L} \delta(\mathbf{r}-\mathbf{r}') + \mathcal{P} G(\mathbf{r}-\mathbf{r}'), \quad (13)$$

where $\mathcal{L}$ is a constant dyad which depends on the shape of the exclusion domain chosen to define the principal value and $\mathcal{P}G$ coincides with expression (5) for $\mathbf{r}-\mathbf{r}' \neq 0$. For a spherical exclusion domain $\mathcal{L} = -\frac{1}{3}\hat{1}$, where $\hat{1}$ is the unit dyad, and we may thus rewrite inside the volume $V$,

$$E(\mathbf{r}') = \frac{3}{\varepsilon_e + 2} E^0(\mathbf{r}') + 3K^2 \beta_e \oint d\mathbf{r}'' \Pi(\mathbf{r}'') G(\mathbf{r}' - \mathbf{r}'') E(\mathbf{r}''), \quad (14)$$

where the dimensionless constant $\beta_e$ is defined by

$$\beta_e = \frac{\varepsilon_e - 1}{\varepsilon_e + 2} \quad (15)$$

and where we have used the notation $\oint$ for the principal value integral with spherical exclusion domain around the singularity:

$$\int d\mathbf{r}' \mathcal{P} G(\mathbf{r}-\mathbf{r}') = \oint d\mathbf{r}' G(\mathbf{r}-\mathbf{r}') \\ := \lim_{a \to 0} \int_{|\mathbf{r}-\mathbf{r}'| > 2a} d\mathbf{r}' G(\mathbf{r}-\mathbf{r}'). \quad (16)$$

Plugging the iteration series of Eqs. (14) into (12) we obtain the Born series for the far-field scattered by the volume $V$:

$$E_s^{(n)}(\mathbf{r}) = (3K^2\beta_e)^n \int d\mathbf{r}_1 \Pi_L(\mathbf{r}_1) G(\mathbf{r}-\mathbf{r}_1) \\ \times \oint \dots \oint d\mathbf{r}_n \Pi_L(\mathbf{r}_n) G(\mathbf{r}_{n-1,n}) E^0(\mathbf{r}_n). \quad (17)$$

The first term in the series ($n=1$) is the so-called (distorted) Born approximation.

### D. Identification

Comparing the generic terms in Eqs. (9) and (17) and keeping in mind the definition (16) we see that there is a possible identification at arbitrary order $n$ between the two series in the limit of small spheres ($a \to 0$) if

$$\rho_n(\mathbf{r}_1, \dots, \mathbf{r}_n) = \frac{1}{V^n} \Pi(\mathbf{r}_1) \cdots \Pi(\mathbf{r}_n) \quad (18)$$

and

$$\alpha_s N = 3\beta_e V. \quad (19)$$

Introducing the particle density

$$f = \frac{4\pi N a^3}{3V} \quad (20)$$

and making use of Eqs. (1) and (2), the last condition can be reformulated as

$$\beta_e = f \beta_s \left( 1 + \frac{2}{3} i (K a)^3 \beta_s \right). \quad (21)$$

The first condition is realized whenever the particle positions are uniform and independent of the test volume $V$. The second relation is precisely the MG formula. In terms of effective permittivity it reads

$$\epsilon_e = 1 + 3f \frac{\beta_s}{1 - \beta_s f} \left( 1 + \frac{2i}{3} (K a)^3 \frac{\beta_s}{1 - \beta_s f} \right). \quad (22)$$

Note that the usual MG expression does not contain the imaginary term depending on $(K a)^3$, for it is derived in the static framework, where no radiative correction applies [$\alpha_s = \alpha_0$ in Eq. (1)]. However, we will still refer to this formula as the “MG formula,” as it is the same equation (19) in terms of polarizability.

### E. Discussion

The above derivation shows that the MG mixing rule gives in principle the exact expression of the coherent field at all orders in the Born series, provided the particles are of dipolar character ($K a \ll 1$) and have uncorrelated positions ($\rho_n = \text{const}$). These last two assumptions are at the root of a frequent confusion in the literature.

(i) The independence of particles positions does not imply the absence of an electromagnetic interaction or multiple scattering. Therefore the vague statement of “independent particles,” which is often called on as a basic assumption for MG mixing rule, should refer to spatial correlation only and not to the coupling between scatterers.

(ii) The particles radiating like dipoles do not imply a quasistatic regime. In particular, the test volume or the interparticle distance need not be small with respect to the wavelength. In the subsequent numerical computations we provide examples where there is a pronounced dynamical effect.

The previous analysis has shown that any deviation from the MG formula for the effective permittivity of an aggregate of Rayleigh scatterers can only be due to a nonuniform $n$-particle PDF. Note, however, that the single-scattering term in the iteration series depends only on the one-particle PDF, which is merely the filling ratio, and can be matched exactly with the Born approximation for the homogeneous sphere. Hence, the MG mixing rule for an homogeneous ensemble of small scatterers is always accurate in the absence of multiple scattering, regardless of the correlations between
particle positions. Another consequence of the derivation is the occurrence of a complex value for the effective permittivity (22), even for nonabsorbing particles ($\epsilon_s$ real). The classical MG formula is usually presented as the real part of this expression, and the introduction of an imaginary correction taking into account the finite-size of scatterer is in general presented as an extended theory. Hence the basic MG mixing rule can also predict some attenuation inside the materials which is due to the diffusion process itself.

At this point it is time to comment on some simplifying assumptions that have been used so long. First, the identification that has been made between the coherent scattered field of the aggregate and the field scattered by a homogeneous sphere relies on the crucial assumption that the Born series converge. However, the Born series happen to converge for a small test volume, regardless of the density, as is shown in the Appendix. This volume need only be statistically relevant—that is, much larger than the typical inclusion size and interparticle distance. Second, we have taken the embedding medium to be a vacuum. In the case of a nonvoid host matrix, the free-space Green's function has to be replaced by the volume Green's function. The term-by-term identification in the series still holds, even though the Green's function is not explicit. Finally, the particle can be assumed nonspherical. This amounts to changing $\alpha_s$ to some other polarizability and taking an exclusion domain of the same shape³³ as the particle in the extraction of the Green's function singularity [Eq. (16)].

### III. VALIDATION PROCEDURE

To validate an homogenization procedure one generally compares the different cross sections of the random media with that of the homogenized medium.

The extinction, scattering, and absorption cross sections (e.g., Ref. 34) of a scatterer with incident wave number $K_0$ and polarization $E_0$ are, respectively, defined by

$$ \sigma_e = \frac{4\pi}{K} \mathrm{Im}\{(\tilde{E}_s(K) \cdot E_0)\}_{K=K_0}, \quad (23) $$

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

$$ \sigma_a = \sigma_e - \sigma_s, \quad (25) $$

where $\tilde{E}_s(K)$ is the far-field scattering amplitude in the remote direction $K[E_s=\tilde{E}_s(K)e^{iKR}/R]$. In our example, the homogenized volume will be chosen to be spherical, so that the different cross sections can rigorously and easily by obtained with the Mie solution (see any classical textbook—e.g. Ref. 34). In addition to this rigorous solution, we will also systematically compute the Born solution, which is the single-scattering approximation for the homogeneous medium.

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

<figure><img src="image_2.png" /><figcaption>FIG. 1. Principle of the homogenization in a spherical test volume. We compare the coherent and incoherent scattering cross sections of the random medium with the scattering, extinction, and absorption cross sections of the homogenized volume.</figcaption></figure>

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

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

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

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

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

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 validation of the homogenization procedure relies on the identification (with obvious notations)

$$ \sigma_s[\text{homogeneous}, \epsilon_e] = \sigma_s^{\text{coh}}[\text{random}], \quad (29) $$

$$ \sigma_a[\text{homogeneous}, \epsilon_e] = \sigma_s^{\text{incoh}}[\text{random}], \quad (30) $$

$$ \sigma_e[\text{homogeneous}, \epsilon_e] = (\sigma_s^{\text{coh}} + \sigma_s^{\text{incoh}})[\text{random}]. \quad (31) $$

The principle of the homogenization procedure is depicted in Fig. 1.

## IV. NUMERICAL RESULTS

### A. Random medium generation

The previous analysis has shown that the MG mixing rule is in principle exact if the particle positions are perfectly uncorrelated—that is, if all the n-point probability distribution functions are constant ($\rho_n$=const for all n). Such a distribution of particles cannot exclude overlap, and therefore the condition under which the MG mixing rule is mathematically true can never be realized for hard spheres. However, it
<figure><img src="image_3.png" /><figcaption>FIG. 2. Two random composite media with same inclusion density but different correlation structure. In the top figure (medium 1), the probabilities of occupation of two neighbor sites on the lattice are independent. In the bottom figure (medium 2), the probabilities of occupation of neighbor sites are strongly correlated.</figcaption></figure>

is possible to design a composite system where the particle positions are prescribed but independent from one another, making them almost uncorrelated. This can be realized by placing the particles at the nodes of a regular lattice (within a spherical test volume) with a random probability of occupation ($p$). The elementary cells are cubes of side $2a$, so that contact is allowed between the particles and the volume density is given by $f=\pi/6p$. This model can represent a perturbed periodic structure such as a photonic crystals with defaults. This random medium will hereafter be referred to as medium 1. To evidence the dramatic impact of strong spatial correlations on the validity of the MG mixing rule, we have generated a random medium (medium 2) with identical density but strong interdependence between neighboring sites. This random assembly is generated in the following way: an empty initial site is chosen at random on the lattice, and a random walk is performed from there in the cardinal directions. A particle is placed on each visited site. The random walk is stopped and restarted elsewhere whenever an already occupied site is reached. It can be verified that the obtained density of scatterers is spatially homogeneous. Figure 2 shows a realization of the different media at equal density.

## B. Field computations

The scattering and exciting fields in the aggregate are obtained by solution of the Foldy-Lax equations (4), which can

be solved by iteration. The lattice structure of the aggregate gives the interaction matrix $G_{ij}$ a Töplitz structure (that is, depending on $i-j$ only) in each space direction. This allows a fast computation of the successive iterates by use of a three-dimensional fast Fourier transform (FFT) (see, e.g., Ref. 35). The computational effort is independent of the particle density and depends only on the size of the embedding medium. The system is solved for each realization of a random medium. The coherent and incoherent scattering cross sections are then obtained by a Monte Carlo average over a large number of realizations. The scattering and absorption cross sections of the corresponding homogeneous medium with MG index are obtained from the Mie solution. The incident wave vector $K_0$ will be chosen along a principal direction (z), and the incident polarization $E_0$ will be chosen horizontal (along x) or vertical (along y).

## C. Numerical experiments

The spherical test volume has been chosen comparable to the incident wavelength to disqualify any possible quasistatic approximation. A numerical compromise was found between the computational time and the statistical significance of this embedding volume which was set to a diameter $2L=64a$. The Monte Carlo simulation has been performed over 600 realizations. The particle density ranges from low to high filling ratios $0<f<0.42$ ($f=0.52$ is the full lattice), and the permittivity has been set to three typical values, from weak to strong: $\epsilon_s=2.25, 3.2$, and 16. The size parameter has been set to $Ka=0.1$, to exclude eventual size effects of the inclu-sions. For all these parameters, the iterative solutions have been checked to converge. We will now detail three numeri-cal experiments that exemplify the main claims of this paper.

### 1. First experiment

We have computed the coherent scattered field for medium 1 (MC1) at high filling ratio ($f=0.41$) and strong inclusion contrast ($\epsilon=3.2$). Figure 3 shows the differential scattering cross section $d\sigma/d\Omega$ for both the aggregate (MC1) and the homogeneous sphere (Mie), in the incidence plane $(x,z)$ and vertical polarization ($E_0=E_0y$). The (polar) scattering angle ranges from 0 (specular) to 100. The angular range [100, 180] has been discarded to make the figure more readable, and the negative angles ([-180, 0]) are completed by symmetry. The single-scattering solutions have been given to show the importance of multiple scattering. This plot shows two important features: (i) The agreement between the random (MC1) and homogenized media (Mie) is excellent over the whole range of scattering angles. (ii) There is a strong multiple scattering effect.

The conclusion is identical in horizontal polarization. We have also checked that there is a strong difference with the quasistatic solution (which is obtained by taking the static Green's function $[K^2G \rightarrow (3\hat{r}\hat{r}-1)/(4\pi r^3)]$ in Foldy-Lax equations).

### 2. Second experiment

To check the influence of the particle correlation function we have computed the coherent scattered field for medium 1
<figure><img src="image_4.png" /></figure>

FIG. 3. First experiment. Differential coherent scattering cross section for random medium 1 (MC1, thick line) and comparison with the differential cross section of the homogenized medium with the same volume (Mie solution) and MG effective permittivity (squares). There is good agreement between the random and homogenized media. The single-scattering solution (MC1, thin line) and the Born approximation for the homogeneous sphere are also given to show the importance of multiple scattering.

(MC1) and medium 2 (MC2) at low density ($f=0.05$), but strong inclusion contrast ($\epsilon=16$). The former condition ensures a strongly nonuniform two-point PDF for medium 2, while the latter ensures the occurrence of multiple scattering. Since the density of both random media is the same, the effective permittivity given by the MG mixing rule is the same. The corresponding differential cross section is shown in Fig. 4. From this plot we can draw the following conclusions: (i) The prediction of the MG mixing rule is very good

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

FIG. 4. Second experiment. Differential coherent scattering cross section for the weakly (MC1) and strongly (MC2) correlated random media with same density and comparison with the differential cross section of the homogenized test volume with MG effective permittivity (Mie). Only the weakly correlated medium MC1 is correctly homogenized. The Born approximation for the homogeneous medium with MG effective permittivity is also given. It coincides with the single-scattering approximation of both MC1 and MC2 (which are not shown here) and is therefore unable to discriminate these different media.

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

FIG. 5. Third experiment. Incoherent scattering cross section at increasing density for the weakly (MC1) and strongly (MC2) correlated random media and comparison with the absorption cross section of the homogeneous sphere with MG effective permittivity.

for MC1, but not for MC2. (ii) There is a pronounced multiple-scattering effect.

We have also checked that MC1 and MC2 have the same single-scattering cross section. This should be the case as soon as both media have the same density, for the single-scattering approximation depends on this last quantity only. This shows once more that the accuracy of the MG mixing rule is determined by the *n*-point PDF's.

### 3. Third experiment

We have computed the total extinction and incoherent scattering cross sections of medium 1 and 2 with $\epsilon=3.2$ for increasing densities and compared with the extinction and absorption cross sections of the corresponding homogeneous sphere (Fig. 5 for the absorption). The extinction of the aggregate is computed from the forward-scattering amplitude via the optical theorem [Eq. (23)]. The fluctuations due to the averaging procedure are within 0.02%.

(i) The predicted extinction is extremely accurate for MC1 (actually within 3% of the relative error). This accuracy was to be expected from the first experiment since scattering is the dominant part in the extinction.

(ii) The predicted extinction is less accurate for MC2 (about a 4% error), but still very satisfying, even at high filling ratios, contrarily to what was observed in the second experiment. This is due to a weaker dielectric constant for the particles ($\epsilon=3.2$), leading to a smaller contribution of multiple scattering. This shows that the total extinction cross section is very robust to the introduction of correlation between particles.

(iii) The agreement of the absorption curve with the MG prediction is acceptable for MC1 until a filling ratio of 20%, where the relative error is about 20%. The absorption obtained for MC2, however, is 10 times the order of magnitude of MG. This shows that absorption is extremely sensitive to the correlation functions of the random aggregate (much more than extinction).
<figure><img src="image_7.png" /><figcaption>FIG. 6. Influence of the size of scatterers on the extinction cross section. A comparison is given between the classical theories (MG and Bruggeman) and the extended theories (EMG and EB), which bring a size correction to the effective permittivity. The configuration is that of the first experiment.</figcaption></figure>

## D. Discussion

(i) Together with the MG theory, the oldest and most popular mixing rule is the Bruggeman mixing rule$^{36}$

$$f\frac{\epsilon_s - \epsilon_{Brugg}}{\epsilon_s + 2\epsilon_{Brugg}} + (1-f)\frac{1 - \epsilon_{Brugg}}{1 + 2\epsilon_{Brugg}} = 0, \quad (32)$$

which is widely believed to be superior to the former at high filling ratio. This seems contradictory with the excellent results obtained with MG on the test case MC1 at high particle density. We have therefore also computed the Mie extinction cross section of the homogenized sphere with Bruggeman effective constant for the configuration of the first experiment ($\epsilon_s=3.2, f=0.41$); see Fig. 6. The use of the Bruggeman mixing rule results in a significant deterioration of the homogenization procedure.

(ii) The next question which may arise from our numerical experiments is whether the finite-size effect of scatterers has been properly taken into account. We have chosen to include a radiative correction for the particle polarization [Eq. (1)] based on the definition of Draine.$^{28}$ However, several finite-size corrections have been proposed in the literature for particle polarization and, consequently, the MG formula. One might wonder whether the previous results are sensitive to the type of correction which is chosen. There are many such extended theories, and we will not cite or test them all (we refer to, e.g., Ref. 37 for a review). One classical approach is the one by Lakhtakia$^{30}$ using an integral equation formalism. Assuming the electric field to be constant within each particle, he obtains a modified polarizability for small spheres and, correspondingly, the so-called extended Maxwell-Garnett (EMG) mixing rule (we express here the relative permittivity with respect to vacuum)

$$\epsilon_e = \frac{3 + 2fb(1, \epsilon_s; Ka)}{3 - fb(1, \epsilon_s; Ka)}, \quad (33)$$

where

<figure><img src="image_8.png" /><figcaption>FIG. 7. Same as Fig. 6 for the absorption.</figcaption></figure>

$$b(\epsilon_2, \epsilon_1; Ka) = \frac{\epsilon_1 - \epsilon_2}{1 + (1 - \epsilon_1/\epsilon_2) \left[ \frac{2}{3} (1 - iKa\sqrt{\epsilon_2}e^{iKa\sqrt{\epsilon_2}} - 1) \right]} \quad (34)$$

Reproducing Bruggeman's approach with the same finite-size correction b for the polarizability of small spheres, Lakhtakia$^{38}$ has also proposed a so-called extended Bruggeman (EB) approach, where the effective permittivity is given by

$$fb(\epsilon_{EB}, \epsilon_s; Ka) + (1-f)b(\epsilon_{EB}, 1; Ka) = 0. \quad (35)$$

It can be easily checked that the difference between $\epsilon_{MG}$ and $\epsilon_{EMG}$ is $O((Ka)^2)$ for the real part and $O((Ka)^4)$ for the imaginary part. For the size parameter we have chosen ($Ka=0.1$), the difference with the classical MG formula based on Eq. (1) is invisible for the extinction and negligible for the absorption. A comparison is shown in Figs. 6 and 7. The same holds for the Bruggeman approach and its extended version.

(iii) There are a certain number of advanced effective medium theories based on multiple-scattering theory (see Ref. 13 and references therein), which involve the distribution functions of particles and yield an improved estimate of the imaginary part of the effective permittivity. The most employed are the so-called quasi crystalline approximation (QCA) and the quasi crystalline approximation with coherent potential (QCA-CP), which are obtained via renormalization techniques and consist in keeping terms involving double-scattering processes only in the diagrammatic series. In their low-frequency form, these approximate theories provide an explicit analytical form of the effective constant in term of the second moment of the pair correlation function $g(r)$. However, this last function is in general unknown or difficult to estimate. Therefore these theories are essentially used in the context of a hard-sphere distribution, for which the Perkus-Yevick approximation provides a precise and simple expression for the moments of the pair-correlation function. It is interesting to note that the QCA differs from the MG formula by an additive term in the imaginary part:
$$ \epsilon_{QCA} - \epsilon_{MG} \propto f \int dr \, r^2 [g(r) - 1]. \quad (36) $$

This shows that the MG formula is as accurate as the QCA in the ideal case g=1 and/or at low density f, and that the MG formula deteriorates together with an increase of the integral on the right-hand side (assuming the QCA is the “reference” solution, even though it involves essentially double scattering).

## V. CONCLUSION

The objective of this paper was to point out a simple and rigorous way to derive the Maxwell-Garnett expression in a way that accounts for the density fluctuations of the particles and for the presence of multiple scattering. Our aim was to set up clearly the limits of validity of this expression by stressing the importance of the particle correlations and the role of multiple scattering. We have shown that, if single scattering is the dominant mechanism, the coherent and incoherent scattering cross sections of the finite random media are accurately described by the scattering and absorption cross sections of the homogenized test volume with MG permittivity. This is true whatever the statistical distribution of the particles as long as it is uniform. On the other hand, if multiple scattering is important, the effective permittivity depends strongly on the n-point probability distribution function of the aggregate. Yet, if the particle positions can be considered uncorrelated, the effective permittivity reduces to the MG expression. Thus, in this particular case, the MG mixing rule remains valid in a regime where it is usually not expected to hold—that is, for a high density of scatterers with strong electromagnetic interaction.

## ACKNOWLEDGMENTS

The present work was partially supported by a grant of the O.N.E.R.A. Toulouse. Many thanks also go to Hervé Tortel for helping us with the Mie computations.

## APPENDIX: THE CONVERGENCE OF THE BORN SERIES

The Foldy-Lax equations (4) can be rewritten in a matrix form

$$ \bar{E} = \bar{E}_0 + \bar{A}\bar{E}, \quad (A1) $$

with obvious notations. We define the norm of the 3N-dimensional vector $\bar{E}$ by

$$ \|\bar{E}\| = \sup_{j=1,...,N} \|E_j\|_{\infty}, \quad (A2) $$

where $\|\cdot\|_{\infty}$ is the supremum norm over the coordinates and the operator norm of $\bar{A}$ by

$$ \|\bar{A}\| = \sup_{\|E\|=1} \|\bar{A}\bar{E}\|. \quad (A3) $$

Using the fact that the particles cannot interpenetrate the following rough estimation can be easily obtained:

$$ \|\bar{A}\| \le \frac{N|\alpha|K^3}{2\pi a} \left( 1 + \frac{3}{(2Ka)^2} + \frac{3}{(2Ka)^3} \right). \quad (A4) $$

This upper bound is proportional to the volume and density of scatterers and thus can be made arbitrary small (in particular smaller than one) by decreasing the former or the latter, with other parameters held fixed.

<footnote n="1">¹J. C. Maxwell-Garnett, Philos. Trans. R. Soc. London, Ser. A **203**, 385 (1904).</footnote>

<footnote n="2">²M. Born and E. Wolf, *Principles of Optics*, 3rd ed. (Pergamon, Oxford, 1965).</footnote>

<footnote n="3">³A. Sihvola, J. Electromagn. Waves Appl. **15**, 715 (2001).</footnote>

<footnote n="4">⁴Ari Sihvola, *Electromagnetic Mixing Formulas and Applications* (The Institution of Electrical Engineers, New York, 1999).</footnote>

<footnote n="5">⁵G. B. Smith, J. Phys. D **10**, L39 (1977).</footnote>

<footnote n="6">⁶D. Stroud and F. P. Pan, Phys. Rev. B **17**, 1602 (1978).</footnote>

<footnote n="7">⁷G. A. Niklasson, C. G. Granqvist, and O. Hunderi, Appl. Opt. **20**, 26 (1981).</footnote>

<footnote n="8">⁸P. Chylek and V. Srivastava, Phys. Rev. B **27**, 5098 (1983).</footnote>

<footnote n="9">⁹M. Lax, Phys. Rev. **85**, 621 (1952).</footnote>

<footnote n="10">¹⁰W. Lamb, D. M. Wood, and N. W. Ashcroft, Phys. Rev. B **21**, 2248 (1980).</footnote>

<footnote n="11">¹¹L. Tsang and J. A. Kong, J. Appl. Phys. **51**, 3465 (1980).</footnote>

<footnote n="12">¹²V. A. Davis and L. Schwartz, Phys. Rev. B **31**, 5155 (1985).</footnote>

<footnote n="13">¹³L. Tsang and J. A. Kong, *Scattering of Electromagnetic Waves, Advanced Topics* (Wiley, New York, 2001).</footnote>

<footnote n="14">¹⁴F. Claro and R. Rojas, Phys. Rev. B **43**, 6369 (1991).</footnote>

<footnote n="15">¹⁵Liang Fu, Pedro B. Macedo, and Lorenzo Resca, Phys. Rev. B **47**, 13818 (1993).</footnote>

<footnote n="16">¹⁶Ruben G. Barrera, Pedro Villasenor-Gonzalez, W. Luis Mochan, and Guillermo Monsivais, Phys. Rev. B **41**, 7370 (1990).</footnote>

<footnote n="17">¹⁷J. M. Perrin and P. L. Lamy, Astrophys. J. **364**, 146 (1990).</footnote>

<footnote n="18">¹⁸W. C. Chew, J. A. Friedrich, and R. Geiger, IEEE Trans. Geosci. Remote Sens. **GE-28**, 207 (1990).</footnote>

<footnote n="19">¹⁹Ruibao Tao, Zhe Chen, and Ping Sheng, Phys. Rev. B **41**, 2417 (1990).</footnote>

<footnote n="20">²⁰S. Datta, C. T. Chan, K. M. Ho, and C. M. Soukoulis, Phys. Rev. B **48**, 14936 (1993).</footnote>

<footnote n="21">²¹Gorden Videen and Petr Chýlek, Opt. Commun. **158**, 1 (1998).</footnote>

<footnote n="22">²²Anna Spanoudaki and Rolf Pelster, Phys. Rev. B **64**, 064205 (2001).</footnote>

<footnote n="23">²³B. Abeles and J. I. Gittleman, Appl. Opt. **15**, 2328 (1976).</footnote>

<footnote n="24">²⁴Sung-Ik Lee, Won Noth, Kevin Cummings, and J. R. Gaines, Phys. Rev. Lett. **55**, 1626 (1985).</footnote>

<footnote n="25">²⁵G. Koh, IEEE Trans. Geosci. Remote Sens. **GE-30**, 184 (1992).</footnote>

<footnote n="26">²⁶Ludmilla Kolokolova and Bo A. D. Gustafson, J. Quant. Spectrosc. Radiat. Transf. **70**, 611 (2001).</footnote>

<footnote n="27">²⁷Feng Wu and Keith W. Whites, IEEE Trans. Antennas Propag. **49**, 1174 (2001).</footnote>

<footnote n="28">²⁸B. T. Draine, Astrophys. J. **333**, 848 (1988).</footnote>

<footnote n="29">²⁹Wiliam T. Doyle, Phys. Rev. B **39**, 9852 (1989).</footnote>
³⁰A. Lakhtakia, Optik (Stuttgart) **91**, 134 (1992).

³¹L. L. Foldy, Phys. Rev. **67**, 107 (1945).

³²J. D. Jackson, *Classical Electrodynamics*, 3rd ed. (Wiley, New York, 1999).

³³A. D. Yaghjian, Proc. IEEE **68**, 248 (1981).

³⁴Craig F. Bohren and Donald R. Huffman, *Absorption and Scattering of Light by Small Particles* (Wiley, New York, 1983).

³⁵J. J. Goodman, B. T. Draine, and P. J. Flatau, Opt. Lett. **16**, 1198 (1991).

³⁶D. A. G. Bruggeman, Ann. Phys. **24**, 636 (1935).

³⁷R. Ruppin, Opt. Commun. **182**, 273 (2000).

³⁸A. Lakhtakia, Optik (Stuttgart) **96**, 25 (1994).