On applicability of differential mixing rules for statistically homogeneous and isotropic dispersions

The classical differential mixing rules are assumed to be independent effective-medium approaches, applicable to certain classes of systems. In the present work, the inconsistency of differential models for macroscopically homogeneous and isotropic systems is illustrated with a model for the effective permittivity of simple dielectric systems of impenetrable balls. The analysis is carried out in terms of the compact group approach reformulated in a way that allows one to analyze the role of different contributions to the permittivity distribution in the system. It is shown that the asymmetrical Bruggeman model (ABM) is physically inconsistent since the electromagnetic interaction between previously added constituents and those being added is replaced by the interaction of the latter with recursively formed effective medium. The overall changes in the effective permittivity due to addition of one constituent include the contributions from both constituents and depend on the system structure before the addition. Ignoring the contribution from one of the constituents, we obtain generalized versions of the original ABM mixing rules. They still remain applicable only in a certain concentration ranges, as is shown with the Hashin-Shtrikman bounds. The results obtained can be generalized to macroscopically homogeneous and isotropic systems with complex permittivities of constituents.


I. INTRODUCTION
A big variety of analytical methods and approaches have been developed to study electrophysical characteristics of disperse systems and mixtures [1][2][3][4][5]. However, it is often unclear which one is applicable to a given system [1,[6][7][8]. Even if a particular approach can be used for one type of systems, it can be unapplicable to other similar systems. Of the most used, but arguable approaches are differential models [9][10][11][12][13][14].
The first differential mixing rule was developed by Bruggeman for the effective permittivity of mixtures [9,10]; it is now known as the asymmetrical Bruggeman model (ABM). Later, the ABM was generalized by Hanai and other authors [11][12][13][14] to obtain the complex permittivity of mixtures with different types of inclusions (the so-called Bruggeman-Hanai or Maxwell-Wagner-Hanai model).
It should be noted that the original Hanai's generalization is based on the Maxwell-Wagner model [11,15] and thus incorporates the interfacial polarization effects (known as Maxwell-Wagner polarization). Generalizations of the differential approach to bi-anisotropic systems were developed by Lakhtakia and co-workers [16,17] (incremental and differential Maxwell-Garnett formalisms).
The term "asymmetrical" means that within the ABM, the system's constituents are divided into "the inclusions" (filler particles) and "the matrix" (host medium) 1 . * andrey.k.semenov@gmail.com 1 Obviously, this division is not rigorous: as the volume concentra-To derive the differential equation for the effective permittivity, the following procedure is envisaged. Given a particular system, whose effective characteristics are formed by the host and already added inclusions, suppose that an infinitesimal portion of inclusions is added to it. This addition causes infinitesimal changes in the effective characteristics, including the permittivity, of the system. Therefore, another new portion of inclusions will be added to a medium with a new effective permittivity. So, as the desired system is built by successive additions of infinitesimal portions of inclusions, each new portion is added to a medium with its own effective permittivity, different from that of the preceding medium. Correspondingly, the portions added at different steps contribute differently to the effective permittivity formation. The differences between these contributions can be neglected for diluted systems with low dielectric contrast, where the system's constituents interact weakly.
For other concentration ranges the ABM is, strictly speaking, inapplicable, but we are unaware of any rigorous analysis of this matter.
To get rid of the indicated limitations, each addition of a new portion of inclusions should be analyzed by taking into account the previously added inclusions [18]. This feature is intrinsic to the symmetrical Bruggeman model (SBM) [19], where all constituents of the system are treated alike, and no differential equation is involved.
Both SBM and ABM methods belong to the class of effective-medium approaches where one or all the constituents of the system are embedded into some effective medium. However, the suggested ways for modeling this medium are different in these methods. In the SBM, the effective medium is formed by all the constituents (including the real host). In differential models, the effective medium is formed recursively, according to the Maxwell-Garnett rule [20] (or Maxwell-Wagner rule [15] if the permittivities are complex-valued), by successive addition of small portions of inclusions to the current effective medium, starting from the pure matrix (for details of the recursive procedure see, for instance, [16,17]). Nevertheless, ignoring the above-indicated restriction on the applicability ranges, differential models are widely used for electric spectroscopy studies of water/oil emulsions [11,[21][22][23][24], soils [6], sands [25] and rocks [13,26], and biological samples [7,27], where the standard SBM does not work properly.
The goal of this research is to scrutinize the internal inconsistency and the ranges of validity of the classical differential models for systems of particles with real-(the ABM and its modifications) and complex-valued (the Maxwel-Wagner-Hanai model and its modifications) permittivities. In order to do this and avoid insignificant, within the scope of the research, specific effects (such as absorption, interfacial Maxwell-Wagner polarization), we develop a generalized differential approach to the effective quasistatic permittivity of macroscopically homogeneous and isotropic dielectric mixtures and then apply it to the simplest system of impenetrable ("hard") balls embedded in a uniform host medium. The results obtained are, in fact, general and applicable as long as the conditions of macroscopic homogeneity and isotropy are fulfilled.
The working model is based on the compact groups approach (CGA) [5,[28][29][30]. The "compact groups" are macroscopic regions in a system with typical linear sizes L that are much smaller than the wavelength of probing field λ in the system: L ≪ λ. Such groups are pointlike with respect to the field, but preserve the effective properties of the system. This fact allows one to effectively estimate the long-wavelength many-particle contributions to the average electric field and induction in macroscopically homogeneous and isotropic systems. In addition, based on the CGA, some relations, essential for electrodynamic homogenization and usually set ab initio [4,31,32], can be substantiated by using the boundary conditions for the electric field (see A) or the Hashin-Shtrikman variational principle [33] (see [34]). The analysis within the CGA actually reduces to simple modeling of the dielectric permittivity distribution in the system.
Very recently [34], the CGA was used to describe the effective dielectric response of dispersions of graded impenetrable and hard-core-penetrable-shell particles. The validity of the approach was demonstrated in [34] by contrasting its results with existing rigorous analytic results and computer simulations for dispersions of hard dielec-tric spheres with power-law permittivity profiles, and by processing experimental data on the effective dielectric response of nonconducting polymer-ceramic composites. Earlier, the CGA was efficiently applied to dispersions of particles with complex permittivities to describe electric percolation phenomena in composites of core-shell particles [5], two-step electrical percolation in nematic liquid crystals filled with multiwalled carbon nanotubes [35], and effective parameters of suspensions of nanosized insulating particles [36]. Finally, the idea of compact groups was also used by Sushko to evaluate the effects of multiple short-range reemissions between particles on the mean free path and the transport mean free path of photons in concentrated suspensions [37] and to discover the 1.5 molecular light scattering in fluids near the critical point [38,39]; the results were supported by extensive experimental data.
The above arguments show that the CGA is wellsubstantiated, flexible, and efficient. For these reasons, it was chosen as a basis for achieving the stated goals. The in-depth analysis of the CGA for the model under consideration can be found in [28][29][30]34].

II. GENERAL THEORETICAL BACKGROUND
The effective permittivity ε is determined as the proportionality coefficient between the average induction D and the average electric field E : where ǫ 0 is the electric constant, ε(r) = ε f + δε(r) is the local value of permittivity in the system, and the angular brackets stand for statistical averaging. For infinite systems, the latter is equivalent to volume averaging, according to the ergodic hypothesis [40]. The CGA allows one to find these averages in the long-wavelength limit. For macroscopically homogeneous and isotropic systems [28][29][30] The functional form of δε(r) is constructed according to the system under consideration. For a system of hard balls, with permittivity ε 1 , embedded in a host matrix, with permittivity ε 0 , ε(r) is equal to ε 1 in the regions occupied by the balls, and ε 0 otherwise. Correspondingly, δε(r) takes the form whereχ 1 (r) is the characteristic function of the region occupied by all the balls. The parameter ε f is found in A and in [34]; it equals the effective permittivity: ε f = ε.
In order to find ε, one should sum up the series in (2), (3). For (4), this gives the following equation for ε [29]: where c is a volume fraction of the inclusions phase. According to the CGA, the long-wavelength results (2), (3), and ε f = ε are valid for any macroscopically homogeneous and isotropic system. Thus, they can be applied to any distribution δε satisfying these conditions and are equivalent to the relation (see A or [34]) The functional form of δε(r) can vary significantly. However, it can always be categorized into two classes, depending on the way in which the system's constituents are treated: (1) symmetrical models (e.g. the SBM), where all the constituents are treated equally; (2) asymmetrical models (e.g. the ABM), where one phase is treated as "the inclusion" phase and the other as "the matrix" phase. Some simplest forms of δε and the corresponding mixing rules are as follows: 1) The SBM where δε formally coincides with (4) at ε f = ε. This means that each inclusion and the matrix are embedded in an effective medium with looked-for permittivity ε. The homogenization condition ε f = ε is actually the basic assumption of the model. The SBM-type distributions lead to the same result (5) for ε. However, it should be noted that the original SBM equation deals only with systems of hard spheres. The constituents' polarizabilities are identified with their individual polarizabilities in the effective medium, and the matrix is assumed to polarize as a single particle does [1]. These two suggestions, used also for systems of non-spherical particles, are contradictory [18].
In terms of the CGA, the form (4) is the most physically reasonable, since it represents the local permittivity ε(r) in the system. As was already mentioned, the latter is required to be macroscopically homogeneous and isotropic, while the actual form of the inclusions is insignificant. This fact significantly distinguishes the CGA from the original SBM.
2) The ABM can also be reproduced easily. Assume that the effective permittivity ε of the system at a certain amount of inclusions is known (see fig. 1(a)). If a small portion of inclusions with the characteristic function ∆χ 1 (r) (χ 1 · ∆χ 1 = 0) is added (the enclosed area in fig. 1(a)), the effective permittivity will change by ∆ε ( fig. 1(b)). The current effective medium with permittivity ε is considered as the matrix for the new inclusions that is void of any inclusions and determined by the characteristic function (1 −χ 1 − ∆χ 1 ). In terms of δε, this assumption can be written as Figure 1. Schematic representation of the ABM differential algorithm: (a) addition of a portion of new particles with concentration ∆c/(1−c) in the particle-void region to the current effective medium with permittivity ε (the lighter area) leads to (b) the formation of a new effective medium with permittivity ε + ∆ε, which serves as the matrix for the next portions of inclusions. Therefore, the previously added portions interact electrically with the new one only via the effective medium (comprising the particles shown darker).
where only the terms of the first orders of smallness in ∆χ 1 (in the meaning of its average) and ∆ε are retained; ε in (6) also should be changed to ε+∆ε, however it is not necessary for this particular form of δε. Substituting (7) into (6), taking into account the ergodic hypothesis, and using the orthogonality condition (1 −χ 1 − ∆χ 1 )∆χ 1 = 0 for the characteristic functions, averaging over the entire system in (6) can be split into averaging over the region occupied by the matrix and averaging over the region occupied by a new portion of inclusions: where again only the terms of the first order of smallness were retained. Changing to the infinitesimal values dε and dc gives the differential equation The point c = 1 is a critical one; the solution of (8) should satisfy there the condition ε = ε 1 . The ABM mixing rule is obtained by integrating the left side of (8) with respect to c from zero to c and the right side with respect to ε from ε 0 to ε: This equation has the same form as the one for the complex permittivity (the integration is performed then using Morera's and Cauchy's theorems [11]). As was mentioned in Introduction, this approach is applicable for low concentrations of inclusions (the upper index l in (7) is used to signify this fact).
In a similar way, the high-concentration rule can be obtained. The inclusions are now considered as "the host medium" and the host medium as "the inclusions", with the characteristic functionχ 0 = (1 −χ 1 ). A portion of "the inclusions" with the characteristic function ∆χ 0 = −∆χ 1 is embedded into that part of the current effective medium which is void of "the inclusions"; its characteristic function is ( Substituting (10) into (6) and taking the necessary integrals give the desired rule: Note that (11) is used rarely, in contrast to the original low-concentration ABM [13].

III. DIFFERENTIAL SCHEME WITHIN THE CGA
In this section, we develop a general procedure for building differential mixing rules based upon the CGA. The above low-(9) and high-concentration (11) rules turn out to be obtainable from the general differential equation under certain simplifications. In other words, these ABM rules are approximate and of limited usefulness. The general differential equation makes it possible not only to obtain their improved modifications, but also investigate the validity limits for the latter.
Substituting (16) into (6) and changing to the infinitesimal values, the differential equation is obtained. Actually, this is the differential form of (5). It is convenient to use (17) to derive new low-and highconcentration modifications of the ABM rules. Consider first the low concentration limit, where c and (ε 0 − ε) can be assumed to be of the same order of smallness as ∆c and ∆ε are. Then δε (h) ABM and the terms in the second brackets in (17) are of the second order of smallness and can be neglected. Correspondingly, δε CGA is determined only by the first and third terms in (16): and only the terms in the first brackets in (17) must be retained. This gives the differential equation This equation can also be derived by directly substituting (18) into (6), retaining the terms of the first order of smallness, and passing to the infinitesimal increments dc and dε. An analogous procedure for the high-concentration limit gives In general, equations (18) and (20) differ considerably from their ABM counterparts (8) and (10), but reduce to them provided the term δε CGA is neglected. This is possible if: (1) ε 0 ≈ ε and ε 1 ≈ ε, respectively; (2) the concentration of the constituent being added is small; (3) |ε 1 − ε 0 | is small as well. It follows that the original ABM mixing rules are, in general, physically inconsistent. In practice, one can attempt to apply them only to diluted low-contrast systems.
The equations (19) and (21) are improved differential equations, which partially take into account, through δε CGA , the interaction between the constituent. The integration of them results in the following mixing rules for low and high concentrations of the filler, respectively: Compared to the ABM, we expect these rules to be more accurate and valid for wider concentration regions. However, as based on the assumptions (18) and (20), they are still approximate. To prove this fact, consider the Hashin-Shtrikman (HS) upper ε + and bottom ε − bounds [33] It is easy to see that the rules (22) and (23) fail to satisfy these bounds (see Fig. 2). Indeed, consider the rule (22) for a system with ε 1 ≫ ε 0 . For c > (1 − e −1/2 ) ≈ 0.393, ε → ε 1 , which is higher than the upper HS bound (24) for the same concentrations. In the region of low concentrations (22) is closer to (5) than (9) and falls in between the HS bounds. Similarly, for (23) and c < e −2 ≈ 0.135, ε → ε 0 , which is lower than the HS bottom bound (25). For arbitrary ε 1 and ε 0 the concentrations where the HS bounds are violated depend on the contrast ε 0 /ε 1 . Figure 2 illustrates the situation for a system with ε 1 /ε 0 = 10 2 . It is interesting to note that the original ABM rules (9) and (11) satisfy the HS bounds. According to the above discussion, this fact is not indicative of the superiority of the ABM rules (9) and (11) over their modifications (22) and (23), but reflects the changes in the interplay between δε  ABM (r), and δε CGA (r) in the formation of ε that occur as c is changed. In other words, a simple extrapolation of the simplifications used within the differential method for one narrow concentration range fails to incorporate the effects essential for the formation of ε in other concentration ranges.
It should be noted that the above results quantitatively support the well-known qualitative arguments [18,26]  that at high concentrations, both the ABM and Maxwell-Wagner-Hanai models do not fully take into account interparticle polarization effects. They also explain why it is necessary to modify the classical differential models, or even introduce additional adjustable parameters, in order to extend their ranges of applicability [44,45]. And they are in accord the final-element calculations [46] which show that at small concentrations, the small changes of the effective permittivity due to the addition of a new portion of inclusions are greater than those predicted by the differential mixing procedure.

IV. CONCLUSION
The differential approach to the effective permittivity ε of a system implies that infinitesimal changes in the concentrations of constituents lead to infinitesimal changes in ε; the result is a differential equation for ε. In the present work, we develop a general differential scheme for dielectric systems of impenetrable balls based upon the compact groups approach (CGA). For this purpose, the CGA is formulated in a way that allows one to analyze the role of different contributions to the model permittivity distribution in a system. The analysis of these contributions and corresponding differential equations reveals that: 1. The low-and high-concentration mixing rules of the classical asymmetrical Bruggeman model (ABM) are reproducible within our model under the suggestion that the electromagnetic interaction between the inclusions already contained in the system with those being added can be replaced by the interaction of the latter with the current effective medium. Therefore, the classical ABM mixing rules are, in general, physically inconsistent and applicable only for diluted (with respect to one of the constituents) systems with low dielectric contrast between the constituents.
2. The overall changes in ε due to addition of an infinitesimal portion of one constituent include the contributions from both constituents (inclusions and the host medium) and depend on the state of the system before the addition. Ignoring the contribution from one of the constituents, we obtain generalized versions of the original ABM mixing rules.
3. The new generalized differential mixing rules are, again, applicable only in certain concentration ranges because beyond those they do not satisfy the Hashin-Shtrikman bounds. This means that different mechanisms are responsible for the formation of ε in different concentration ranges. Simple extrapolation of the results obtained for a certain concentration range cannot incorporate all the effects essential for the formation of ε in the whole concentration range.
The results obtained can be generalized to macroscopically homogeneous and isotropic systems with complex permittivities of the constituents 2 . However, it should be emphasized that the novelty of this article is not a derivation of new differential mixing rules, but the quantitative proof of a general statement that differential mixing rules are approximate and applicable only for narrow ranges of concentration and dielectric contrast. Attempts to go beyond those ranges will be, strictly speaking, misleading because of uncontrollable mistakes inherent to the differential method.
ACKNOWLEDGEMENT I thank Dr. M. Ya. Sushko for helpful discussions and encouragement and Prof. A. Lakhtakia for drawing my attention to incremental and differential Maxwell-Garnett formalisms. 2 It should be noted that, when applied to a system of particles with complex permitivities, the CGA is capable of giving quantitative estimates of such specific effects as interfacial polarization. Indeed, the CGA generalizes the ABM which, being a generalization (known as the Maxwell-Wagner-Hanai model) of the Maxwell-Wagner model (see Introduction) to complex-valued permittivities, automatically incorporates these effects. We plan to present these results elsewhere.

APPENDIX
Appendix A: Derivation of equation (6) The first step is to present the equation ∆E(r) − ∇(∇E(r)) + k 2 0 ǫ 0 ε f E(r) = −k 2 0 ǫ 0 δε(r)E(r) (A1) for a wave propagating in the inhomogeneous medium, with local permittivity ε(r) = ε f + δε(r), in the integral form Here: ∇ is the del operator; ∆ is the Laplace operator; the integral is taken over the volume V of the system under consideration; E 0 (r) = E 0 exp (i √ ε f k 0 · r) is a probing field with amplitude E 0 and wave vector k = √ ε f k 0 in the medium; ε f is the permittivity of an unknown background medium, in which each constituent (including the host medium) is embedded;T is the electromagnetic field propagator (the Green's tensor for (A1)). It can be shown [28,38,47] thatT can be associated with the tensor T, such that for "sufficiently good" bounded scalar functions ψ(r). In the long-wave limit (|k| → 0), the components of T satisfy the relation where δ(r) is the Dirac delta-function, and δ αβ is the Kronecker delta. The part T (1) characterizes the effects of multiple reemissions within compact groups [28]. Substituting (A3) into (A2), making elementary algebraic manipulations, and averaging statistically, we obtain the relation The local deviation δε here corresponds to macroscopic compact groups [28,30], not single particles (as in microscopic approaches, such as [4,31,32]). For macroscopically isotropic and homogeneous systems, the statistical average in the integrand depends only on |r − r ′ |. Because of this symmetry and a specific angular dependence of T (2) αβ , the second term in the right-hand side of (A4) vanishes, and (A4) reduces to E(r) = ξE 0 , ξ = 3ε f 3ε f + δε(r) . (A5) The average displacement D(r) is found using definition (1) and applying the above symmetry reasoning to the T (2) αβ -involving integral. The result is: where it was taken into account that ξ + η = 1.
Using (1) and taking into account (A5), (A6), and (A7), we obtain ε − ε f = (ε + 2ε f )η. (A8) In order to find the unknown ε f , we need one more relation between ε and ε f . We derive it using the equations between the electric field in vacuum and those in the effective and background media: whence, in view of (A7), ε − ε f = εη.
This and (A8) give the relation Of its two roots ε f = 0 and ε f = ε, only the latter is physically meaningful. It corresponds to the well-known Bruggeman-type homogenization and gives the equality η| ε f =ε = 0, that is, equation (6). It should be emphasized that for macroscopically homogeneous and isotropic systems, considered in the longwavelength limit, the homogenization condition ε f = ε is independent of the functional form of δε. The same result was obtained in [34] by combining the CGA with the Hashin-Shtrikman variational theorem [33]. This condition is the key relation postulated in many theories, including the strong-property-fluctuation theory [4], where it is used to get rid of the secular terms in the Born series for the renormalized electric field.