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

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 1 ͑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.6][7][8] 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 14,15 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.
While the MG model is in general observed to deteriorate at high filling ratio of inclusions, it has been shown for certain configurations 22,23,27 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. 13ur 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 ⑀ e .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.

A. Description of the random medium
We consider an aggregate of N nonoverlapping dielectric identical particles, located at random positions r j within a volume V and illuminated by a monochromatic radiation of wave number K =2 / .The particles have relative permittivity ⑀ s 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, 28 Some other dynamical corrections to the CM formula have been proposed. 29,30For 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: 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 9,31 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 0, G 0 ͑r͒ is given by where 1 ˆis the unit dyad and r ˆr ˆis the projector on the direction 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 The system ͑4͒ can be solved by iteration, provided the corresponding series converges.The scattered field outside the spherical volume can then be rewritten: where is the contribution of the nth iteration to the scattered field.The first term ͑n =1͒ is the single-scattering approximation.
The ensemble average ͗E s ͑r͒͘ is the so-called coherent scattered field.Averaging Eq. ͑7͒ we obtain the corresponding iteration series, with generic term: where n denotes the n-point probability distribution function ͑PDF͒ of particle positions.In this averaging procedure we have used N −1Ӎ N and have neglected the occurrence of repeated entries in the n-uple ͑r i 1 , … , r i n ͒, so that the configuration space is actually described by the n-point PDF.This is justified when N ӷ n-that is, for a large number of particles.

C. Effective homogeneous medium
We now consider an homogeneous medium of complex relative permittivity ⑀ 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 curlcurlE͑r͒ − ⑀͑r͒K 2 E͑r͒ = i 0 J, ͑10͒ where is the harmonic time pulsation, J is the source of current that produces the incident field, and ⑀͑r͒ = ⑀ e if r V and ⑀͑r͒ = 1 elsewhere.The previous equation may be conveniently rewritten as a free space-equation with additional source terms: where ⌸ is the characteristic function of the volume ͓⌸͑r͒ =1 if r V , ⌸͑r͒ = 0 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:

͑12͒
Due to its singularity on the diagonal, the Green's kernel is decomposed into a singular part and a principal value P ͑Ref.33͒: where L is a constant dyad which depends on the shape of the exclusion domain chosen to define the principal value and PG coincides with expression ͑5͒ for r − rЈ 0. For a spherical exclusion domain L =− 1 3 1 ˆ, where 1 ˆis the unit dyad, and we may thus rewrite inside the volume V,

͑14͒
where the dimensionless constant ␤ e is defined by and where we have used the notation ͛ for the principal value integral with spherical exclusion domain around the singularity: Plugging the iteration series of Eqs.͑14͒ into ͑12͒ we obtain the Born series for the far-field scattered by the volume V: The first term in the series ͑n =1͒ is the so-called ͑dis-torted͒ 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 → 0͒ if and Introducing the particle density and making use of Eqs.͑1͒ and ͑2͒, the last condition can be reformulated as

͑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 Note that the usual MG expression does not contain the imaginary term depending on ͑Ka͒ 3 , for it is derived in the static framework, where no radiative correction applies ͓␣ s = ␣ 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 ͑Ka Ӷ 1͒ and have uncorrelated positions ͑ n = 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 ͑⑀ 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 ␣ s to some other polarizability and taking an exclusion domain of the same shape 33 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 where E ˜s͑K͒ is the far-field scattering amplitude in the remote direction K͓E s = 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 singlescattering 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͒:

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

͑27͒
Taking the ensemble average of the last equation leads to the identity

͑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 no-tations͒ The principle of the homogenization procedure is depicted in Fig. 1.

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 ͑ 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 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 = /6 p.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: ⑀ s = 2.25, 3.2, and 16.The size parameter has been set to Ka = 0.1, to exclude eventual size effects of the inclusions.For all these parameters, the iterative solutions have been checked to converge.We will now detail three numerical experiments that exemplify the main claims of this paper.

First experiment
We have computed the coherent scattered field for medium 1 ͑MC1͒ at high filling ratio ͑f = 0.41͒ and strong inclusion contrast ͑⑀ = 3.2͒.Figure 3 shows the differential scattering cross section d / d⍀ for both the aggregate ͑MC1͒ and the homogeneous sphere ͑Mie͒, in the incidence plane ͑x , z͒ and vertical polarization ͑E 0 = E 0 y͒.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 2 G → ͑3r ˆr ˆ−1͒ / ͑4r 3 ͔͒ in Foldy-Lax equations).

Second experiment
To check the influence of the particle correlation function we have computed the coherent scattered field for medium 1 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.
͑MC1͒ and medium 2 ͑MC2͒ at low density ͑f = 0.05͒, but strong inclusion contrast ͑⑀ =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 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 singlescattering 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.

Third experiment
We have computed the total extinction and incoherent scattering cross sections of medium 1 and 2 with ⑀ = 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 ͑⑀ = 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͒.

D. Discussion
͑i͒ Together with the MG theory, the oldest and most popular mixing rule is the Bruggeman mixing rule 36 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 ͑⑀ 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. 28However, 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͒ Reproducing Bruggeman's approach with the same finitesize 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͑⑀ EB ,⑀ s ;Ka͒ + ͑1 − f͒b͑⑀ EB ,1;Ka͒ = 0. ͑35͒ It can be easily checked that the difference between ⑀ MG and ⑀ 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 doublescattering 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: ⑀ QCA − ⑀ MG ϰ f ͵ dr r 2 ͓g͑r͒ − 1͔.͑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 with obvious notations.We define the norm of the 3N-dimensional vector E ¯by Using the fact that the particles cannot interpenetrate the following rough estimation can be easily obtained: 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.

FIG. 1 .
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.

FIG. 3 .
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.

FIG. 6 .
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.