Geometrical Description in Binary Composites and Spectral Density Representation

Abstract

In this review, the dielectric permittivity of dielectric mixtures is discussed in view of the spectral density representation method. A distinct representation is derived for predicting the dielectric properties, permittivities ε, of mixtures. The presentation of the dielectric properties is based on a scaled permittivity approach, $\xi =({\epsilon }_{\mathsf{e}}-{\epsilon }_{\mathsf{m}}){({\epsilon }_{\mathsf{i}}-{\epsilon }_{\mathsf{m}})}^{-1}$, where the subscripts e, m and i denote the dielectric permittivities of the effective, matrix and inclusion media, respectively [Tuncer, E. J. Phys.: Condens. Matter 2005, 17, L125]. This novel representation transforms the spectral density formalism to a form similar to the distribution of relaxation times method of dielectric relaxation. Consequently, I propose that any dielectric relaxation formula, i.e., the Havriliak-Negami empirical dielectric relaxation expression, can be adopted as a scaled permittivity. The presented scaled permittivity representation has potential to be improved and implemented into the existing data analyzing routines for dielectric relaxation; however, the information to extract would be the topological/morphological description in mixtures. To arrive at the description, one needs to know the dielectric properties of the constituents and the composite prior to the spectral analysis. To illustrate the strength of the representation and confirm the proposed hypothesis, the Landau-Lifshitz/Looyenga (LLL) [Looyenga, H. Physica 1965, 31, 401] expression is selected. The structural information of a mixture obeying LLL is extracted for different volume fractions of phases. Both an in-house computational tool based on the Monte Carlo method to solve inverse integral transforms and the proposed empirical scaled permittivity expression are employed to estimate the spectral density function of the LLL expression. The estimated spectral functions for mixtures with different inclusion concentration compositions show similarities; they are composed of a couple of bell-shaped distributions, with coinciding peak locations but different heights. It is speculated that the coincidence in the peak locations is an absolute illustration of the self-similar fractal nature of the mixture topology (structure) created with the LLL expression. Consequently, the spectra are not altered significantly with increased filler concentration level—they exhibit a self-similar spectral density function for different concentration levels. Last but not least, the estimated percolation strengths also confirm the fractal nature of the systems characterized by the LLL mixture expression. It is concluded that the LLL expression is suitable for complex composite systems that have hierarchical order in their structure. These observations confirm the finding in the literature.

1. Introduction

Electrical properties of composite materials have attracted researchers to seek relationships between overall composite properties and intrinsic properties of the parts forming the mixture (constituents) and their spatial arrangement in the mixture [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. Mixture formulas based on analytical and effective medium approaches were developed such that for various arrangements of inclusions, predicting the dielectric properties of composites was plausible [15,16,17,18].

A deep understanding of dielectric mixtures would be of great value for (i) calculating the dielectric constant of a mixture composed of substances with known dielectric constants, (ii) calculating the dielectric constant of the second component of a two-component mixture when the dielectric constants of the mixture and the first component are known [1], and (iii) estimating the morphology of a two-component mixture when the dielectric constants of the mixture and each of the components are known [19,20,21].

During the 1970s Fuchs [22,23] illustrated that the dielectric permittivity of mixtures could be expressed as a summation of depolarization factors. This approach took into account that each particle would yield a different depolarization factor depending on its shape. The beauty was that the summation term would yield the concentration of the individual shapes. In the late 1970s, Bergman showed that one can separate the geometrical contributions from the pure dielectric response of a composite if and only if the dielectric properties of the constituents are known [24,25,26]. Milton corrected errors in Bergman’s original derivation [27,28,29], and later Golden and Papanicolaou [30,31] gave the rigorous derivation for the spectral representation theory. If the Fuchs and Bergman approaches are studied in detail, both correspond to the same formalism, which is spectral density representation [19,20,23,24,32,33,34].

Notice that the present author has illustrated similarities between the dielectric relaxation and dielectric response of dielectric mixtures using the spectral density representation; the origin of similarities is very significant for the comprehension of the physics of dielectrics [21].

The concept of having cognition about the structure of composites, how the phases are arranged, is very useful in materials design and characterization, because special materials can be manufactured with a knowledge of structure-property relationships. We should not forget that the functionality in nature is created with structure. For regular arrangements of phases, there exist equations based on theoretical calculations of simple enough geometries [15,18]. However, disordered and fractal structures are abundant in nature [35,36], and a comprehension of materials properties with disordered or fractal structures has been a challenge for researchers for some decades. The fractal geometry or systems indicating hierarchical order has been one of the interesting topics in applied and theoretical (mathematical) physics [35,37,38,39,40]. As an example, the electrical properties of metal aggregates in insulating matrix media were studied extensively [36,41,42,43,44,45,46,47,48,49]. In these studies either the dimension of the electrical network or system, or the resonance frequency of the electrical impedance was used as a measure to indicate the fractal dimensions.

When real systems are taken into account, the structural information is, on the other hand, usually obtained by optical/microscopic techniques, this information is later analyzed to estimate the fractal dimensions. While if one applies the spectral density approximation to electrical data of a mixture, structural information could be obtained—whether the mixture has regular, disordered or fractal structural characteristics or not. We should consider the application of the spectral density to the mixture immittance (Immittance term is used to express both impedance and admittance in a combined from.) data as optical or X-ray techniques such that structural characteristics of the mixture are determined with the analysis—the resolved depolarization factors corresponds to the arrangements of phases [19,20,49].

One can also utilize a dielectric mixture formula to model the electrical properties of the composite system in hand, such that the model contains structural information, e.g., there exist effective medium theories for composites with spherical and ellipsoidal inclusions [3,50,51]. In this review, we employ the spectral density representation as the general representation approach for composites to resolve the geometrical description of a model system described by the Landau-Lifshitz/Looyenga (LLL) effective medium formula [52,53], or in other words we challenge the physical significance of the LLL expression.

The LLL expression was extensively used to describe the dielectric properties of dispersive systems composed of powders or exhibiting porous structure [3,54,55,56,57,58,59,60,61,62,63,64,65,66,67]. It was even shown [55] that the LLL formula was more reliable when mixtures contained strongly dissipative particles and was compared to others like Maxwell Garnett (MG) [68,69], Bruggeman [70], etc. (see for example Refs. [1,2,3,4,5,6,7,8,9,10,11,12] for other formulas). First, it is presented that the spectral density representation can in fact be written in a novel, elegant form that can be implemented in already existing dielectric data analysis techniques [71,72,73,74,75]. Later, the simplified notation is used with a numerical procedure which is derived to solve inverse problems [19,20,76,77] to resolve the spectral densities in the LLL [52,53] dielectric mixture expression. The significance of the presented approach on dielectric mixtures is also discussed.

The paper is organized as follows. The spectral density representation for a binary mixture is presented in §2. In addition, similarities between dielectric relaxation in dielectrics and dielectric permittivity of binary mixtures are illustrated. The dielectric data representation is described in §3., and hints are given for analyzing impedance data of mixtures. The numerical method for solving the inverse integral is also presented explicitly for interested readers in §4. The numerical data generation and the LLL expression are presented in §5. The results obtained by the inverse integral solution and by the proposed conventional dielectric dispersion expression are compared in §6. The conclusions are in contained §7.

2. Spectral Density Representation

In spectral density representation analysis of binary mixtures, the dielectric permittivity of a heterogeneous (effective) medium is expressed as [24,25,26,27,28,29,30,31,33,78,79]

$$\begin{array}{c}\hfill {\epsilon }_{\mathsf{e}}={\epsilon }_{\mathsf{m}}\left\{1+q\left[A\left(\frac{{\epsilon }_{\mathsf{i}}}{{\epsilon }_{\mathsf{m}}}-1\right)+{\int }_0^1\frac{\mathsf{G}\left(x\right)\mathrm{d}x}{{\left(\frac{{\displaystyle {\epsilon }_{\mathsf{i}}}}{{\displaystyle {\epsilon }_{\mathsf{m}}}}-1\right)}^{-1}+x}\right]\right\}\end{array}$$

(1)

where ${\epsilon }_{\mathsf{e}}$, ${\epsilon }_{\mathsf{m}}$ and ${\epsilon }_{\mathsf{i}}$ are the complex dielectric permittivity of the effective, matrix and inclusion media, respectively; q and x are the concentration of inclusions and the spectral parameter, respectively. The function $\mathsf{G}\left(x\right)$ is the spectral density function (SDF) and possesses information about the topological description of the mixture. Eq. (1) can be arranged in a simplified form, cf. § 7, as follows:

$$\begin{array}{c}\hfill \xi ={\xi }_s+q{\int }_0^1\frac{\mathsf{G}\left(x\right)\mathrm{d}x}{1+{\epsilon }_{\mathsf{m}}^{-1}{\mathsf{\Delta }}_{\mathsf{im}}x}\end{array}$$

(2)

where ${\mathsf{\Delta }}_{\mathsf{im}}={\epsilon }_{\mathsf{i}}-{\epsilon }_{\mathsf{m}}$ and ξ is the complex and frequency dependent “scaled” permittivity [21]

$$\begin{array}{c}\hfill \xi =\frac{{\epsilon }_{\mathsf{e}}-{\epsilon }_{\mathsf{m}}}{{\epsilon }_{\mathsf{i}}-{\epsilon }_{\mathsf{m}}}\end{array}$$

(3)

The constant ${\xi }_s$ in Equation 2 is complex and depends on the concentration and structure of the composite; its real part is related to the so-called “percolation strength” [33,80]. The mathematical properties and conditions that SDF satisfies are presented in §7. [26,33,78,79,81].

Equation (3) is a very similar expression to the distribution of relaxation times (DRT) representation of a broad dielectric dispersion (relaxation) [21,76,82,83,84,85,,86,87,88,89]

$$\begin{array}{c}\hfill \epsilon \left(ı\omega \right)={\epsilon }_{\infty }+\mathsf{\Delta }\epsilon {\int }_0^{\infty }\frac{\mathsf{G}\left(\tau \right)\mathrm{d}\tau }{1+ı\omega \tau }\end{array}$$

(4)

where ε, ${\epsilon }_{\infty }$ and $\mathsf{\Delta }\epsilon$ are the complex permittivity, permittivity at optical frequencies and dielectric strength, respectively; and $ı\equiv \sqrt{-1}$. The quantities ω and τ are the angular frequency and relaxation time, respectively. The distribution function for the relaxation times is $\mathsf{G}\left(\tau \right)$. Comparison of Equations (2) and (4) demonstrates that both the DRT and the scaled permittivity of SDF are actually the same. However, the new complex parameter ϖ in SDF $(\varpi \equiv {\epsilon }_{\mathsf{m}}^{-1}{\mathsf{\Delta }}_{\mathsf{im}})$ is similar to the pure complex frequency $ı\omega$ in DRT representation, and the real number constant ${\epsilon }_{\infty }$ is a complex number ${\xi }_s$ in SDF representation.

The significance of ${\xi }_s$ is that, just like ${\epsilon }_{\infty }$, it does not indicate any dispersion in the dielectric spectra, meaning that it does not contribute to the dielectric relaxation. One can also look at these two parameters as regions with no charge or dipolar movements. In addition, the spectral parameter x is similar to the relaxation time τ in the DRT. Finally, the dielectric strength $\mathsf{\Delta }\epsilon$ in the DRT representation is like the concentration of inclusions in the SDF.

Because of these similarities between the SDR and DRT approaches, methods developed for dielectric data analysis [71,73,74,75] can be applied to the scaled permittivity ξ of SDF [19,20,21,34,49]. For example, one of the most employed dielectric dispersion expressions, the Havriliak-Negami empirical expression [90], can be used to analyze the scaled complex dielectric permittivity data of a mixture [34],

$$\begin{array}{c}\hfill \xi \left(\varpi \right)={\xi }_s+\frac{q}{{[1+{\left(\varpi x\right)}^{\alpha }]}^{\beta }}\end{array}$$

(5)

where α and β are parameters of a general distribution function [73,90], and ϖ is the scaled complex frequency. When $\alpha =\beta =1$, a Debye-type relaxation is observed in the dielectric dispersion representation [91]. In the case of spectral density representation, however, the Maxwell Garnett approximation is obtained for $\alpha =\beta =1$ and $x=(1-q)/d$, since $\mathsf{G}\left(x\right)=\delta [x-(1-q)/d]$ [19,20,21,78,79,81], where d is the dimensionality of the system.

One can, therefore, in principle write a new, more general empirical mixture formula [34] by isolating the dielectric permittivity ${\epsilon }_{\mathsf{e}}$ of the composite in Equation (3) and substituting it in Equation (5),

$$\begin{array}{c}\hfill {\epsilon }_{\mathsf{e}}={\epsilon }_{\mathsf{m}}+{\mathsf{\Delta }}_{\mathsf{im}}\left\{{\xi }_s+\frac{q}{{[1+{\left(\varpi x\right)}^{\alpha }]}^{\beta }}\right\}\end{array}$$

(6)

Finally, note that the fractional expression inside the curly parentheses in Equation (6) can be exchanged with any of the dielectric dispersion relations in the literature [71,74,90,92,93,94,95,96].

3. Representation of Dielectric Data

The dielectric function for a d-dimensional (or composite with arbitrarily shaped inclusions) is defined as follows with Maxwell Garnett (MG) expression for a composite [68,69]

$$\begin{array}{c}\hfill {\epsilon }_{\mathsf{e}}^{\mathrm{MG}}(\omega ;{\epsilon }_{\mathsf{m}},{\epsilon }_{\mathsf{i}},q,d)={\epsilon }_{\mathsf{m}}+\frac{{\epsilon }_{\mathsf{m}}dq{\mathsf{\Delta }}_{\mathsf{im}}}{(1-q){\mathsf{\Delta }}_{\mathsf{im}}+d{\epsilon }_{\mathsf{m}}}\end{array}$$

(7)

The dielectric data of the composite can be expressed in one of the four immittance representations [72,74,83]: (i) the complex resistivity $\rho \left(\omega \right)$, (ii) the complex modulus $M\left(\omega \right)\equiv ı\omega {\epsilon }_0\rho \left(\omega \right)$; (iii) the complex permittivity $\epsilon \equiv {\left[M\left(\omega \right)\right]}^{-1}$, and (iv) the complex conductivity $\sigma \left(\omega \right)\equiv ı\omega {\epsilon }_0\epsilon \left(\omega \right)\equiv {\left[\rho \left(\omega \right)\right]}^{-1}$. When we are dealing with frequency dependent dielectric properties of composites, the effective conductivity of the composite can sometimes influence the imaginary part of the dielectric function—hindering the dielectric losses due to the interfacial polarization—as shown in inset (b) of Figure 1. In such cases it is more appropriate to use the complex resistivity representation (plot) as shown in the Argand diagrams in log-log and linear scales in Figure 1 and 1c, respectively.

Figure 1. The Argand diagram of resistivity for a Maxwell Garnett composite in log-log scale, the effective permittivity is calculated with ${\epsilon }_{\mathsf{e}}={\epsilon }_{\mathsf{e}}^{\mathrm{MG}}(\omega ;{\epsilon }_{\mathsf{m}},{\epsilon }_{\mathsf{i}},0.3,3)$, with ${\epsilon }_{\mathsf{m}}=2+{10}^{-12}{\left(ı{\epsilon }_0\omega \right)}^{-1}$ and ${\epsilon }_{\mathsf{i}}=10+{10}^{-10}{\left(ı{\epsilon }_0\omega \right)}^{-1}$. The arrows indicate the direction of increasing angular frequency ω. Insets (a) and (b) are the real and the imaginary parts of permittivities; ${\epsilon }_{\mathsf{e}}$, ${\epsilon }_{\mathsf{m}}$ and ${\epsilon }_{\mathsf{i}}$ are the permittivities for the matrix (– – –), inclusion (- - - - -) and effective (——–) media. Inset (c) is the Argand diagram of resistivity in linear-linear scale. Logarithmic scale is base 10.

Figure 1. The Argand diagram of resistivity for a Maxwell Garnett composite in log-log scale, the effective permittivity is calculated with ${\epsilon }_{\mathsf{e}}={\epsilon }_{\mathsf{e}}^{\mathrm{MG}}(\omega ;{\epsilon }_{\mathsf{m}},{\epsilon }_{\mathsf{i}},0.3,3)$, with ${\epsilon }_{\mathsf{m}}=2+{10}^{-12}{\left(ı{\epsilon }_0\omega \right)}^{-1}$ and ${\epsilon }_{\mathsf{i}}=10+{10}^{-10}{\left(ı{\epsilon }_0\omega \right)}^{-1}$. The arrows indicate the direction of increasing angular frequency ω. Insets (a) and (b) are the real and the imaginary parts of permittivities; ${\epsilon }_{\mathsf{e}}$, ${\epsilon }_{\mathsf{m}}$ and ${\epsilon }_{\mathsf{i}}$ are the permittivities for the matrix (– – –), inclusion (- - - - -) and effective (——–) media. Inset (c) is the Argand diagram of resistivity in linear-linear scale. Logarithmic scale is base 10.

The ohmic conductivity ${\sigma }_{\mathsf{e},dc}$ (or resistivity ${\rho }_{\mathsf{e},dc}$) of the material can be estimated from the complex resistivity Argand plot. However, once the conductivity contributions are cleared from the immittance or dielectric data, using the estimated resistivity value in Figure 1 as $\omega \to 0$, the pure dielectric dispersion (permittivity) would be obtained. In addition, the high frequency dielectric permittivity $\epsilon (\omega \to \infty )\equiv {\epsilon }_{\infty }$ can be further subtracted from the data to obtain the pure dielectric polarization (susceptibility χ) of the composite as presented in Figure 2; $\chi ={\chi }^{\prime }-ı{\chi }^{\prime \prime }\equiv \epsilon -{\epsilon }_{\infty }+ı{\sigma }_{dc}{\left({\epsilon }_0\omega \right)}^{-1}$. The imaginary part of χ has a peak around $\omega \sim 1^{+}\mathrm{rad}{\mathrm{s}}^{-1}$. The linear scale plot of the susceptibility χ is a semicircular curve, cf. Figure 2b, as seen in Figure 1c.

If we now consider the frequency dependent properties of the scaled permittivity ξ for the considered MG composite, the real part and the imaginary parts are similar to the dielectric permittivity ${\epsilon }_{\mathsf{e}}$ of the composite, cf. inset in Figure 3, which shows the Argand diagram of the scaled permittivity ξ; observe that the increasing frequency is in the opposite direction when compared to the Argand plot of the susceptibility in Figure 2b. The real part of ξ is a mirror image of $\Re \left({\epsilon }_{\mathsf{e}}\right)$. Unlike the imaginary part of ${\epsilon }_{\mathsf{e}}$ (due to ohmic losses), $\Im \left(\xi \right)$ shows a clear peak around $\omega \sim 1^{-}\mathrm{rad}{\mathrm{s}}^{-1}$. The shift in the origin position in the Argand diagram in the inset of Figure 3 is related to the percolation strength ${\xi }_s$, which is close to zero, and the concentration of the inclusions.

Figure 2. (a) Dielectric susceptibility χ as a function of angular frequency ω; the real and the imaginary parts are presented with solid (——–) and dashed (– – –) lines, respectively. The arrow indicates the direction of increasing angular frequency ω. (b) The Argand diagram (Cole-Cole plot) of susceptibility. Logarithmic scale is base 10.

Figure 2. (a) Dielectric susceptibility χ as a function of angular frequency ω; the real and the imaginary parts are presented with solid (——–) and dashed (– – –) lines, respectively. The arrow indicates the direction of increasing angular frequency ω. (b) The Argand diagram (Cole-Cole plot) of susceptibility. Logarithmic scale is base 10.

There are similarities between the susceptibility and the scaled permittivity plots when the same frequency ω axis is used; however, note that the scaled frequency ϖ for the scaled permittivity is a complex quantity. We therefore illustrate the dependence of the real angular frequency ω as a function of ϖ in a three dimensional curve-plot in Figure 4. In addition, the real and the imaginary parts of the scaled permittivity ξ are shown in Figure 5 as a function of ϖ. As shown in the figure, the actual dependence of ξ on ϖ is more complicated than ε on ω. On the contrary, this dependence can be used to estimate the spectral density function for a given system, which is explicitly given with a numerical procedure in the next section.

Figure 3. The real ($\Re \left(\xi \right)$; ——–) and the imaginary ($\Im \left(\xi \right)$; – – –) parts of scaled permittivity ξ in the Maxwell Garnett approximation. Inset is the Argand diagram for ξ with the arrow showing the direction of increasing angular frequency ω. The inclusions are spherical $d=3$, and $q=0.3$. Logarithmic scale is base 10.

Figure 3. The real ($\Re \left(\xi \right)$; ——–) and the imaginary ($\Im \left(\xi \right)$; – – –) parts of scaled permittivity ξ in the Maxwell Garnett approximation. Inset is the Argand diagram for ξ with the arrow showing the direction of increasing angular frequency ω. The inclusions are spherical $d=3$, and $q=0.3$. Logarithmic scale is base 10.

Figure 4. Dependence of the real frequency ω as a function of complex “scaled” frequency ϖ in the spectral representation. Logarithmic scale is base 10.

Figure 4. Dependence of the real frequency ω as a function of complex “scaled” frequency ϖ in the spectral representation. Logarithmic scale is base 10.

Figure 5. Three-dimensional line plots of (a) the real $\Re \left(\xi \right)$ and (b) the imaginary $\Im \left(\xi \right)$ parts of the “scaled” permittivity for a Maxwell Garnett mixture. The inclusions are spherical $d=3$, and $q=0.3$.

Figure 5. Three-dimensional line plots of (a) the real $\Re \left(\xi \right)$ and (b) the imaginary $\Im \left(\xi \right)$ parts of the “scaled” permittivity for a Maxwell Garnett mixture. The inclusions are spherical $d=3$, and $q=0.3$.

4. Numerical Estimation of Spectral Density Function

The derived spectral density expression in Equation 2 is a Bolter equation [97], which is a special form of the Fredholm integral equations [98]. Such equations are usually considered to be ill-conditioned because of their non-unique solutions. However, the approach used here and presented several times elsewhere [19,73,76,77,89] leads to unique solutions by means of a constrained least-squares fit and the Monte Carlo integration methods. Some other approaches to solving the spectral density function are also suggested in the literature [33,78,79,80,81,99,100,101,102]. The presented numerical method to solve inverse integral transforms has been previously previously been used in different problems [19,20,72,73,76,77,89]. In this particular approach, the integral in Equation (3) is first written in a summation form over some number of randomly selected and fixed $x_n$-values, $x_n\in [0,1]$, where n is less than the total number M of experimental (known) data points in the complex scaled permittivity, ξ,

$$\begin{array}{c}\hfill \xi ={\xi }_s+\sum \frac{{\mathsf{g}}_n}{1+\varpi x_n}n\le M\end{array}$$

(8)

This converts the non-linear problem in hand to a linear one with ${\mathsf{g}}_n$ being the unknowns, weights of the randomly selected $x_n$ values. In the present notation $\mathsf{g}$ is $q\mathsf{G}$. Later, a constrained least-squares algorithm is applied to get the corresponding $\mathsf{g}$-values and ${\xi }_s$,

$$\begin{array}{c}\hfill min\sum {\left[\underline{\xi }-\mathfrak{K}\mathfrak{g}\right]}^2\mathrm{and}\mathfrak{g}\ge 0\end{array}$$

(9)

where $\mathfrak{K}$ is the kernel-matrix,

$$\begin{array}{c}\hfill \mathfrak{K}=\left(\begin{array}{cccc}1& {\mathbf{K}}_{11}& {\mathbf{K}}_{12}& \dots \\ 1& {\mathbf{K}}_{21}& {\mathbf{K}}_{22}& \dots \\ 1& {\mathbf{K}}_{31}& {\mathbf{K}}_{32}& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right)\end{array}$$

(10)

Here ${\mathbf{K}}_{ij}={[1+{\varpi }_i{x}_j]}^{-1}$; index i runs on the angular frequency points $i=1,\cdots ,M$; and index j runs on the randomly selected x values, $j=1,\cdots ,n$. The parameters $\underline{\xi }$ and $\mathfrak{g}$ in Equation (9) are column vectors, respectively, the scaled permittivity calculated from the experimental (known) data and the searched spectral density,

$$\begin{array}{ccc}\hfill \underline{\xi }=\left(\begin{array}{c}\xi \left[\varpi \left({\omega }_1\right)\right]\\ \xi \left[\varpi \left({\omega }_2\right)\right]\\ \xi \left[\varpi \left({\omega }_3\right)\right]\\ \xi \left[\varpi \left({\omega }_4\right)\right]\\ ⋮\end{array}\right)& \mathrm{and}& \mathfrak{g}=\left(\begin{array}{c}{\xi }_s\\ {\mathsf{g}}_1\\ {\mathsf{g}}_2\\ {\mathsf{g}}_3\\ ⋮\end{array}\right)\hfill \end{array}$$

(11)

In our numerical procedure, we perform many minimization steps with fresh, new sets of randomly selected $x_j$-values. Here we have adopted $2^{12}$ minimization steps. The ${\mathsf{g}}_j$-values and ${\xi }_s$ obtained are recorded in each step, which later build up the spectral density distribution $\mathsf{g}$ and a distribution for the percolation strength, ${\xi }_s$. For a large number of minimization loops, the x-axis actually becomes continuous—the Monte Carlo integration hypothesis—contrary to regularization methods [80,101]. The number M of data points chosen is 24, and the number N of unknown $\mathsf{g}$-values is 22. In the analysis presented below, the total number of randomly selected x values is $2^{12}\times 22=90112$.

Application of the numerical procedure to the Maxwell Garnett expression, impedance data of a porous rock-brine mixture and two-dimensional “ideal” structures, and the accuracy in the numerical inversion, have previously been presented elsewhere [19,20]. The estimated spectral density functions for the MG expression were delta sequences [103] as expected, without any significant percolation component, because of the estimated concentration q, cf. Equation (20). In the next section, we apply the numerical procedure to the Landau-Lifshitz/Looyenga [52,53] expression to better understand the nature of dielectric mixtures, which obey this relation.

5. Landau-Lifshitz/Looyenga Expression

Landau and Lifshitz [52] and Looyenga [53] independently, using different approaches, developed an expression for dielectric mixtures, that implies the additivity of cube roots of the permittivities of mixture constituents when taken in proportion to their volume fractions (see Refs. [3,54,55,56,57,58,59,60,61,62,63,64,65,66,67] for examples).

$$\begin{array}{c}\hfill {\left[{\epsilon }_{\mathsf{e}}^{\mathrm{LLL}}(\omega ;{\epsilon }_{\mathsf{m}},{\epsilon }_{\mathsf{i}},q)\right]}^{1/3}=(1-q){\epsilon }_{\mathsf{m}}^{1/3}+q{\epsilon }_{\mathsf{i}}^{1/3}\end{array}$$

(12)

This expression is used extensively in the literature for powdered materials and optical properties of material mixtures [3,56]. In the following calculations, we choose the same values for the dielectric functions of the phases as before, cf. Figure 1. The concentration q of the inclusion phase is varied between 0.1 and 0.9 in the simulations.

The extracted spectral functions are like distributions, and they are analyzed by means of comparing them with a known distribution. We apply the Lévy statistics [104,105,106,107], which are widely used for interacting systems in different research fields [19,89,106,107,108,109,110,111,112]. The Lévy stable distribution is a natural generalization (approximation) of the normal (Gaussian), Cauchy or Lorenz and Gamma distributions. It is used when analyzing sums of independent, identically distributed, random variables by a diverging variance. Its characteristic function is expressed as

$$\begin{array}{c}\hfill \mathsf{L}(x;\mathfrak{A},\mu ,\gamma ,\zeta )=\mathfrak{A}|exp\{-\left|\zeta (x-\mu ){|}^{\gamma }\right\}|\end{array}$$

(13)

Here, γ is the characteristic exponent ($\gamma >0$), μ is the localization parameter, ζ is the scale parameter and $\mathfrak{A}$ is the amplitude. The special forms of Equation (13) are the Gaussian [$\mathsf{L}(x;A,\mu ,2,\zeta )$], the Lorentz or Cauchy [$\mathsf{L}(x;A,\mu ,1,\zeta )$] and Gamma [$\mathsf{L}(x;A,\mu ,1/2,\zeta )$] distributions. Different forms of probability density functions for Lévy statistics exist, and we adopted a stable distribution used in the literature [104,105,106]. We omitted the imaginary parts in the characteristic function because of their insignificance in the results.

Figure 6. The Argand diagram of resistivity for a Landau-Lifshitz/Looyenga composite in log-log scale; the effective permittivity is calculated by ${\epsilon }_{\mathsf{e}}={\epsilon }_{\mathsf{e}}^{\mathrm{LLL}}(\omega ;{\epsilon }_{\mathsf{m}},{\epsilon }_{\mathsf{i}},0.3,3)$, with ${\epsilon }_{\mathsf{m}}=2+{10}^{-12}{\left(ı{\epsilon }_0\omega \right)}^{-1}$ and ${\epsilon }_{\mathsf{i}}=10+{10}^{-10}{\left(ı{\epsilon }_0\omega \right)}^{-1}$. The arrows indicate the direction of increasing angular frequency ω. Insets (a) and (b) are the real and the imaginary parts of permittivities; ${\epsilon }_{\mathsf{e}}$, ${\epsilon }_{\mathsf{m}}$ and ${\epsilon }_{\mathsf{i}}$ are the permittivities for the matrix (– – –), inclusion (- - - - -) and effective (——–) media. Inset (c) is the Argand diagram of resistivity in linear-linear scale. Logarithmic scale is base 10.

Figure 6. The Argand diagram of resistivity for a Landau-Lifshitz/Looyenga composite in log-log scale; the effective permittivity is calculated by ${\epsilon }_{\mathsf{e}}={\epsilon }_{\mathsf{e}}^{\mathrm{LLL}}(\omega ;{\epsilon }_{\mathsf{m}},{\epsilon }_{\mathsf{i}},0.3,3)$, with ${\epsilon }_{\mathsf{m}}=2+{10}^{-12}{\left(ı{\epsilon }_0\omega \right)}^{-1}$ and ${\epsilon }_{\mathsf{i}}=10+{10}^{-10}{\left(ı{\epsilon }_0\omega \right)}^{-1}$. The arrows indicate the direction of increasing angular frequency ω. Insets (a) and (b) are the real and the imaginary parts of permittivities; ${\epsilon }_{\mathsf{e}}$, ${\epsilon }_{\mathsf{m}}$ and ${\epsilon }_{\mathsf{i}}$ are the permittivities for the matrix (– – –), inclusion (- - - - -) and effective (——–) media. Inset (c) is the Argand diagram of resistivity in linear-linear scale. Logarithmic scale is base 10.

In Figure 6, the simulated dielectric permittivity with Equation (12) is presented at the complex resistivity level. The actual dielectric data is shown in the insets in Figure 6a,b. Compared to the MG expression in Equation 7, shown in Figure 1, the LLL expression does not show a knee-point as the complex resistivity decreases with increasing frequency. In addition, the linear scale Argand diagram of the LLL expression in Figure 6c is not a perfect semicircle as the MG one. The dielectric susceptibility after subtraction of the frequency independent parameters ${\epsilon }_{\infty }$ and σ is presented in Figure 7. The losses, ${\chi }^{\prime \prime }$, are non-symmetrical for the LLL expression, which indicates that the actual dielectric response can be modeled by a non-Debye dielectric dispersion, e.g., the Havriliak-Negami expression; this is also visible in the Argand diagram of the susceptibility, cf. Figure 7b.

Figure 7. (a) Dielectric susceptibility χ as a function of angular frequency ω; the real and the imaginary parts are presented with the solid (——–) and dashed (– – –) lines, respectively. The arrow indicates the direction of increasing angular frequency ω. (b) The Argand diagram (Cole-Cole plot) of susceptibility. Logarithmic scale is base 10.

Figure 7. (a) Dielectric susceptibility χ as a function of angular frequency ω; the real and the imaginary parts are presented with the solid (——–) and dashed (– – –) lines, respectively. The arrow indicates the direction of increasing angular frequency ω. (b) The Argand diagram (Cole-Cole plot) of susceptibility. Logarithmic scale is base 10.

The reconstructed complex dielectric permittivity and the scaled permittivity are presented in Figure 8. It is important to mention that the fitting is performed at the scaled permittivity level. While numerically calculating the spectral function, the randomly selected spectral parameters x are picked between ${10}^{-3}$ and 1 in logarithmic scale; the pre-distribution is log-linear (for details see [76]). In such an integration limit as in Equation (2), the distribution or the functional contributions of x-values lower than ${10}^{-3}$ are included in the percolation strength ${\xi }_s$.

Figure 8. (a) The real and (b) the imaginary parts of the dielectric permittivity calculated with Equation (12) for $q=0.3$; the permittivities of the matrix and the inclusion phases are ${\epsilon }_{\mathsf{m}}=2+{10}^{-12}{\left(ı{\epsilon }_0\omega \right)}^{-1}$ and ${\epsilon }_{\mathsf{i}}=10+{10}^{-10}{\left(ı{\epsilon }_0\omega \right)}^{-1}$. (c) Argand diagram of scaled permittivity without the percolation strength contribution. The inset in (a) shows the dielectric susceptibility after the subtraction of ohmic conductivity and the permittivity at high frequencies. The inset in (c) shows the Argand plot of ${\xi }_s$ for each point; the values are very narrowly distributed. The lines with points are the simulated data of Equation 12, and the symbols (∘) are the dielectric response calculated with the estimated spectral density function from the proposed numerical algorithm. There is very good agreement between the simulated and the analyzed data sets. Logarithmic scale is base 10.

Figure 8. (a) The real and (b) the imaginary parts of the dielectric permittivity calculated with Equation (12) for $q=0.3$; the permittivities of the matrix and the inclusion phases are ${\epsilon }_{\mathsf{m}}=2+{10}^{-12}{\left(ı{\epsilon }_0\omega \right)}^{-1}$ and ${\epsilon }_{\mathsf{i}}=10+{10}^{-10}{\left(ı{\epsilon }_0\omega \right)}^{-1}$. (c) Argand diagram of scaled permittivity without the percolation strength contribution. The inset in (a) shows the dielectric susceptibility after the subtraction of ohmic conductivity and the permittivity at high frequencies. The inset in (c) shows the Argand plot of ${\xi }_s$ for each point; the values are very narrowly distributed. The lines with points are the simulated data of Equation 12, and the symbols (∘) are the dielectric response calculated with the estimated spectral density function from the proposed numerical algorithm. There is very good agreement between the simulated and the analyzed data sets. Logarithmic scale is base 10.

The estimated spectral density functions for two concentration levels, $q=0.2$ and $q=0.3$, are presented in Figure 9. There are six visible peaks, labeled one to six from left to right. (We do not take peak 0 into account in the analysis.) The estimated integrals of portion of the bell-shaped distributions are presented on the graphs. Each peak is analyzed by the Lévy distribution as mentioned before, and the solid lines show the individual distributions, cf. Figure 9. It is remarkable that the estimated distributions are similar in form but shifted up a bit with increase in concentration q. This behavior is an illustration of the self-similar fractal nature of the considered composite system in the LLL expression, in which the topological arrangement does not change significantly with increased concentration; see shifts in SDF of regular lattices with increasing concentrations in Ref. [20]. To support this statement, the spectral density functions of mixtures described with the LLL expression for nine different concentrations are shown in Figure 10. The data are shifted with constant steps in concentration for clarity; however, as observed the spectral function amplitudes do not increase proportionally with increasing inclusion concentration, cf. Table 1 and cf. Figure 11. The difference between peak positions is constant in logarithmic scale $|P_j-P_{j-1}|\approx 0.47$, cf. Figure 12.

Figure 9. Spectral density functions of the Landau-Lifshitz/Looyenga equation for two concentrations; (a) $q=0.2$ and (b) $q=0.3$. The solid lines (——–) are Lévy distributions adapted to the estimated peaks. The numerical data is presented with symbol, error bars. The numbers on the peaks denote the significant peaks. Logarithmic scale is base 10.

Figure 9. Spectral density functions of the Landau-Lifshitz/Looyenga equation for two concentrations; (a) $q=0.2$ and (b) $q=0.3$. The solid lines (——–) are Lévy distributions adapted to the estimated peaks. The numerical data is presented with symbol, error bars. The numbers on the peaks denote the significant peaks. Logarithmic scale is base 10.

Table 1 lists not only the positions but also the shape parameters of the Lévy expression for the six peaks resolved for each concentration. The peak positions μ are not varied with increasing concentration of the inclusions phase. The amplitude $\mathfrak{A}$ changes with increasing concentration for each peak. The scale parameter ζ and the characteristic exponent γ indicate some relation to concentration, but their exact relations are not sought. In Figure 11, we show the change in the amplitude for four of the peaks with increase in the concentration of the inclusions. It is remarkable that around $q=0.6$ the behaviors of the spectral density functions change; the increase in the amplitude of the peaks with increasing concentration starts to decrease for increasing concentration as $q>0.6$. The ratio of the three selected peaks to peak 1 are shown in the inset of 11, and a similar activity is also observed in the ratio. Since the ratio between the amplitudes of the selected peaks does not indicate a simple linear relation to concentration q, the topological description of the system cannot be qualitatively investigated. However, as mentioned previously, the location and form of the spectral density functions greatly resemble each other, indicating that they are in fact related to the self-similar hierarchical nature of the composites expressed by LLL expression.

Table 1. Fit parameters of the Lévy distribution, Equation Equation 13, for the six peaks estimated by the numerical algorithm.

Peak	q	$\mathfrak{A}$	μ	γ	ζ	Peak	q	$\mathfrak{A}$	μ	γ	ζ
1	0.1	0.01	-2.47	1.58	6.18	2	0.1	0.06	-2.01	0.89	15.84
	0.2	0.02	-2.49	1.71	6.57		0.2	0.11	-2.02	1.10	12.74
	0.3	0.03	-2.50	1.66	7.29		0.3	0.18	-2.02	1.09	12.83
	0.4	0.04	-2.50	1.75	7.83		0.4	0.22	-2.02	1.20	11.81
	0.5	0.04	-2.50	1.94	7.53		0.5	0.26	-2.03	1.30	11.43
	0.6	0.04	-2.49	2.06	7.81		0.6	0.28	-2.03	1.33	11.46
	0.7	0.04	-2.49	1.98	7.86		0.7	0.28	-2.03	1.22	11.84
	0.8	0.04	-2.48	2.09	7.59		0.8	0.24	-2.03	1.26	11.63
	0.9	0.03	-2.48	1.79	7.20		0.9	0.16	-2.02	1.15	11.54
3	0.1	0.01	-1.62	2.97	8.23	4	0.1	0.23	-1.22	0.58	35.98
	0.2	0.02	-1.62	2.70	8.26		0.2	0.54	-1.24	0.52	49.50
	0.3	0.03	-1.61	2.21	9.34		0.3	0.82	-1.24	0.52	53.04
	0.4	0.03	-1.61	2.27	9.10		0.4	1.04	-1.24	0.56	49.16
	0.5	0.04	-1.61	2.51	8.84		0.5	1.20	-1.24	0.57	48.10
	0.6	0.04	-1.60	2.22	9.05		0.6	1.30	-1.24	0.60	47.57
	0.7	0.03	-1.61	2.58	8.91		0.7	1.26	-1.24	0.60	48.15
	0.8	0.02	-1.61	2.34	8.69		0.8	1.05	-1.24	0.62	45.98
	0.9	0.01	-1.62	2.43	8.32		0.9	0.58	-1.24	0.62	41.94
5	0.1	0.20	-0.62	1.85	10.89	6	0.1	0.22	-0.16	0.67	15.10
	0.2	0.36	-0.62	1.90	10.51		0.2	0.31	-0.14	0.77	10.98
	0.3	0.54	-0.62	1.70	11.15		0.3	0.38	-0.14	0.81	10.58
	0.4	0.62	-0.62	1.72	10.89		0.4	0.40	-0.14	0.80	9.84
	0.5	0.69	-0.61	1.58	10.90		0.5	0.37	-0.14	0.97	9.08
	0.6	0.65	-0.61	1.63	10.25		0.6	0.33	-0.14	1.02	8.87
	0.7	0.56	-0.61	1.66	9.57		0.7	0.28	-0.14	1.08	8.93
	0.8	0.39	-0.60	1.80	8.46		0.8	0.18	-0.14	1.33	8.22
	0.9	0.26	-0.60	1.50	9.33		0.9	0.10	-0.14	1.34	8.45

Table 1. Fit parameters of the Lévy distribution, Equation Equation 13, for the six peaks estimated by the numerical algorithm.

Peak	q	$\mathfrak{A}$	μ	γ	ζ	Peak	q	$\mathfrak{A}$	μ	γ	ζ
1	0.1	0.01	-2.47	1.58	6.18	2	0.1	0.06	-2.01	0.89	15.84
	0.2	0.02	-2.49	1.71	6.57		0.2	0.11	-2.02	1.10	12.74
	0.3	0.03	-2.50	1.66	7.29		0.3	0.18	-2.02	1.09	12.83
	0.4	0.04	-2.50	1.75	7.83		0.4	0.22	-2.02	1.20	11.81
	0.5	0.04	-2.50	1.94	7.53		0.5	0.26	-2.03	1.30	11.43
	0.6	0.04	-2.49	2.06	7.81		0.6	0.28	-2.03	1.33	11.46
	0.7	0.04	-2.49	1.98	7.86		0.7	0.28	-2.03	1.22	11.84
	0.8	0.04	-2.48	2.09	7.59		0.8	0.24	-2.03	1.26	11.63
	0.9	0.03	-2.48	1.79	7.20		0.9	0.16	-2.02	1.15	11.54
3	0.1	0.01	-1.62	2.97	8.23	4	0.1	0.23	-1.22	0.58	35.98
	0.2	0.02	-1.62	2.70	8.26		0.2	0.54	-1.24	0.52	49.50
	0.3	0.03	-1.61	2.21	9.34		0.3	0.82	-1.24	0.52	53.04
	0.4	0.03	-1.61	2.27	9.10		0.4	1.04	-1.24	0.56	49.16
	0.5	0.04	-1.61	2.51	8.84		0.5	1.20	-1.24	0.57	48.10
	0.6	0.04	-1.60	2.22	9.05		0.6	1.30	-1.24	0.60	47.57
	0.7	0.03	-1.61	2.58	8.91		0.7	1.26	-1.24	0.60	48.15
	0.8	0.02	-1.61	2.34	8.69		0.8	1.05	-1.24	0.62	45.98
	0.9	0.01	-1.62	2.43	8.32		0.9	0.58	-1.24	0.62	41.94
5	0.1	0.20	-0.62	1.85	10.89	6	0.1	0.22	-0.16	0.67	15.10
	0.2	0.36	-0.62	1.90	10.51		0.2	0.31	-0.14	0.77	10.98
	0.3	0.54	-0.62	1.70	11.15		0.3	0.38	-0.14	0.81	10.58
	0.4	0.62	-0.62	1.72	10.89		0.4	0.40	-0.14	0.80	9.84
	0.5	0.69	-0.61	1.58	10.90		0.5	0.37	-0.14	0.97	9.08
	0.6	0.65	-0.61	1.63	10.25		0.6	0.33	-0.14	1.02	8.87
	0.7	0.56	-0.61	1.66	9.57		0.7	0.28	-0.14	1.08	8.93
	0.8	0.39	-0.60	1.80	8.46		0.8	0.18	-0.14	1.33	8.22
	0.9	0.26	-0.60	1.50	9.33		0.9	0.10	-0.14	1.34	8.45

Figure 10. Spectral density functions of the Landau-Lifshitz/Looyenga equation for nine different concentrations, $q=\{0.1,\cdots ,0.9\}$. The data are shifted for better comparison. Six peaks are resolved between $x={10}^{-3}$ and $x=1$. Interestingly, the positions of the most probable spectral parameters of the six are not altered with increasing concentration. Logarithmic scale is base 10.

Figure 10. Spectral density functions of the Landau-Lifshitz/Looyenga equation for nine different concentrations, $q=\{0.1,\cdots ,0.9\}$. The data are shifted for better comparison. Six peaks are resolved between $x={10}^{-3}$ and $x=1$. Interestingly, the positions of the most probable spectral parameters of the six are not altered with increasing concentration. Logarithmic scale is base 10.

Figure 11. Spectral density function values at constant x; $log\left(x\right)=-2.5$ (★), $log\left(x\right)=-2$ (•), $log\left(x\right)=-1.275$ (▾) and $log\left(x\right)=-0.675$ (▪). The inset illustrates the ratios for $log\left\{\mathsf{g}\right[log\left(x\right)=-2.5\left]\right\}/log\left\{\mathsf{g}\right[log\left(x\right)=-2\left]\right\}$ (□), $log\left\{\mathsf{g}\right[log\left(x\right)=-2.5\left]\right\}/log\left\{\mathsf{g}\right[log\left(x\right)=-1.275\left]\right\}$ (▿) and $log\left\{\mathsf{g}\right[log\left(x\right)=-2.5\left]\right\}/log\left\{\mathsf{g}\right[log\left(x\right)=-0.675\left]\right\}$ (∘). Logarithmic scale is base 10.

Figure 11. Spectral density function values at constant x; $log\left(x\right)=-2.5$ (★), $log\left(x\right)=-2$ (•), $log\left(x\right)=-1.275$ (▾) and $log\left(x\right)=-0.675$ (▪). The inset illustrates the ratios for $log\left\{\mathsf{g}\right[log\left(x\right)=-2.5\left]\right\}/log\left\{\mathsf{g}\right[log\left(x\right)=-2\left]\right\}$ (□), $log\left\{\mathsf{g}\right[log\left(x\right)=-2.5\left]\right\}/log\left\{\mathsf{g}\right[log\left(x\right)=-1.275\left]\right\}$ (▿) and $log\left\{\mathsf{g}\right[log\left(x\right)=-2.5\left]\right\}/log\left\{\mathsf{g}\right[log\left(x\right)=-0.675\left]\right\}$ (∘). Logarithmic scale is base 10.

Figure 12. Peak positions, $\mu =log\left(\overline{x}\right)$. The solid line is a linear fit with $\mu =0.470P_j-2.99$, where j is the peak number in Table 1.

Figure 12. Peak positions, $\mu =log\left(\overline{x}\right)$. The solid line is a linear fit with $\mu =0.470P_j-2.99$, where j is the peak number in Table 1.

Finally, the statistical analysis of the concentration calculated from the integral of the spectral density function (integral expression on the right-hand side of Equation 2) and the value of the percolation strength ${\xi }_s$ at each Monte Carlo step are summarized in Figure 13 and Figure 14. The number distribution of the integral of the spectral density function calculated at each Monte Carlo cycle is not centered at the actual concentration q taken but deviated. The deviation indicates that there is a percolation path network structure, cf. Figure 13. As expected from the definition of the spectral density function, Equation (20), addition of the estimated concentration $\overline{q}$ and percolation strength $\overline{{\xi }_s}$ yield numerical values very close to the actual concentration q. The value $\overline{q}$ is calculated from the distribution. It is striking that the applied numerical method is capable of estimating the concentration of the filler material when it is not known in advance.

Figure 13. The number density distribution of estimated $\int \mathsf{g}\left(x\right)\mathrm{d}x\equiv \overline{q}$ for two concentration q levels, $q=\{0.3,0.4\}$. The dashed vertical lines are the expectation values of $\overline{q}$, which are $0.255\pm 0.065$ and $0.310\pm 0.096$ for concentration levels $0.3$ and $0.4$, respectively. Logarithmic scale is base 10.

Figure 13. The number density distribution of estimated $\int \mathsf{g}\left(x\right)\mathrm{d}x\equiv \overline{q}$ for two concentration q levels, $q=\{0.3,0.4\}$. The dashed vertical lines are the expectation values of $\overline{q}$, which are $0.255\pm 0.065$ and $0.310\pm 0.096$ for concentration levels $0.3$ and $0.4$, respectively. Logarithmic scale is base 10.

Figure 14. The expected values of constants, $\overline{q}$ and $\overline{{\xi }_s}$, in the spectral density approach. The integral of $\mathsf{g}\left(x\right)$ yields the concentration of the inclusions $\overline{q}$, denoted by open squares (□). The expected value of percolation strength is $\overline{\xi }$ denoted by open circles (∘). Both values are also presented with error bars. Addition of the two expected values $\overline{{\xi }_s}+\overline{q}$, denoted by stars (★), leads to the actual concentration q as in Equation (20), which is denoted by the solid line (——–). In the inset, the number distributions of ${\xi }_s$, percolation strength, estimated in each Monte Carlo cycle for $q=0.3$ and $q=0.4$ are illustrated with the most expected $\overline{{\xi }_s}$ being $0.041\pm 0.006$ and $0.083\pm 0.015$, respectively, for the two concentrations.

Figure 14. The expected values of constants, $\overline{q}$ and $\overline{{\xi }_s}$, in the spectral density approach. The integral of $\mathsf{g}\left(x\right)$ yields the concentration of the inclusions $\overline{q}$, denoted by open squares (□). The expected value of percolation strength is $\overline{\xi }$ denoted by open circles (∘). Both values are also presented with error bars. Addition of the two expected values $\overline{{\xi }_s}+\overline{q}$, denoted by stars (★), leads to the actual concentration q as in Equation (20), which is denoted by the solid line (——–). In the inset, the number distributions of ${\xi }_s$, percolation strength, estimated in each Monte Carlo cycle for $q=0.3$ and $q=0.4$ are illustrated with the most expected $\overline{{\xi }_s}$ being $0.041\pm 0.006$ and $0.083\pm 0.015$, respectively, for the two concentrations.

6. Application of the Havriliak-Negami Expression

Previously, it has been stated that the scaled permittivity ξ in Equation (3) can be expressed as in the conventional form as in the case of dielectric relaxation [21], Equation (5). To illustrate and verify this statement, we apply a complex nonlinear least-squares curve fit algorithm to the scaled permittivity of the LLL expression, denoted as ${\xi }^{\mathrm{LLL}}$ below. The error in the curve fitting procedure is used to quantify the fitness of the model function in Equation (5). The error is calculated as the sum of the relative error at each point as follows,

$$\begin{array}{c}\hfill \mathfrak{E}=\sum {\left[\frac{\Re \left({\xi }^{\mathrm{LLL}}\right)-\Re \left({\xi }^{\mathrm{HN}}\right)}{\Re \left({\xi }^{\mathrm{LLL}}\right)}\right]}^2+{\left[\frac{\Im \left({\xi }^{\mathrm{LLL}}\right)-\Im \left({\xi }^{\mathrm{HN}}\right)}{\Im \left({\xi }^{\mathrm{LLL}}\right)}\right]}^2\end{array}$$

(14)

Here ${\xi }^{\mathrm{HN}}$ is the model expression of Equation (5). The fit results are listed in Table 2, where the model values for the spectral parameter x, concentration q and percolation strength ${\xi }_s$ are presented with over-lines as in the previous section.

Table 2. Fit parameters of the Havriliak-Negami [90] expression to scaled permittivity, Equation (5). The error $\mathfrak{E}$ is calculated with Equation (14).

q	α	β	$\overline{x}$	$\overline{q}$	$\overline{{\xi }_s}$	$\mathfrak{E}$
0.1	1.034	0.438	0.536	0.095	0.000	1.19$\times {10}^{-3}$
0.2	0.822	0.555	0.514	0.200	0.009	7.28$\times {10}^{-6}$
0.3	0.716	0.635	0.464	0.297	0.032	2.45$\times {10}^{-4}$
0.4	0.850	0.454	0.522	0.345	0.066	4.27$\times {10}^{-6}$
0.5	0.800	0.470	0.494	0.391	0.128	2.06$\times {10}^{-5}$
0.6	0.828	0.422	0.493	0.397	0.216	4.08$\times {10}^{-6}$
0.7	0.822	0.409	0.475	0.371	0.342	5.20$\times {10}^{-6}$
0.8	0.815	0.399	0.456	0.300	0.510	6.41$\times {10}^{-6}$
0.9	0.810	0.390	0.438	0.179	0.727	7.60$\times {10}^{-6}$

Table 2. Fit parameters of the Havriliak-Negami [90] expression to scaled permittivity, Equation (5). The error $\mathfrak{E}$ is calculated with Equation (14).

q	α	β	$\overline{x}$	$\overline{q}$	$\overline{{\xi }_s}$	$\mathfrak{E}$
0.1	1.034	0.438	0.536	0.095	0.000	1.19$\times {10}^{-3}$
0.2	0.822	0.555	0.514	0.200	0.009	7.28$\times {10}^{-6}$
0.3	0.716	0.635	0.464	0.297	0.032	2.45$\times {10}^{-4}$
0.4	0.850	0.454	0.522	0.345	0.066	4.27$\times {10}^{-6}$
0.5	0.800	0.470	0.494	0.391	0.128	2.06$\times {10}^{-5}$
0.6	0.828	0.422	0.493	0.397	0.216	4.08$\times {10}^{-6}$
0.7	0.822	0.409	0.475	0.371	0.342	5.20$\times {10}^{-6}$
0.8	0.815	0.399	0.456	0.300	0.510	6.41$\times {10}^{-6}$
0.9	0.810	0.390	0.438	0.179	0.727	7.60$\times {10}^{-6}$

The fit results are shown in Figure 15 for the $q=0.3$ case in a form similar to Figure 8. There are actually no particular differences between the two methods, except that the numerical techniques based on the Monte Carlo algorithm are capable of resolving individual peaks, as shown in the comparison graph in Figure 16. It is not clear, for example, from the Havriliak-Negami approach that the system indicates a self-similar fractal-like structure. One should note that the Havriliak-Negami distribution, cf. § 7. [73], presented below, was spread evenly over spectral parameter values larger than one, $x>1$, which is not possible in the spectral density representation. However, when analyzing the data it is very convenient and trivial to implement available curve fitting programs. Similar to Figure 14, the estimated concentration $\overline{q}$ and percolation strength $\overline{{\xi }_s}$ from the parametric analysis satisfy the condition in Equation 20 as shown in Figure 17.

7. Conclusions

In this review, we first presented the significance of the spectral density representation and then derived an expression for dielectric mixtures that resembles the distribution of relaxation times representation in dielectric relaxation phenomenon. In the derivation we used the spectral density representation of Ghosh and Fuchs [33]. It is shown that existing knowledge on dielectric relaxation theory can be applied to dielectric properties of composites. To confirm the hypothesis, both a similar method for estimating the distribution of relaxation times and an extensively used empirical formula to express dielectric relaxation are employed to estimate the spectral density functions of composites simulated with the LLL expression. The numerical method, based on the Monte Carlo technique, estimated a couple of peaks that did not change their location in the spectra with increased concentration of inclusions. This static behavior of the spectra indicates that there exists a hierarchical structural order in the composite. As a result, we infer that the LLL expression is proper for systems with a self-similar fractal nature, such as composites with colloid aggregates and porous materials. We have explicitly shown why the LLL expression could be applied to describe the dielectric properties of powdered and porous systems. Last but not least, the findings are significant to confirm the structure of composite systems, whose dielectric permittivities are described in the literature with the LLL expression.

Figure 15. (a) The real and (b) the imaginary parts of the dielectric permittivity calculated with Equation (12) for $q=0.3$ and the modeled response obtained with the application of the Havriliak-Negami expression to the scaled permittivity. The permittivities of the matrix and the inclusion phases are the same as in Figure 8. (c) Argand diagram of scaled permittivity without the percolation strength contribution. The inset in (a) shows the dielectric susceptibility after the substruction of ohmic conductivity and the permittivity at high frequencies. The lines with points are the simulated data, and the symbols (∘) are the data estimates of the Havriliak-Negami expression, Equation 5. There is very good agreement between the simulated and the analyzed data sets. Logarithmic scale is base 10.

Figure 15. (a) The real and (b) the imaginary parts of the dielectric permittivity calculated with Equation (12) for $q=0.3$ and the modeled response obtained with the application of the Havriliak-Negami expression to the scaled permittivity. The permittivities of the matrix and the inclusion phases are the same as in Figure 8. (c) Argand diagram of scaled permittivity without the percolation strength contribution. The inset in (a) shows the dielectric susceptibility after the substruction of ohmic conductivity and the permittivity at high frequencies. The lines with points are the simulated data, and the symbols (∘) are the data estimates of the Havriliak-Negami expression, Equation 5. There is very good agreement between the simulated and the analyzed data sets. Logarithmic scale is base 10.

Figure 16. Comparison of the distributions $sfg$ obtained by the Havriliak-Negami approach (– – –) and the novel numerical approach (——–). Observe that the Havriliak-Negami distribution also considers spectral parameter x values larger than 1, $log\left(x\right)>0$. Logarithmic scale is base 10.

Figure 16. Comparison of the distributions $sfg$ obtained by the Havriliak-Negami approach (– – –) and the novel numerical approach (——–). Observe that the Havriliak-Negami distribution also considers spectral parameter x values larger than 1, $log\left(x\right)>0$. Logarithmic scale is base 10.

Figure 17. The most probable values of constants $\overline{q}$ and $\overline{{\xi }_s}$ from application of the Havriliak-Negami expression. The integral of $\mathsf{g}\left(x\right)$ yields the concentration of the inclusions $\overline{q}$, denoted by open squares $(\square )$. The percolation strength is $\overline{\xi }$, denoted by open circles (∘). The sum of the two constants, $\overline{{\xi }_s}+\overline{q}$, denoted by stars (★), leads to the actual concentration, which is denoted with the solid line (——–).

Figure 17. The most probable values of constants $\overline{q}$ and $\overline{{\xi }_s}$ from application of the Havriliak-Negami expression. The integral of $\mathsf{g}\left(x\right)$ yields the concentration of the inclusions $\overline{q}$, denoted by open squares $(\square )$. The percolation strength is $\overline{\xi }$, denoted by open circles (∘). The sum of the two constants, $\overline{{\xi }_s}+\overline{q}$, denoted by stars (★), leads to the actual concentration, which is denoted with the solid line (——–).