A mixing formula accounting for inversion of matrix structure

Known mixing models are analyzed with the aim to retrieve permeability of metal inclusions from the measured constitutive parameters of a binary composite. The application-oriented models are interpreted in terms of inclusion shape-factor and percolation threshold, which are two measurement-fitted parameters. A model that accounts for the inversion of the Maxwell Garnett matrix structure is proposed. The structure inversion point is close to the percolation threshold, and the inversion takes place within a transition filling range that is a third fitting parameter. The proposed model is compared with the effective medium model in terms of the complex susceptibility calculated as the function of filling and of frequency and in terms of Bergman-Milton shape-factor distribution charts. The model validity is illustrated by treatment of the measured microwave constitutive parameters of a composite filled with carbonyl nickel.Known mixing models are analyzed with the aim to retrieve permeability of metal inclusions from the measured constitutive parameters of a binary composite. The application-oriented models are interpreted in terms of inclusion shape-factor and percolation threshold, which are two measurement-fitted parameters. A model that accounts for the inversion of the Maxwell Garnett matrix structure is proposed. The structure inversion point is close to the percolation threshold, and the inversion takes place within a transition filling range that is a third fitting parameter. The proposed model is compared with the effective medium model in terms of the complex susceptibility calculated as the function of filling and of frequency and in terms of Bergman-Milton shape-factor distribution charts. The model validity is illustrated by treatment of the measured microwave constitutive parameters of a composite filled with carbonyl nickel.


I. INTRODUCTION
The interest in high-frequency and microwave properties of composites and their interpretation grows steadily. Most of natural and synthetic materials are composites. The study of their constitutive parameters may give valuable information on the structure of these composites, which is important to various disciplines such as medicine (estimation of bone porosity from impedance measurements 1 ), chemistry, 2 seismology, and petrophysics [3][4][5][6] and to high-frequency applications (antennas, electromagnetic compatibility materials, and microwave absorbers [7][8][9][10]. The research aim is to find the mixing model (the formulas relating the properties of a binary mixture with those of its components and their volume fractions p) suitable to retrieve the microwave permeability of metal inclusions from the measured constitutive parameters of composites filled with metal powders. The wider the filling range of composites under study, the higher susceptibility of inclusions may be retrieved accurately by any model (see Sec. II). The difference between mixing models shows itself at high filling. For a reliable selection of a valid model, the filling range should exceed the percolation threshold; therefore, the model should account for possible transformation of composite structures with filling increase.
To simplify the comparison, sorting, and possible application of mixing models to various disciplines, we analyze mixing formulas in terms of the normalized susceptibility χ. Therefore, the approach and formulas may be readily adjusted to treatment of any "generalized conductivity" of composites, i.e., permittivity ε, permeability μ, conductivity σ, diffusion coefficient, and Young's module. 11 The susceptibility of component i is normalized vs the "generalized conductivity" of the binder marked as component 2. Particularly, for permittivity, χi 2 = εi/ε 2 − 1; permeability, χi 2 = μi/μ 2 − 1; and an impermeable binder, μ 2 = 1. The higher the normalized susceptibility, the higher the contrast between the corresponding properties of mixture components.
Microwave electromagnetics reveals additional demands for model consistency as the formulas are applied to describe the dynamic susceptibility while all mixing models are quasistatic and omit the scaling factor (i.e., the size of inclusion or of the inclusion cluster normalized to the wavelength). Usually, the factor is accounted for by attributing the effective susceptibility to an inclusion. The effective susceptibility accounts for skinning large particles 15,20 and inductance of long fibers. 13,14 Calculations of effective susceptibility do not depend on the mixing model and fall out of the research scope.
The developed model is verified by comparison with the models derived for matrix and symmetric mixtures, by analysis of Bergman-Milton geometrical spectral functions (see Sec. V), by analysis of high-frequency susceptibility spectra specific for each mixing model, and by comparison with the measured complex susceptibility of composites filled with carbonyl nickel. 15 The number of published mixing models is difficult to count, so to clarify the model grouping, we have to define the principles of formula selection. The formulas should describe the complex constitutive parameters and their frequency dispersion for the mixture, while the normalized susceptibility of inclusion can vary from zero to infinity. The formulas should account for inclusion shape (defined by its shape-factor N) and for the percolation threshold (the critical filling pc where inclusions agglomerate into infinite-long clusters). The formulas should include the minimal number of physically clear parameters fitted from the measured data.
We do not consider the numerical computation algorithms (an overview of computational electromagnetics approach to composites is given in Ref. 19) that do not result in an explicit formula and do not consider the special formulas applicable to particular inclusion shapes, 12 filling ranges, and particular contrast of component properties. 21 The published mixing models approximate the inclusion shape by ellipsoids. Depending on volume fraction p, one or both components may be considered as prolate or oblate ellipsoids of rotation; the ratio of the ellipsoid axis to that of rotation may be close to zero (flakes), to unity (spheres), or to infinity (fibers). The limitation on distortion from the spherical shape is imposed by the approximate account for mixture anisotropy: the effective filling (i.e., the fraction of inclusions interacting with the exciting field) account is accurate for particles with N ≈ 1/3 and for ellipsoids with the shape-factor N close to 0 and 0.5 only. The analyzed models describe isotropic composites for spherical inclusions (N = 1/3) only. In the case of fibers (N = 0), the models describe the anisotropic one-dimension composite, whereas in the case of flakes (N = 1/2), the models describe the planar two-dimensional composite. For the isotropic fiber-filled composite, the effective filling factor is equal to 1/3 of the real volume fraction of fibers p eff = p fiber /3; similarly, for the isotropic flake-filled composite, p eff = p flake /2. We use the above approximate account for mixture anisotropy, 15 as the rigorous one is yet to be developed.
To calculate the properties of a composite based on the corresponding properties and volume fractions of its components (the direct mixing problem), one needs to know that the mixing model is valid for the composites under study. The knowledge of the valid mixing model is necessary to relate the measured constitutive parameter of a composite to the constitutive parameter of one of its components (the inverse mixing problem) as well.
There are many reviews devoted to the analysis and sorting of the known mixing models. [16][17][18] The mixing formulas are usually compared based on physical assumptions the corresponding models are built upon, while the number of fitting parameters falls out of considerations. The fitting parameters are obtained from the measured data and are usually related to the effects neglected by the rigorous first-principles models. Therefore, the verification and validity range of the reviewed mixing models are undefined.
The model verification is much more complicated than a simple comparison with a set of experimental data. An invalid model may be fit to interpolate the dependence of a measured parameter on filling, being at the same time nothing more than an empirical formula with some formal parameters. An empirical formula has an undefined applicability range and is inaccurate for describing a wide range of component properties. Dealing with the problem to retrieve the microwave permeability of metal inclusions, we see that the microwave permeability (and consequently the property contrast of mixture components) varies due to frequency dispersion from unity to several hundreds or higher depending on the particular permeable alloy. The fitting parameters of an invalid formula vary with the frequency similar to the permeability contrast, while the parameters of a valid mixing model should not depend on the frequency. In the measurements of dynamic susceptibility, the fitted parameters vary with the frequency for any model because of measurement errors as well as because of the model inaccuracy. Therefore, selection between the valid and invalid models is formally the problem of mathematical statistics (the higher is the variation of mixing parameters with the change of susceptibility contrast for the same number of fitting parameters, the less probable is the tested model), which is of low relevance for practical needs.
The more practicable approach to the model verification is proposed in Ref. 15. The dependence of the inclusion permeability μ incl on frequency f is retrieved with a tested model treating both permittivity and permeability of composites, measured in a wide frequency range. The samples are filled with different volume fractions p of permeable powder; the filling factor p varies from zero to maximum, limited by the preparation technology. The less is the variation of the retrieved permeability dependence on frequency μ incl (f ) with the filling factor p change, the closer is the tested model to reality and the closer is the retrieved inclusion permeability μ incl to its real value.
To develop the above verification procedure, we need at first to systematize and to review critically the known mixing models to be able to outline their limitations and applicability range.
Literature study reveals more than 50 mixing rules. The number of original formulas is much smaller since the same formula is derived several times, e.g., the Lichtenecker and the Legg and Given formulas 41,56 are identical, the Maxwell-Garnett formula written for ellipsoidal inclusion is named the Ollendorff formula, 55 the Bruggeman formula 22 is named the Polder-van Santen, Das-DasGupta, or Rao formulas. [23][24][25] Some formulas are changed while being cited or are rewritten for a particular case, e.g., the Wiener formula is commonly written as Here, χ 12 = ε 1 /ε 2 − 1 and χ 32 = ε 3 /ε 2 − 1 are susceptibilities of component 1 with volume fraction p 1 and of the mixture (named component 3), respectively. The susceptibilities are normalized to permittivity ε 2 of component 2 with volume fraction p 2 = 1−p 1 . Equation (1) describes an anisotropic composite where the fibers or layers are parallel to the external field. The formula is similar to capacity calculation for the parallel connection of two capacitors. The first capacitor presents component 1 with permittivity ε 1 and volume fraction p 1 ; the second capacitor presents component 2 ARTICLE scitation.org/journal/adv with permittivity ε 2 (with the normalized susceptibility equal to zero: χ 22 = ε 2 /ε 2 − 1 = 0) and volume fraction p 2 . Below, we quote the original paper 26 in detail as the Wiener formula in its general form [Eq. (2)] is used further to develop a novel mixing model that accounts for structure inversion.
The Wiener formula [Eq. (2)] describes the permittivity dependence of anisotropic binary mixture ε 3 on volume fraction pi and permittivity εi of the i-component. Compared to the original form, 26 Eq. (2) is rewritten in terms of susceptibility χi 2 = εi/ε 2 − 1 normalized to permittivity ε 2 of component 2. The right-hand side of Eq. (2) is generalized for ellipsoidal inclusions, (2) Here, factor U depends on the local permittivity ε local in the inclusion vicinity; indices 1, 2, and 3 refer to component 1, component 2, and the mixture itself (component 3), respectively. Note that in contrast to other mixing formulas, the local permittivity ε local may differ from permittivity of mixture ε 3 or of its components ε 1 and ε 2 and may depend on the shape-factor N of the inclusion.
It is easy to see that if U = 0, Eq. (2) turns into Eq. (1). If U = ∞, Eq. (2) turns to Eq. (3), which describes the same anisotropic composite as Eq. (1), but here the component layers are oriented perpendicular to the external field, Equation (3) is similar to capacity calculation for the series connection of the same capacitors as for Eq. (1). For an anisotropic composite, Eqs. (1) and (3) describe the highest and the lowest permittivity limits similar to Hashin-Shtrikman static 27 or Bergman dynamic 28 limits for an isotropic mixture.
Wiener has defined the value of U-factor only for the extreme shapes of the inclusion, such as infinite fiber (N = 0, U = ∞ or N = 1 2 , U = ε host ), infinite plane (N = 0, U = ∞ or N = 1, U = 0), and sphere (N = 1 3 , U = 2ε local ), but it is easy to generalize Wiener's considerations upon the shape of ellipsoid of rotation: U = (N −1 − 1)ε local . Therefore, the local susceptibility χ local2 normalized to the permittivity of component 2 may be defined by the following equation: Equation (4) Therefore, Eq. (2) with Eqs. (4) and (5) (2)] is that the increase in U-factor from zero to infinity describes the gradual transformation from the parallel to series circuit formed of mixture components. This feature is used (see Sec. III) to develop a novel mixing formula with the inversion of matrix structures.
Further analysis of mixing models is limited to anisotropic binary mixtures where one of the components is treated as an inclusion with a definite shape approximated by ellipsoid of rotation; the second component is regarded as a matrix or binder and its shape is undefined.
Mixture anisotropy is defined by the shape of ellipsoids with the long axis oriented parallel to the external field; the random orientation of ellipsoids in the case of isotropic mixture is accounted for by the effective filling factor that includes the fraction of ellipsoids with the longest axis oriented parallel to the external field. 11 The account is rough for slightly oblate or slightly prolate ellipsoids, but it is accurate enough for practical cases of fibers and flakes. 12,13,15,30 The concepts of matrix and inclusion may be relative as the component attribution to the binder or filler may depend on filling, e.g., for components with similar mechanical properties. A mixture may be considered similar to a binary AB system with limited solubility of component A in B and B in A in the solid state. 31 At low filling (pa → 0), the mixture is similar to a single-phase solid solution of component A in B, where the inclusions (A-component) of known shape are placed in a binder medium (i.e., the matrix mixture described by the Maxwell Garnett model 32 ). Close to the percolation threshold, the mixture is "symmetric" and resembles a two-phase zone of solid solutions AB and BA in a phase diagram. Within the "two-phase zone," both components are equal and have definite shapes (the Bruggeman 22 and the Sihvola 33 models). At the inclusion concentration close to unity (A-component volume fraction pa → 1), the matrix structure is inverted, resembling the solid-state solution of component B in A: the inclusion component turns into the binder medium with an indefinite shape and the matrix component is localized in the pores with definite shape (i.e., the inverted Maxwell Garnett model 32 ). The filling range of matrix and symmetric zones depends on many factors (mutual wettability of components, mixing technique, surfactant additions, and mechanical properties) that fall out of idealized first-principles considerations lying in the basis of the first known mixing models. These factors are taken into account at least partially by the parameters fitted from the measured dependence of composite susceptibility on the filling. Therefore, the comparison of mixing models should account for the physically rigorous model the formula is derived for as well as for a number and physical interpretation of the free parameters fitted from the experimental data.
We try to consider the mixing models that account for the majority of effects relating the volume fractions and component properties to the properties of a composite. Still, there is a list of important parameters that are missed by all models, e.g., the granulometric distribution of inclusion size and shape, the contribution of inclusion oxide shell, or vice versa the contact resistance of neighboring inclusions. The variation of the above parameters can affect the formal parameters of a mixing model ARTICLE scitation.org/journal/adv
The article structure is as follows: analysis of the known mixing models and development of the mixing formula with 3 fitting parameters, comparison of the Bergman-Milton spectral density functions for the proposed and for the known 2-parameter formulas, and application of the proposed formula to retrieve the permeability of Ni-spheres from the measured constitutive parameters of a composite.

II. ANALYSIS OF ARCHIVED MIXING FORMULAS FOR BINARY COMPOSITES
The known mixing models are compared below to select the most general ones, to reveal their limitations, and to find the approach amending them.
The original mixing formulas satisfying the principles listed above are presented in Tables I-III. Table grouping is performed based on the number of free parameters, on the composite type structure (the table columns), and on the formula derivation principle. The right column (Tables I and II) presents the formulas obtained by differentiation for the case of a low contrast between the component permittivity; the low-contrast mixing models are grouped separately, as they are low sensitive to the composite structure.
To simplify the formula comparison and analysis, all formulas are rewritten here in terms of susceptibility χ normalized to the constitutive parameter (e.g., permittivity, permeability, and diffusion coefficient) of component 2 nominally considered for the binder with undefined shape. The indices 1, 2, and 3 refer to inclusion, binder, and mixture, respectively; thus, χ 32 and χ 12 are the mixture and the inclusion susceptibilities, respectively, normalized to that of the binder.

Matrix mixtures Symmetric mixtures
Odelevskiy 1947 Sihvola 1999 The models with nonspherical inclusions are simplified for anisotropic composites with collateral inclusions. Table I presents the oldest first-principles formulas derived for composites with the spherical inclusions; the models do not have free parameters and do not account for any deviation from the idealized hypothetical structures they are derived for. The Maxwell-Garnett and Bruggeman models have a priori limited application and accuracy, as even the spherical inclusions may aggregate into clusters with effective shape that differs from the sphere. The above formulas present, to our viewpoint, historical interest only; nevertheless; they were rederived multiple times. [34][35][36][37][38] The differential based exponential 39,40 and logarithmic 41,56 models are widely used in petrophysics as the semiempirical Archie and Topp formulas, 42,43 despite the low-contrast limitation of the component susceptibility. The exception from the low-contrast limitation is the asymmetric Bruggeman model 22 generalized by Hanai 44 for dynamic properties. The asymmetric direct and inverse formulas are derived 45 by the differentiation of the low-filled and high-filled symmetric Bruggeman formula. 22 Table II presents the formulas modified to account for the inclusion shape. The exception is the Musal-Hahn-Bush 46 formula, where the fitted parameter P fit is a formal one without the distinct physical sense. Analysis of the dynamic susceptibility shows that the series application of matrix and symmetric mixing formulas 46,47 leads to unexplainable peaks of absorption; 48 therefore, the formula is excluded from considerations.
The inclusions are usually approximated as ellipsoids of rotation; the ellipsoid polarizability is described by the depolarizationor shape-factor N [Eq. (4)]. In contrast to the models derived for isotropic mixtures with spherical inclusions, the formulas in Table II describe the susceptibility of anisotropic composites with ellipsoids oriented parallel to the electric field. These shapes and orientations define the anisotropy of the composite. The practicable composites are usually filled with either very prolate ellipsoids (needles or fibers) or vice versa very oblate ellipsoids (disklike scales or flakes). Both shapes have one of the depolarization factors much smaller than the other two factors, i.e., the polarizability along the fiber axis or along the flake plane is much higher than that along other dimensions. This difference justifies the effective filling account for random orientation of elliptic inclusions in isotropic composites.
The generalized Maxwell Garnett 11,21,32 and the symmetric Bruggeman 22 models are the best known in this group. The Maxwell Garnett formula can be written in two forms. The well-known one [Eq. (6)] describes the direct matrix mixture, The inverse Maxwell Garnett formula [Eq. (7)] describes the mixture where the roles of host and guest components are inversed (the inclusion component is assumed to be the host one 50 ), The problem with the generalized symmetric Bruggeman model [Eq. (8)] is that the generalization of the original formula onto elliptical inclusions is shown to be physically contradictive (both components can have the same shape-factor N only in the case of spheres). 30,44,49 As a result, the percolation threshold here is defined totally by the shape-factor of the inclusion (pc = N), The less known coherent potential model 45,46,52 gives the dependence of properties on the filling, which is close to the symmetric Bruggeman model, and the percolation threshold here is also defined by the inclusion shape-factor N. Difference of Eq. (9) from Eq. (8) is in the percolation threshold, which appears to be lower than the shape-factor [pc = N/(1 + N)], The differential based exponential and logarithmic models amended by Boyle 50 and Rother 41 to account for higher polarizability of ellipsoids than that of spheres (Table II, (Table I). The limitation leads to the zero-filling percolation threshold evident for permittivity calculation with resistive fillers. Obviously, the exponential and power-based models are unable to describe the dielectric properties of composites with metal powders.
The last differential based model is the Bruggeman-Hanai 44 model amended by Bussian 51 for ellipsoidal inclusions. The formula inherits the percolation threshold behavior from the symmetric Bruggeman model. 22 Similar to the direct and inverse Maxwell Garnett formulas, there are two asymmetric Bruggeman-Bussian formulas. The first one [Eq. (10)] is obtained by susceptibility differentiation for a low-filled composite, and the second formula [Eq. (11)] is obtained by susceptibility differentiation for a high-filled composite, ARTICLE scitation.org/journal/adv The asymmetric Bruggeman-Bussian model inherits the drawbacks of the symmetric model and "ignores the system structure before the addition of one constituent." 45 The drawback of the models (Table III) that account for inclusion shape and display the percolation threshold effect (the structure transformation of a composite at the threshold) presented is a neglect of the component interaction; consequently, the percolation threshold pc is defined here by the inclusion shape only. Table III presents the mixing models with two fitting parameters. The first one is the inclusion shape-factor N accounted for in Eqs. (6)- (11), and the second parameter is the percolation threshold pc that depends on a number of factors besides the inclusion shape that are difficult to enumerate. The simplest model is the Odelevskiy formula, 11,15,20 This formula is, in fact, the Maxwell-Garnett formula [Eq. (6)] renormalized for the filling range lower than the percolation threshold (the filling where the structure is still that of the matrix mixture). The formula drawbacks are the limited filling range p < pc, the assumption of a regular structure with identical interaction of all inclusions, and the low-contrast asymptotic that diverges from the rigorous LLL formula 39,40 in the threshold vicinity. Above the threshold (p > pc), the mixture is assumed 11 to be described by the popular version of the Wiener formula [Eq. (1)].
The approach proposed by Odelevskiy is the sectional stitching of the matrix and symmetric models with the stepwise switch between different structures at the percolation threshold pc. Close to the threshold, the accuracy is obviously poor as the structure transformation is gradual in reality.
The more advanced Sihvola model 33 provides the formula (Table III) where the first free parameter N is the inclusion shapefactor and the second free parameter a Sihvola defines as a formal structure parameter (0 < a < ∞). We have shown earlier 15 that the formal parameter a is related to the critical filling as a = N(pc −1 − 1). Therefore, the original formula (Table III) can be rewritten in the same terms of N and pc that have clear physical meaning. To simplify the analysis and comparison with other formulas, the formula (Table III) is reduced for the case of anisotropic composite with regularly oriented inclusions to the form The Sihvola formula [Eq. (13)] is the most general of the formulas listed in Tables I-III Comparison of the Sihvola formula [Eq. (13)] with the generalized symmetric Bruggeman formula [Eq. (8)] shows that the difference is in the renormalized critical filling pc = N/(a + N). It may be readily shown that in the case of the Bruggeman formula, pc = N for both components, while the joint solution of Eq. (7) and Eq. (13) shows that at the fillings p < pc, the shapeless binder may be described by the shape-factor N 2 , which is equal to the critical filling N 2 = pc. Therefore, the Sihvola model is indeed the symmetric Bruggeman model with a rigorous account for the ellipsoidal inclusions.
Selection between the Odelevskiy and the Sihvola models [Eqs. (6) and (8)] from the measured susceptibility data is a difficult task, 15 as the difference between the models becomes visible at the filling close to and above the percolation threshold pc. At high fillings, the sample preparation is an engineering problem (the filling increase above pc leads to unregulated porosity and the binary mixture turns actually into a three-component composite) and the contribution of unaccounted for interfacial effects (contact resistance and insulation shell) becomes high.
We omit from consideration the known mixing models with three fitting parameters, namely, the exponential effective medium model proposed by McLachlan 53 and Rozanov. 54 The reasons are the analysis of high-frequency susceptibility spectra (the McLachlan formula does not describe the frequency dispersion of the permittivity specific for resistive composites) and the vague physical meaning of two parameters related with critical indices. The model 54 is applicable to inclusions close to spherical only.
Looking at the above classification tables is that it is possible to distinguish three groups of mixing models (Tables I-III); each group describes a single structure that, in fact, can be formed within some limited filling range. The drawback of the direct Maxwell Garnett, inverted Maxwell Garnett, Odelevskiy, Bruggeman, and Sihvola models is that each model is applicable to a single mixture structure out of the three possible structures (see Sec. I, analogy with binary phase diagrams). The visual similarity between the low-filled direct matrix-type and symmetric-type models is superficial; the attempt to describe the permittivity and permeability (the parameters with different contrast of component properties) for the same composite shows the difference. 15,20 Nevertheless, this visual similarity reveals the fundamental limitation on the value of retrieved susceptibility χ 12 inherent to all diluted composites: the lower the filling p and the higher the shape-factor N, the lower the inclusion susceptibility can be retrieved from χ 32 -measurements accurately: χ 12 < (1 − p) −1 N −1 (see, e.g., the permeability retrieval error in Ref. 20).
Therefore, the wider the filling range where the mixture structure is close to that of the model, the more accurate the retrieved susceptibility of inclusions. The problem is that the models valid in a wide filling range should describe several mixture structures and account for structure transformation. Such models are still undeveloped.

III. THE UNIFIED MIXING FORMULA
It is possible to develop the Odelevskiy idea of stepwise stitching 11 of the matrix and of the symmetric structures at p = pc more accurately. In this case, Eq. (6) should be pc-renormalized for the filling range 0 < p/pc < 1; the difference from Eq. (12) is that χ 12 should be substituted with χ 12 × pc. The overcritical filling p > pc should be ARTICLE scitation.org/journal/adv described with Eq. (14) instead of Eq. (2). However, the advantages of stepwise stitching are dubious as the joined structures are too different and the accuracy will be poor close to the stitching point; therefore, the structure transformation should be gradual, while the stitched structures should be the limiting ones (should exist at p → 0 and p → 1). The proposed model is a development of the Odelevskiy sectional joint of two structures at p = pc and of the Sihvola unified mixing approaches; the development is based on the Wiener concept of the local host medium. The advantage of the Wiener formula [Eq. (2)] is that it is readily turned into all the considered formulas [except for the differential-based Bruggeman-Bussian ones, Eqs. (10) and (11)] that account for a wide range of component property contrast and of inclusion shapes. Specifically, if Nχ 12 → 0 Eq. (2) is simplified to Eq. (1) that is a onedimensional LLL formula; 39 if the binder is assumed to be the local media (χ local2 = 0), Eq.
(2) turns into the direct Maxwell Garnett formula with the percolation threshold equal to unity (pc = 1); if the filler is assumed to be the local media (χ local2 = χ 12 ), Eq. (2)  The above transformations of the Wiener formula show that the published attempts to account for inclusion interaction are developed for some imaginary mixture structures and lead to the listed above dependences of the effective local susceptibility χ local2 on component properties, inclusion shape, and critical filling that are difficult to interpret. Similar to Ref. 53, we do not try to derive rigorously another idealized mixing model but consider that empirical combinations of the known mixing models 46,47 are physically unjustified. We assume that a binary composite has two limiting mixture structures, namely, the direct [Eq. (6)] and the inverted [Eq. (7)] matrix structures formed at low content and at high content of component 1, respectively. All known mixing models ignore the inversion of matrix structures.
Developing the unified mixing approach, we stitch these limiting structures with a single formula. Contrary to the Odelevskiy model with the stepwise stitch of different structures [Eq. (12) for p < pc and Eq. (1) for p > pc] at the percolation threshold filling pc, we consider a gradual transition between the limiting structures. The symmetriclike structure is assumed to be formed in a transition filling range Δpc between the direct and the inverted matrix structures; the local medium susceptibility in this range should be a mean value between the host susceptibilities of the limiting matrix mixtures. In contrast to the effective media models, 22,33 we consider the local susceptibility χ local2 to differ in a wide filling range from that of the mixture χ local2 ≠ χ 32 , but χ local2 should go up similar to the effective media models with filling p increase. At p ≪ pc, the local medium is the matrix component χ local2 = 0; at p ≈ pc, the local medium should be close to the mixture χ local2 ≈ χ 32 , but at p ≫ pc, the local medium is formed accounting for the filling factor by the former inclusions χ local2 ≈ p × χ 12 . The filling range of structure transformation Δpc depends on many factors falling out of idealized rigorous models. Therefore, we have to select the susceptibility averaging function of filling p that changes gradually the local susceptibility from χ local2 = 0 at p → 0 to χ local2 = χ 12 at p → 1 with the transition range Δpc and the inversion point pc as free parameters fitted from the measured data.
To account for the structure inversion, we take the generalized mean (other averaging functions exclude the percolation phenomenon) between the direct χ local2 = (1 − p) × χ 22 = 0 and the inverted χ local2 = p × χ 12 matrix mixtures as the susceptibility averaging function to calculate the local susceptibility, The weight factor w is a function of the distance from the structure inversion point |pc − p| and describes here the gradual transition between the direct and inverted structures. For w = 0, the transition is stepwise and abrupt; for w = 1, the matrix structures are not formed at all, as the transition zone Δpc occupies the whole (0 < p < 1) filling range.
The filling p attributed to a particular sample is commonly a mean value for the sample under study. The mixture structure at the scale related to the distance between neighboring inclusions is defined by the local filling that fluctuates around the mean value p. The mean value of fluctuation amplitude may be described by the fluctuation factor δ defined by a mixing technology (0 < δ < ∝). The closer the mixture to a pure component (p → 0 or p → 1), the lower the fluctuation factor affects the mixture structure; the product (1 − p)p accounts for the decrease.
Below, we assume that the RMS deviation | p c −p p(1−p) | 2 from the inversion filling pc is defined by the Gaussian distribution φ, where the fluctuation factor δ is the standard deviation, Therefore, the structure transition range Δpc is a statistical value depending on the structure inversion filling pc and fluctuation factor δ. The transition range Δpc may be defined as the filling range where the local filling deviates from the inversion point pc with 95% probability Δpc ≈ 2δ(1 − pc)pc.
The weight factor w in Eq. (14) is the cumulative distribution of deviations from the inversion point pc,

ARTICLE scitation.org/journal/adv
Here, t is the integration variable. Note that Eq. (16a) is obtained for a Gaussian distribution for filling deviations from the inversion point pc, and other distribution functions may be applied as well. Moreover, the cumulative distribution for w(p) may be substituted with a simpler function satisfying the limitations: w(0) = 0, w(1) = 1, and w(pc) = 0.5, which approximates the structure transition, e.g., Accounting for the factor normalizing the depolarization coefficient and the percolation threshold ( N p c 1−p c 1−N ) (the factor is the same as for the Sihvola formula), we can rewrite Eq. (14) as Substituting the local susceptibility χ local2 into the Wiener formula [Eq. (2)] by the value obtained with Eq. (17), we get the mixing formula that accounts for the inversion of matrix structure, Here, χ 32 and χ 12 are mixture and inclusion susceptibilities, respectively, and p is the volume fraction of inclusions. The shape of inclusion is described by the depolarization coefficient N similar to other formulas, pc is the structure inversion filling, and t is the integration variable. The weight factor w(p) describes the smoothness of structure transformation due to filling fluctuations and defines the width of transition range Δpc.
Simulation  (18)] the structure inversion filling pc identically to the percolation threshold filling for Eqs. (12) and (13). The physical meanings of structure inversion and percolation threshold are close to each other, but here the structure transformation is described with two parameters, and the second parameter δ is related to the width of the transition range Δpc = 2δ(1 − pc)pc. The threshold filling is equal to the structure inversion filling only for Δpc → 0; the wider the transition range Δpc, the higher the structure inversion point compared to the percolation threshold. The latter is the filling where the first infinite cluster of inclusions is formed, while the former is the filling where the inversion of matrix structures takes place.
The validity of the resulting mixing formula [Eq. (18)] is verified further by a numerical comparison (Sec. IV) with the Sihvola model [Eq. (13)] at the fillings close to the percolation threshold by analysis (Sec. V) of the Bergman-Milton geometrical spectral function 36 and by comparison (Sec. VI) with the measured data for dynamic susceptibility of a metal-filled composite.

IV. FILLING AND FREQUENCY DISPERSIONS OF COMPLEX SUSCEPTIBILITY FOR THE MATRIX-INVERSION AND FOR THE SIHVOLA MODELS
We assume the Sihvola formula [Eq. (13)] to be the most rigorous to describe the symmetric-type mixtures; therefore, for p → pc, Eq. (18) should give susceptibility values close to the ones calculated with Eq. (14). The formula comparison is performed via computer simulation. The contrast of component properties is selected high enough (χ 12 = 0 + i × 100) to reveal the effect of inclusion shape on the mixture susceptibility. The calculated dependence of the mixture susceptibility on filling for the Sihvola and for the proposed formulas is shown in Fig. 1.
The inclusion shape-factor N = 0.24 and the percolation threshold pc = 0.33 (the structure inversion filling for the developed model) are close to the values obtained for sphere-filled composites. 15,20 The filling fluctuation factor δ = 0.5 is taken to obtain a reasonably wide transition zone Δpc with the susceptibility curve close to that of the Sihvola model.
The calculated dependences of the real χ ′ 32 and imaginary χ ′′  (13) is that χ ′′ 32 (p)shows [ Fig. 1(b)] three regions: within the direct matrix zone, χ ′′ 32 is lower than that for the Sihvola formula; within the inverse matrix zone, χ ′′ 32 exceeds the value for the Sihvola formula, while within the transition zone, the slope of χ ′′ 32 (p) is higher than that for Eq. (13). The values of χ ′′ 32 are approximately the same near pc for both models. The filling factor hodographs (Argand diagrams) of complex susceptibility χ 32 for the same composite with χ 12 = 0 + i × 100 are shown in Fig. 2. The dots in the hodographs [ Fig. 4(a)] mark the pcvalue that corresponds to the percolation threshold for the Sihvola model (gray circle) and to the structure inversion point (black circle) for the proposed model. The left graph displays calculation results for Eq. (18) (black lines) and for Eq. (14) (the gray line). In the case of the Sihvola model, the width and height of χ ′ 32 -peak are related to the contrast of component properties only; in the case of Eq. (18), the peak shape is defined by the δ-value, which is an experimentally fitted parameter. The higher the susceptibility χ 12 , the lower δ should be for the better agreement with Eq. (14) and the closer the critical filling [pc in Eq. (14)] for the Sihvola model is to the structure inversion point for the developed model [pc in Eq. (18)]. Figure 2 The difference between the Gauss weight factor [Eq. (16a)] and the simplified approximation [Eq. (16b)] becomes visible for δ > 4 where the transition range Δpc covers the filling 0 < p < 1 and the matrix structures are not formed at all. Therefore, the weight factor w(p) calculations with Eq. (16b) being much simpler than with Eq. (16a) do not lead to accuracy loss.
The susceptibility dependence on frequency for the Sihvola (gray lines) and for the developed (black lines) models is shown in Fig. 3. The gray dashed line corresponds to equal values of the percolation threshold and the inversion filling.
It is easy to see that by fitting the parameters pc and δ, the susceptibility dispersion curves χ 32 (f ) for a single arbitrary filling p can be made identical for both models (the gray solid and black solid curves). The difference reveals itself (Fig. 1) studying a range of fillings, where the parameters are constant for a true model and vary with filling for a false one.

V. ANALYSIS OF BERGMAN-MILTON SPECTRAL FUNCTIONS
There are several methods to calculate the shape-factor distribution density function. 36,40,54 In the case of the Wiener mixing formula [Eq. (2)], the function may be calculated analytically using the Ghosh-Fuchs technique 57 from the composite susceptibility, where the susceptibility of a cluster inclusion χ 12 = −(n + is) −1 is defined by the running shape-factor n and resistance s, Here, N is the shape factor of a free inclusion with resistance s and volume fraction p, and n is the running shape-factor, which is  the argument of the shape-factor distribution density spectral function. Equation (19) describes the formation of inclusion clusters as a function of mixture structure; the most probable [corresponding to the maximum of B(n)] shape-factor of clusters differs from the shape-factor of a free inclusion particle and depends on filling and the sample structure. The Wiener formula [Eq. (2)] is readily turned into other mixing formulas (Sec. III) by substituting the appropriate value for the local susceptibility; therefore, Eq. (19) may be applied to calculate analytically the dependence of shape-factor distribution density on filling for the mixing models derived from Eq. (2).
We calculate the dependence of the shape-factor distribution density for the Sihvola model on filling and on running shape-factor by substituting χ local2 into Eq. (19) The result is presented in Fig. 4 as a 3D plot, as a contour map, and as cross section of the 3D plot at the selected filling factors (p = 0.01, 0.33, and 0.6).
The spectral function of the shape-factor distribution density for the proposed model is calculated similarly by substitution of the local host susceptibility χ local2 in Eq. (19) with Eq. (17) [and Eq. (16a) for the weight factor w(p)].

ARTICLE scitation.org/journal/adv
For the sake of comparison, the shape-factor distribution density for the developed model [Eq. (18)] is presented in Fig. 5 similar to that of the Sihvola model in Fig. 4.
Comparing the Bergman-Milton function maps for the Sihvola (Fig. 4) and for the proposed (Fig. 5) models, it is possible to see that for p < pc, both models show the decrease in effective shapefactor N with filling increase due to cluster formation. The difference is that the cluster shapes are almost identical for a matrix mixture compared to the shape spread for a symmetric mixture.
At p = pc, the infinite long (with N → 0) clusters parallel to the exciting field are formed in both models. The difference is that for the symmetric mixture, the fraction of unclustered inclusions of initial shape is negligible compared to the proposed model.
The increase in filling above pc leads to simultaneous formation of the infinite long clusters (N = 0) and of the shorter fragments widely distributed in shape for the Sihvola model. In the case of the proposed model, the filling increase close to pc leads to the formation of two distinct peaks of effective shape-factor N. The intensive wide peak corresponds to long clusters with N → 0, but it is much wider than for the Sihvola model. The position of the second peak is related to the initial shape factor N of free inclusions and to the value of structure inversion filling pc; the shape-factor is higher than that of free inclusions and the peak is much sharper than for the Sihvola model. The phenomenon may be interpreted as orientation of some clusters perpendicular to the exciting field. The interpretation seems reasonable as the shape distribution here is uniform compared to the shape spread for the Sihvola model.
The map of shape-factor distribution density for the proposed model (Fig. 5) confirms the percolation threshold behavior (see Fig. 1) that is inherent to the symmetric-type effective-media models though the proposed model is a combination of matrix models.

VI. COMPARISON WITH EXPERIMENTAL DATA
Calculations with the Sihvola and with the proposed models [Eqs. (13) and (18)] are compared with the constitutive parameters measured for the paraffin bound composite (ε 2 = 2.3) filled with carbonyl nickel powder; the preparation and measurement techniques are described in Ref. 15. The measurements of static permittivity ε 3 and complex permeability as a function of frequency μ 3 (f ) were performed for a set of samples filled with spherical particles (∅12 mcm) up to about p ≈ 0.35; 15 the measurement of static permittivity at higher filling was problematic. The measured data made it possible to determine by fitting the effective shape-factor of the nickel inclusions (N = 0.2) and the percolation threshold (pc = 0.43).
Here, we approximate the measured dependence of the complex permittivity on frequency ε 3 (f ) for an identical composite filled with carbonyl nickel as close to the percolation threshold as possible. The filling factor (p = 0.375) is limited by the brittleness and porosity of the cold-pressed washer. The values N = 0.2 and pc = 0.43 are taken from Refs. 15, and the fluctuation factor δ is fitted for the measured dependence of dielectric susceptibility on frequency, The measured permittivity for a sample with p = 0.375 of Ni 15 is presented in Fig. 6(a) as a susceptibility (χ 32 = ε 3 /ε 2 − 1) hodograph with a dotted line; the frequency range is 0.15-20 GHz. The gray solid line presents the Sihvola model with N = 0.2, pc = 0.43, and inclusion conductivity σ 1 = 10 6 Ω −1 m −1 ; the black solid line presents Eqs. (18) and (16a) for N = 0.2, pc = 0.43, and δ = 1.
The measured dielectric susceptibility looks to be in good agreement with Eq. (18) than with Eq. (13), but the discrepancy with the measured data here depends strongly on accuracy of the fitted parameters. The formula selection for the filling range p < pc (p = 0.375 < pc = 0.43) is problematic (see Fig. 3). The reason is that the selection between much more different mixing formulas [Eqs. (6), (8), (9), and (13)] for low-filled (p ≪ pc) composites is practically impossible: the closer the mixing models, the higher the filling to reveal the difference. The reliable selection between Eq. (13) and Eq. (18) may be performed with a study of constitutive parameters of composites filled close to and above the structure inversion filling. Both the preparation and the permittivity measurements of composites overfilled with metal powders present experimental problems yet to be solved.

ARTICLE scitation.org/journal/adv
The dependence of permeability on frequency μ 1 (f ) for nickel inclusions [ Fig. 6(b)] is retrieved with Eqs. (18) and (16a) from the measured dependence of magnetic susceptibility on frequency for p = 0.375, N = 0.2, pc = 0.43, and δ = 1. The permeability curves retrieved with the proposed formula are close to those retrieved with Eq. (13), though the low-frequency permeability is about 20% higher than the cited data [the gray curves in Fig. 6(b)]. Permittivity measurements [ Fig. 6(a)] show that Eq. (18) at high filling is closer to the measured data than the Sihvola model; therefore, the microwave permeability of nickel spheres presented in Fig. 6(b) with black curves may be more accurate compared to the cited data. 15 The more accurate the mixing models (the more effects are taken into account), the more difficult the model selection from the measured data, and the wider filling range is necessary to be studied to distinguish the models reliably. The proposed mixing formula [Eq. (18)] predicts the dependence of the composite susceptibility on the filling and on the contrast of component properties close to that of the Sihvola model [Eq. (13)]; the difference reveals itself only at the fillings close to and above the percolation threshold p ≥ pc.

VII. CONCLUSIONS
The mixing formulas are usually applied to solve two mixing problems: the direct one is to calculate a constitutive parameter of a binary composite based on parameters of composite components and their volume fractions. The inverse problem is to retrieve a constitutive parameter of a mixture component from its volume fraction and from the measured constitutive parameters of the mixture. The direct mixing problem is usually an engineering task, while the inverse one arises if the direct measurements of the component susceptibility are impossible. The comparison of mixing models and the developed formula are proposed with the purpose to retrieve the microwave permeability of metals from the measured constitutive parameters of composites.
The analysis of mixing formulas (Tables I-III) for binary composites is based on the number of parameters fitted for the measured dependence of composite susceptibility on the volume fraction of components. The common approach to validate the applied formula by comparison with the measured data is incorrect, as the discrepancy depends on contrast of component properties and on the number of fitting parameters. The analysis points out the way to select the formulas accurately.
The widely used general Maxwell Garnett and symmetric Bruggeman models [Eqs. (6) and (8)] have one fitting parameter (namely, the inclusion shape-factor N) and do not account for the effect of mixture structure on the constitutive parameters. The models are suitable for rough estimation of inclusion parameters for diluted mixtures only. The reason is that the properties of a diluted mixture are defined more by the shape-factor than by the inclusion susceptibility χ 12 −1 ≥ (1 − p)N [see Eqs. (6) and (8)]: the lower 1/χ 12value compared to (1 − p)N, the lower the contribution of inclusion susceptibility χ 12 to the composite susceptibility χ 32 . The limitation may be amended by the increase in filling p, but there the above models fall beyond their accuracy range.
The more advanced Odelevskiy and Sihvola mixing models [Eqs. (12) and (13)] with two fitting parameters (the shape-factor N and the percolation threshold filling pc) are accurate within a wider filling range than Eqs. (6) and (8). The difference between the models reveals itself only for composites filled close to the percolation threshold pc.
The proposed model is a combination of two matrix models. We assume that the composite structure resembles a phase structure of the quenched or annealed binary system, depending on the mixing technology. The model accounts for structure transformations in a composite similar to the ones in a binary solid system; the difference is mainly in the scale of mixed components. The successful attempts to describe the variety of phase diagrams based on rigorous physical models are unknown, but all basic diagram types are listed in handbooks for physical chemistry. 31 The composite is considered as a system with limited solubility of both components: the solid solutions are similar to the direct and inverse matrix mixtures, whereas the heterogeneous mixture of the solutions is similar to a symmetric-type composite.
The direct matrix mixture is transformed in the framework of the proposed model into the inverse matrix mixture. The transformation depends on filling: the symmetric-type structure is considered as a transition between the direct and the inverse matrix structures. In contrast to the Odelevskiy approach, 11 the structure transformation takes place gradually within a transition filling range Δpc. The percolation threshold pc coincides with the structure inversion filling for Δpc → 0 only. The proposed model has three experimentally fitted parameters, namely, the shape-factor N, the structure inversion filling pc, and the filling fluctuation factor δ. The account for gradual transformation within the transition range Δpc(δ) increases the maximal filling where the formula is accurate compared to the two-parameter Odelevskiy and Sihvola models.
The proposed formula is tested by comparison with the measured dependence of complex permittivity and permeability on the frequency for a sample filled (p = 0.375) close to the percolation threshold with carbonyl nickel. The values of the shape-factor and critical filling are fitted for the measured dependence of permittivity on filling. 15 The measured hodograph of dielectric susceptibility appears to be closer to the proposed formula than to the Sihvola model (see Fig. 2), thus justifying the formula. At low fillings (p ≪ pc), the proposed matrix-inversion and the Sihvola formulas are difficult to distinguish (see Figs. 1-3).
The Bergman-Milton spectral function chart for the proposed model (Fig. 5) proves the percolation threshold behavior inherent to the effective-media models (Fig. 4), though the proposed model is a combination of matrix models.
The proposed model [Eq. (18)] assumes a Gaussian distribution of local filling around the inversion filling pc [Eq. (16a)]. Other distributions or simpler functions w(p) approximating the gradual inversion of matrix structure [e.g., Eq. (16b)] may be applied for Eq. (18) as well. The resulting difference in Bergman-Milton spectral functions and in mixture susceptibility for the accurate [Eq. (16a)] and approximated [Eq. (16b)] weight functions w(p) is problematic to distinguish experimentally. The problems with accurate measurements of both the permittivity and filling factor for the composites high-filled with metal powders make the refinement of model details or the increase in the number of fitting parameters barely expedient.
The comparison of shape-factor spectral function charts for the Sihvola (Fig. 4) and for the proposed models (Fig. 5) shows that the latter implies a more regular distribution of cluster shapes both below and above the critical filling pc. At the filling close to pc, the proposed model shows a spread of long clusters N → 0 including the ARTICLE scitation.org/journal/adv infinite ones, while the Sihvola model shows the infinite-long clusters only (the sharp delta-distribution peak at N = 0). The peak of cluster density with the shape factor higher than that of the initial inclusion at p ≫ pc may be attributed to orientation of long clusters perpendicular to the exciting field. The symmetric-type models show at p ≫ pc, a wide and low maximum of the shape-factors partially exceeding that of the initial inclusion. The difference of shape-factor spectral functions becomes obvious at p → 1.
The shape-distribution for the proposed model (Fig. 5) seems more realistic than the distribution for the symmetric-type models (Fig. 4).
The more accurate the mixing model accounts for inclusion interaction, the more general the model is and the closer the retrieved parameters of mixture components are to their real properties. This interaction cannot be calculated rigorously, e.g., for multipole interaction, because the inclusion arrangement is defined by the structure transformations that the composite undergoes with the increase in filling p. In the case of the Odelevskiy and Sihvola models, the effect of structure transformation is accounted for by one fitting parameter related to the threshold filling (pc). The proposed model accounts for transformation of the mixture structures with two fitting parameters (pc and δ), and the second parameter accounts for fluctuations in inclusion arrangement and is related to the transition filling range Δpc(δ). The model is more general than the published ones (Tables I-III) (14)] and asymmetric, or differential [Eqs. (10) and (11)] type models.
The account for transformations of composite structures from direct matrix to symmetric and to inverse matrix types makes the proposed model more general and accurate up to higher filling p compared to the known mixing models. Therefore, the model application to retrieve the inclusion permeability μ 1 from the measured constitutive parameters of a composite may amend the limitation (1 − p)N ≤ (μ 1 − 1) −1 on the retrieved susceptibility and increase the accuracy at low frequencies where the inclusion permeability μ 1 may be high due to the frequency dispersion.