SENTENAC A.: Effective-medium theory for finite-size aggregates

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.


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, 1 many expressions for the effective permittivity of a composite medium have been proposed.The simplest ones are the so-called mixing rules. 2 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), [3][4][5][6][7][8] 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 9 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.
1][12][13][14] The main application is atmospheric optics, where water droplets and aerosols may include microcontaminants that modify their scattering and absorption properties. 15In 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.][18][19][20] 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-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.

DESCRIPTION OF THE RANDOM MEDIUM
We consider an assembly of N nonoverlapping identical spherical particles of radius a and permittivity ⑀ s located at random positions 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 where The n-point distribution functions are usually written in the form n ͑r 1 , . . .,r n ͒ = 1 ͑r 1 ͒ . . . 1 ͑r n ͒g n ͑r 1 , . . .,r n ͒, ͑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 ജ 3 and one contents oneself with the knowledge of the socalled pair-correlation function (PCF), g 2 .The PCF has two essential asymptotic properties: which accounts for the absence of overlap between small spherical particles, and which expresses the mutual decorrelation of particles at large separation.For a spatially homogeneous distribution, the PCF depends only on the separation of particles:

CALCULATION OF THE SCATTERED FIELD A. Extended Born Series
We now assume that the aggregate is illuminated by an incident monochromatic plane wave Here and everywhere an implicit time dependence exp͑−it͒ is assumed.We denote K 0 the incident wave vector and E 0 the (constant) polarization vector.The electrical properties of the random medium are described by a random permittivity function ⑀͑r͒ = ⑀ s if r belongs to a particle and ⑀͑r͒ = 1 otherwise.The electric field satisfies the well-known integral equation 8 where K =2 / is the wavenumber in vacuum and G 0 is the free-space dyadic Green's function (I is the unit dyad): The expanded expression of G 0 away from its singularity r = 0 is given by The Green's tensor can be decomposed into a singular part and a principal value (PV): where L is a constant dyad that depends on the shape of the exclusion domain chosen to define the PV. 21For a spherical exclusion domain L =−1/3I, and we may thus rewrite where we have used the notation ͛ for the principal value integral with spherical exclusion domain around the singularity: Plugging the iteration series of Eq. ( 11) into Eq.( 7) we obtain the following expression for the scattered field E s = E − E inc for r outside the test volume: where we have introduced a rephased Green's tensor and a source term ͓␤͔ that is written as a functional expansion in ␤, where ␤ is the dimensionless function The successive terms of this series are found to be

͑17͒
The series (15) is expected to converge for small values of 3K 2 ␤, that is, for small contrast, the geometry held fixed.

B. Coherent Scattered Field
Denote E s ͑K͒ the far-field amplitude in the remote direction K ˆ: The ensemble average ͗E ˜s͑K͒͘ is the so-called coherent field.Using the far-field expression of the Green's tensor, we obtain easily where the dyad K ˆK ˆis the projector along the direction K ˆ.
Hence the coherent field is simply related to the ensemble average ͓͗␤͔͑r 1 ͒͘.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 ␤͑r͒.Therefore one usually 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 ␤͑r͒, respectively: where B 2 ͑r,rЈ͒ ª ͗␤͑r͒␤͑rЈ͒͘.͑22͒ As we will see, the contribution of the second-order term ͗ 2 ͓␤͔͘ is necessary to predict an imaginary part of the effective permittivity that is not accounted for by the firstorder term.The statistical mean B 1 ͑r͒ is easily calculated, since it is merely proportional to the probability of having a particle at position r: where f is the volume density of the scatterers, The correlation B 2 ͑r , rЈ͒ is given by the probability that the positions r and rЈ lie simultaneously in a particle: B 2 ͑r,rЈ͒ = ␤ s 2 Proba͓r particle,rЈ particle͔.͑26͒ From the calculations in Appendix A we obtain where

͑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: where r 12 = r 1 − r 2 and 1+2 = 1 + 2 .The second term on the right-hand side can be simplified by using the isotropy of ⌽.Performing a preliminary angular integration we may rewrite

͑30͒
where ␥ is related to the average of the Green's tensor over a sphere (S r is the sphere of radius r): The expression of ␥ can be easily derived by using the integration formulas and dS r e −iK•r r ˆr ˆ= 4r 2 ͭͫ sin͑Kr͒ and noting that K ˆ0K ˆ0E 0 = 0. Explicitly, we obtain

͑33͒
Now for small particles ͑Ka Ӷ 1͒ we may use ⌽ Ӎ 1 in the integral (30) and retain the Taylor expansion of ␥ about zero: Altogether this leads to the simplified expression, valid for Ka Ӷ 1:

͑35͒
where ␤ s Ј and ␤ s Љ are the real and imaginary parts of ␤ s , and discarding the O͑␤ s 3 ͒ terms, we may rewrite this in the simple final expression where ‫ء‬ stands for the convolution operator.

HOMOGENIZATION PROCEDURE A. Identification
We will now seek an effective volume with index ␤ e whose scattered field coincides with the coherent scattered field of the random medium; that is, ͓␤ e ͔ = ͓͗␤͔͘.3][24][25] The price to pay for an overall identification is the occurrence of a nonconstant permittivity.Finding 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 ͑ 1 ͓␤ e ͔ = ͗ 1 ͓␤͔͒͘ leads to the equality ␤ e = f␤ s , or 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 ͑ 1+2 ͓␤ e ͔ = ͗ 1+2 ͓␤͔͒͘ is expressed by a more complex equation: Note that this identity is evidently satisfied with g 2 = 1 and ␤ e = ␤ 0 or, in terms of effective permittivity,

͑41͒
It is interesting to note that this is again the effective constant ⑀ MG predicted by the classical MG mixing rule, apart from the imaginary term.This last term can, however, 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 8 for lossless particles ͑␤ s Љ=0͒ which coincides analytically to within the O͑␤ s 2 ͒ term: We will now seek the actual index ␤ e as a multiplicative correction to ␤ 0 inside the volume V: ␤ e ͑r͒ = ͓1 + ͑r͔͒␤ 0 .

͑43͒
Inserting Eq. (43) into Eq.( 40) leads to the following relation for the perturbation : where we have set

͑45͒
Finally the effective permittivity is given by inside the test volume V, with 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.

Large Volumes
Although the Born series is not expected to converge at large volumes, we may take the formal limit V → ϱ in Eq. ( 47), thereby obtaining a constant value For isotropic PCF, the last integral can be turned into the radial integral

Short-Range Correlations
Equation (47) can also be made explicit in the case of short-range and isotropic correlations.Denote by 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, ϳ a, and if furthermore the particles are small ͑Ka Ӷ 1͒, then we may the retain the Taylor expansion (34) of ␥͑u͒ about zero and replace the upper bound by the infinity in integral.After a simple change of variables, we obtain where and g ˜2 is the normalized PCF defined for particles with unit radius, that is g 2 ͑r͒ = g ˜2͑r / a͒.This leads to the following expression for the effective permittivity, to within O͓͑Ka͒ 3 ͔ terms:

͑52͒
which coincides analytically with the low-frequency QCA approximation 8 for lossless particles ͑␤ s

Spatially Averaged Index
In the general case where the correlation is not shortrange 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:

͑53͒
For a spherical test volume with diameter L, this is easily found to be where ⌽ is the intersection volume of normalized spheres [Eq.( 28)].For isotropic PCF again, the last expression can be simplified to

͑55͒
which is consistent with the large-volume formula (49).The corresponding effective permittivity, will be tested in Section 6.

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 ͑⑀ h ͒ of radius L, immersed in vacuum and containing random spherical inclusions (⑀ s , radius a).The distribution of small particles is identical to that described in Section 2 that is extracted from a spatially homogeneous, hardsphere distribution.The dielectric properties of the aggregate are now characterized by the random function ⑀͑r͒ which can take three different values ͑⑀ s , ⑀ h ,1͒ according to whether r is inside one particle, outside the particle but in the test volume V, or outside the test volume, respectively.We denote by ⑀ Mie ͑r͒ the permittivity function in the absence of inclusions (homogeneous sphere ⑀ h ) and by E Mie the reference solution corresponding to the total field produced by this homogeneous sphere when illuminated by incident field E inc .We denote by G Mie the dyadic Green's function of the sphere, whose expression can be found in Tai. 27The electric field integral equation now reads

͑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: 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: where G h is given by Eq. ( 8) with K replaced by the wavenumber K h in the medium ͑K h 2 = ⑀ h K 2 ͒.This second approximation is justified since G 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 Mie ͑r , rЈ͒ for close points ͑r , rЈ͒ is imposed by G h ͑r − rЈ͒.
Following the same identification procedure as in Section 4, we recover successively the MG mixing rule in the host medium, where the EFA in the host medium, and a finite-size QCA where G h = e −iK 0 •r E 0 • ͑G h E 0 ͒ and ␤ 0h is defined by Eq. ( 37) with the replacements ␤ s → ␤ sh and K → K h .A more tractable expression for Eq. ( 64) can again be obtained by considering the spatially averaged index ␤ ¯e and specializing to isotropic PCF g 2 .In that case formula ( 55) is recovered with the replacements K → K h and ␤ 0 → ␤ 0h .

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 twopoint 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. 28 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, 30 for which an exact solution is known. 31,32A 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 33 or the method of moments 18 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 ͱ ⑀Ka, where ⑀ and a are the permittivity and radius of the particle, respectively.In this case the field scattered by the aggregate is given merely by where E j is the exciting field on the jth dipole and ␣ s its polarizability.The exciting fields are mutually related by the system of Foldy-Lax equations: 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 34 : Some other dynamical corrections have been proposed, 35,36 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:

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, 10,14,19,37,38 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 39 of a scatterer at incident wave vector K 0 and incident polarization E 0 are respectively defined by Now the field scattered by a random aggregate can be decomposed into an average (coherent part) and a fluctuating value (incoherent part):

͑72͒
For nonabsorbing inclusions we can equate the extinction and scattering cross section: Taking the ensemble average of the last equation leads to the identity

͑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): 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 homogenous distribution, the second moment of the PCF can be evaluated analytically in the PY approximation 31,32 by which leads to a completely explicit formula for the effective permittivity of Eq. ( 52) predicted by QCA for lossless particles ͑␤ s Љ=0͒.In the QCA-CP approximation, however, an implicit equation must be solved numerically:  The performance of QCA and QCA-CP has been investigated in detail through numerical 17,19,20,29 as well as experimental 40 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 ͑⑀ 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 ⑀ = 2.25, which are very similar to those of ⑀ = 3.2.For the field computations, the incident wave vector and polarization have been taken along principal directions (K ˆ0 = z ˆand E ˆ0 = 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 ͑⑀ ഛ 3.2͒, and very close to one another.This is to be expected since low densities here correspond to large vol-Fig.4. Extinction cross section of the homogeneous medium and average total scattered field of the aggregate for ⑀ s = 3.2.Fig. 5. Same as Fig. 4 with ⑀ s = 16.Fig. 6.Absorption cross section of the homogeneous medium and incoherent scattering cross section of the aggregate for ⑀ s = 3.2.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 ⑀ e [after Eq. ( 55)] as a function of the medium size, for fixed parameters f = 0.3, ⑀ s = 3.2, and Ka = 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.

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 r , 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 Proba͓r particle,rЈ particle͔ = N Proba͓r,rЈ 1͔ + N͑N − 1͒Proba͓r 1,rЈ 2͔.

͑A2͒
and terminates the calculation of the joint probability.

Fig. 1 .
Fig. 1.Random aggregate of particles is obtained by cutting a finite volume inside an infinite distribution.

Fig. 2 .
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).

Fig. 3 .
Fig. 3. Illustration of the validation procedure.An identification is made between the different cross sections of the random and the homogenized medium.