Rigorously solvable model for the electrical conductivity of dispersions of hard-core–penetrable-shell particles and its applications

# Rigorously solvable model for the electrical conductivity of dispersions of hard-core–penetrable-shell particles and its applications

M. Ya. Sushko    A. K. Semenov Department of Theoretical Physics and Astronomy, Mechnikov National University, 2 Dvoryanska St., Odesa 65026, Ukraine
July 20, 2019
###### Abstract

We generalize the compact group approach to conducting systems to give a self-consistent analytical solution to the problem of the effective quasistatic electrical conductivity of macroscopically homogeneous and isotropic dispersions of hard-core–penetrable-shell particles. The shells are in general inhomogeneous and characterized by a radially-symmetrical, piecewise-continuous conductivity profile. The local value of the conductivity is determined by the shortest distance from the point of interest to the nearest particle. The effective conductivity is expressed in terms of the constituents’ conductivities and volume concentrations; the latter account for the statistical microstructure of the system. The theory effectively incorporates many-particle effects and is expected to be rigorous in the static limit. Using the well-tested statistical physics results for the shell volume concentration, this conclusion is backed up by mapping the theory on available 3D random resistor network simulations for hard spheres coated with fully penetrable concentric shells. Finally, the theory is shown to fit experimental data for real composite solid electrolytes. The fitting results indicate that the effect of enhanced electrical conduction is generally contributed to by several mechanisms. These are effectively taken into account through the shell conductivity profile.

###### pacs:
42.25.Dd , 77.22.Ch, 77.84.Lf, 82.70.-y

## I Introduction

The objectives of this paper are threefold: (1) to develop a homogenization theory for the effective quasistatic electrical conductivity of macroscopically homogeneous and isotropic particulate substances and dispersions of particles with the core-shell morphology; (2) to test the theory by comparing its predictions with available results of random resistor network (RRN) simulations; and (3) to exemplify the applicability of the theory to real systems by processing experimental data for composite solid electrolytes (CSEs) prepared by dispersing fine insulating particles into matrix ionic conductors.

The indicated class of composites attracts a special attention due to nontrivial behavior of their . Through the addition of filler particles (for instance, alumina particles, with electrical conductivity ) to matrix ionic conductors (such as polycrystalline metal halides, whose typical electrical conductivities ), of the resulting CSEs can be increased dramatically, by one to three orders of magnitude as compared to . This effect is called enhanced ionic conduction. Since its discovery by Liang Liang1973 () in polycrystalline lithium iodide containing alumina particles, it has been observed in dozens of CSEs (for a detailed bibliography, see reviews Wagner1980 (); Takahashi1989 (); Dudney1989 (); Maier1995 (); Agrawal1999 (); Uvarov2011 (); Kharton2011 (); Gao2016 ()) and composite polymer-based electrolytes (Knauth2008 (); Sequeira2010 ()) to make those promising materials for electrolytic applications.

Experiment also reveals that the maximum conduction enhancement usually occurs as the filler volume concentration reaches values in between 0.1 and 0.4. It is followed by a decrease in as is further increased. Such a nonmonotonic dependence of upon is a challenging problem for homogenization theory, since the existing approaches to two-phase systems, such as the classical Maxwell-Garnet Maxwell1873 (); Maxwell1904 () and Bruggeman Bruggeman1935 (); Landauer1952 () mixing rules, their numerous modifications (see Bohren1983 (); Bergman1992 (); Sihvola1999 (); Tsang2001 (); Torquato2013 (); Milton2004 ()), cluster expansions Torquato1984 (); Torquato1985b () for dispersions of spheres with arbitrary degree of impenetrability, their extensions Torquato1985 () with the Padé approximant technique, and systematic simulations MyroshnychenkoPRE2005 (); MyroshnychenkoJAP2005 (); Myroshnychenko2008 (); Myroshnychenko2010 () of random 2D systems of hard-core–penetrable-shell discs by combining Monte Carlo algorithms and finite element calculations do not exhibit it. The reason is that two-phase models oversimplify the actual microstructure of CSEs and disregard the processes involved.

A typical way out is to model a CSE as a three-phase system and determine by solving a pertinent homogenization problem. The solution is expressed in terms of the geometric and electric parameters of the phases. These parameters are estimated so as to incorporate the relevant physical effects and account for the observed behavior of . Several classes of such models have been proposed.

(i) Cubic lattices of cubic insulators surrounded by highly conductive layers and embedded in a conductive material Wang1979 (); Jiang1995a (); Jiang1995b (). The arrangement of particles on a simple cubic lattice makes it possible to represent the system with a resistor network and to calculate for the entire range of . The results suggest that the conductivity in the layer outside each particle may have a maximum at certain distance away from the surface.

(ii) Three-component resistor models with a matrix represented by normally conducting bonds, the inert randomly distributed (quadratic or cubic) particles by insulating bonds, and the interface region by highly conducting bonds Bunde1985 (); Roman1986 (); Blender1987 (). The models are solved by Monte Carlo simulations or a position-space reorganization technique and exhibit two threshold concentrations of the insulating material. One corresponds to the onset of interface percolation and the second one to a conductor-insulator transition. The effective medium and continuum percolation approaches to these models are discussed in Rojo1988 () and Roman1987 (); Roman1990 (), respectively.

(iii) Random three-phase dispersions of spherical particles comprising hard cores coated with concentric shells, either hard or penetrable (see Fig. 1), of potentially higher conductivity. Such core-shell models better suit the physical conditions in CSEs, but are harder to analyze. Analytical studies, such as Stoneham1979 (); Brailsford1986 (); Nan1991L (); Nan1991 (); Nan1993 (); Wiec1994 (), of them are usually limited to the case of hard shells, involve a sequence of one-particle approximations, and repeatedly use the Maxwell-Garnet Maxwell1873 (); Maxwell1904 () or/and Bruggeman Bruggeman1935 (); Landauer1952 () mixing rules. The case of penetrable shells has been attacked through RRN simulations, such as Siekierski2005 (); Siekierski2006 (); Siekierski2007 () for mono- and Kalnaus2011 () for polysized particles. The essential details of these core-shell model results are scrutinized in Sec. VIII.

In what follows, we derive a self-consistent analytic many-particle solution for of macroscopically homogeneous and isotropic 3D model dispersions of hard-core–penetrable-shell spheres, the shells being, in the general case, electrically inhomogeneous and characterized by radially-symmetrical, piecewise-continuous conductivity profiles (see Fig. 1 for the details of the model). The desired is a functional of the constituents’ conductivities and volume concentrations that satisfies a certain integral relation, rigorous in the static limit. The volume concentrations account for the statistical microstructure of the system.

The derivation is carried out using the compact groups approach (CGA) Sushko2007 (); Sushko2009CompGroups (); Sushko2009AnisPart (); Sushko2017 (). It was originally designed to efficiently take into account many-particle polarization and correlation effects, without an in-depth modeling of those, in concentrated dielectric dispersions. In this paper, elaborating its statistical-averaging version, we (1) generalize the CGA to conducting systems whose constituents have complex permittivities with first-order poles at frequency ; (2) scrutinize, for such systems, the passage to the (quasi)static limit in all terms of the iterative series for the averaged electric field and current; (3) bring new arguments, not restricted to dielectric systems, for the internally-consistent closure of the homogenization procedure and determination of the complex permittivity of the auxiliary host matrix; and (4) propose a technique for dealing with inhomogeneous overlapping regions.

Using the well-tested statistical physics results Rikvold85 (); Lee1988 (); Rottereau2003 () for the shell volume concentration, we then validate the solution obtained by mapping it onto the entire set of available 3D random resistor network simulation data Siekierski2005 (); Siekierski2006 (); Siekierski2007 () for dispersions of hard spheres coated with fully penetrable (electrically uniform or inhomogeneous) concentric shells. The solution is capable of recovering all of these data in the entire ranges of simulated. To our best knowledge, no such an analytic solution has been offered so far even for the simplest case of uniform shells with their conductivity being equal to that of the cores.

Finally, we apply the model to real CSEs. The results of processing experimental data Liang1973 () clearly indicate that the concept of inhomogeneous penetrable shells provides an efficient way for describing the net effect on by different mechanisms. The latter may contribute most significantly in different ranges of . If so, they are accounted for by different parts in the model shell conductivity profile. This fact opens new opportunities for scrutinizing the physics of processes in real composites.

The paper is arranged as follows. Some basic equations and definitions of macroscopic electrodynamics for media with complex permittivities of the constituents are recalled in Sec. II. With those in mind, the CGA is generalized in Sec. III to the problem of the effective quasistatic complex permittivity of macroscopically homogeneous and isotropic dispersions. The governing equation for is expressed in terms of the statistical moments for the local deviations of the permittivity distribution in the dispersion from the complex permittivity of the host in the auxiliary system. By requiring that the CGA and boundary conditions Sillars1937 () for complex electric fields be compatible, is determined in Sec. IV. The calculations of for dispersions of isotropic core-shell particles with electrically homogeneous and inhomogeneous shells are performed in Secs. V and VI, respectively. The resultant equations for are presented in Sec. VII. Their validity is shown in Sec. VIII by mapping their solutions onto extensive RRN simulation data Siekierski2005 (); Siekierski2007 (); Siekierski2006 (). The applicability of the theory to real CSEs Liang1973 () is discussed in Sec. IX. The main results of the paper are summarized in Sec. X.

## Ii Basic equations and definitions

Maxwell’s macroscopic equations for a nonmagnetic heterogeneous medium have the form

 divD=4πρ,divH=0, (1)
 curlE=−1c∂H∂t,curlH=4πcj+1c∂D∂t, (2)

where , , , , , and are, respectively, the electric field, electric displacement, magnetic field, free charge density, free current density, and speed of light in vacuum. The densities and are related by the continuity equation

 ∂ρ∂t+divj=0. (3)

Assuming the linear constitutive equations

 D=εE,j=σE, (4)

where and are the local (real) dielectric constant and electrical conductivity in the medium, let us write out the equations for the spatial field distributions caused by time-harmonic (, being the imaginary unit) probing radiation whose working frequencies are sufficiently small to neglect any dielectric relaxation phenomena.

Differentiating both sides of the first Eq. (1) with respect to time and using Eqs. (3) and (4), we immediately obtain

 divJ=0, (5)

where

 J=−iω4π^εE (6)

and

 ^ε=ε+i4πσω (7)

can be interpreted as the complex current density and the complex permittivity of the medium, respectively. In the quasistatic limit (), Eq. (6) reduces to Ohm’s law, given by the second Eq. (4), and Eq. (5) to the second-order partial differential equation

 ∇⋅(σ∇Φ)=0, (8)

which is typically used to find the electric potential distribution , , in homogenization problems.

In terms of , the second Eq. (2) reads

 curlH=−iωc^εE.

Differentiating this with respect to time and using the first Eq. (2), we find the equation for the electric field distribution:

where is the magnitude of the wave vector of the incident field in vacuum.

## Iii Compact group approach to homogenization of conducting systems

The main points of this approach in application to macroscopically homogeneous and isotropic nonconducting systems are discussed in detail in Sushko2007 (); Sushko2009CompGroups (); Sushko2009AnisPart (); Sushko2017 (). Here, closely following the summary in Sushko2017 (), we outline a generalization of the approach to macroscopically homogeneous and isotropic dispersions comprising conducting constituents, such as conducting dielectrics or imperfectly insulating materials. We assume that for a given , the complex permittivities of all the constituents have structure (7), where and are in general piecewise-continuous and bounded real functions of spacial coordinates; and that its can be calculated based upon the following suggestions:

(1) is equivalent, in its response to a long-wavelength probing field (), to an auxiliary system prepared by embedding the constituents (particles and matrix) of into a uniform host (perhaps, imagined) with some permittivity .

(2) can be viewed as a set of compact groups of both particles and regions occupied by the real matrix. The compact groups are defined as macroscopic regions whose typical sizes are much smaller than the wavelength of probing radiation in , but which yet include sufficiently large numbers of particles to remain macroscopic and retain the properties of the entire .

(3) The complex permittivity distribution in is

 ^ε(r)=^εf+δ^ε(r) (10)

where is the contribution from a compact group located at point . The explicit form of is modeled according to the geometrical and electrical parameters of ’s constituents.

(4) can be found as the proportionality coefficient in the relation

 ⟨J(r)⟩=−iω4π⟨^ε(r)E(r)⟩=−iω4π^εeff⟨E(r)⟩, (11)

where and are the local values of the complex current and electric field, respectively, and the angle brackets stand for the ensemble averaging. In the quasistatic limit, provided and is real-valued, Eq. (11) reduces to the common definition Bergman1992 (); Sihvola1999 (); Torquato2013 () of the effective conductivity:

 ⟨j(r)⟩=⟨σ(r)E(r)⟩=σeff⟨E(r)⟩. (12)

(5) The electric field distribution in obeys the equation

which directly follows from Eqs. (9) and (10). The equivalent integral equation is

 E(r)=E0(r)−∫Vdr′T(|r−r′|)k20δ^ε(r′)E(r′), (14)

where , , and (with ) are, respectively, the incident wave field, its amplitude, and its wave vector in , and is the Green’s tensor of Eq. (13).

(6) The formal solution for and those for , , and are representable in the form of infinite iterative series. For systems whose constituents have the permittivities of form (7), the functions in the integrands and also the function remain bounded even at , where as well. Mathematically, the situation is identical to that for nonconducting systems and can be treated analogously. Namely, in the iterative series for and , each under the integral sign is replaced by its decomposition , derived for a spherical exclusion volume of radius , into a Dirac delta function singular part and a principal value part Ryzhov1965 (); Weiglhofer1989 ():

 ˜Tαβ(r)=13k2δαβδ(r)eikr+14πk2(1r3−ikr2) ×(δαβ−3eαeβ)eikr−14πr(δαβ−eαeβ)eikr, (15)

where is the Dirac delta function, is the Kronecker delta, and is the -component of the unit vector . The contributions to made by the subseries containing, in their integrands, the principal value parts are estimated to be of the order , at most Sushko2007 (), as compared to those made by the subseries with only the Dirac delta function parts. For a finite typical linear size of the system, the former can be decreased below any preset value by taking a sufficiently small . So, passing to the limit and formally replacing each in the integrals for and by its Delta function part, we obtain

 limω→0⟨E(r)⟩=limω→0[1+⟨^Q(r)⟩]E0, (16)
 limω→0⟨J(r)⟩=−ilimω→0ω^εf4π[1−2⟨^Q(r)⟩]E0, (17)

where

 ^Q(r)≡∞∑s=1(−13^εf)s(δ^ε(r))s. (18)

(7) The above results can be obtained without resort to iterative series. Indeed, it follows from Eq. (III) that

 limω→0k20^εf˜Tαβ(r)=τ(1)αβ+τ(2)αβ, (19)

where

 τ(1)αβ=13δαβδ(r),τ(2)αβ=δαβ−3eαeβ4πr3. (20)

Substituting Eqs. (19) and (20) into Eq. (14), making simple algebraic manipulations and statistical averaging, and implying that , we obtain

 (21)
 ⟨J(r)⟩=−iω4π^εf[1+2⟨δ^ε(r)3^εf+δ^ε(r)⟩]E0+i34π∫Vdr′τ(2)(|r−r′|)⟨ω^ε(r)δ^ε(r′)3^εf+δ^ε(r)E(r′)⟩. (22)

For macroscopically isotropic and homogeneous systems, two-point statistical averages depend only on . Due to this symmetry and because of a special form of the angular dependence of , the integrals in Eqs. (21) and (22) vanish. Finally, viewing the expressions in the angle brackets as the sums of infinite geometric series, we arrive at Eqs. (16)–(18).

This consideration is very similar to that used in the strong-property-fluctuation theory (SPFT) Tsang2001 (); Ryzhov1965 (); Ryzhov1970 (); Tamoikin1971 (); Tsang1981 (); Zhuck1994 (); Michel1995 (); Mackay2000 (); Mackay2001 (). However, our theory gives another interpretation to , appeals to the macroscopic symmetry of the entire system instead of the symmetry of correlation functions, and postulates no condition on the stochastic field in order to improve the convergence of the iteration procedure and decide on the value of . Once the latter is determined, the analysis of reduces to modeling , calculating its moments , and finding their sum in Eqs.  (16) and (17).

## Iv Determination of ^εf

If the permittivities of the constituents, that of , and have at the structure (7), then, at least in this limit, it is the Bruggeman-type of homogenization that is compatible with the formalism of the CGA and definition (11). To prove this statement, we first remind that is the amplitude of the probing electric field in the uniform fictitious matrix of permittivity , and is the effective electric field in the homogenized dispersion of permittivity , caused by the same probing field. Next, we recall the boundary condition Sillars1937 ()

 ^ε1E1n=^ε2E2n

for the normal components of complex electric fields at the surface between two conducting dielectrics (or imperfectly insulating materials) with permittivities of type (7). For the surface between the fictitious matrix and the homogenized dispersion, it gives

 (23)

This relation, Eq. (11), and Eqs. (16)–(18) yield, at , the system of equations

 ^εf=^εeff(1+⟨^Q⟩),^εf(1−2⟨^Q⟩)=^εeff(1+⟨^Q⟩). (24)

Since , it follows immediately that

 ^εf=^εeff (25)

and, for this ,

 ⟨^Q(r)⟩=0. (26)

The latter is the desired governing equation for . It is valid in the limit .

## V Statistical moments ⟨(δ^ε(r))s⟩

We consider a dispersion of spherically symmetrical core-shell particles embedded into a uniform matrix. Suppose that the local permittivity value at a point within the dispersion is determined by the distance from to the center of the nearest ball as

 ^ε(r)=⎧⎨⎩^ε1iflR2. (27)

Here is the radius of the core of complex permittivity , is the outer radius of the shell of complex permittivity , and is the complex permittivity of the matrix. Within the CGA Sushko2007 (); Sushko2009CompGroups (); Sushko2009AnisPart (); Sushko2017 (), such a system can be modelled as follows.

Let be the Heaviside step function and () be the characteristic functions of balls centered at point and having radii . Suggesting that and allowing the balls to overlap, consider the complex permittivity distribution of form (10) with

 δ^ε(r) = Π1(r)Δ^ε1+[Π2(r)−Π1(r)]Δ^ε2 (28) +[Π0(r)−Π2(r)]Δ^ε0,

where , and each

 Πi(r)=1−N∏a=1(1−χ(i)a(r)) (29)

is the characteristic (indicator) function of the collection of balls of radius . In the limit , and Eq. (28) leads to the model permittivity distribution (27) for a dispersion of penetrable core-shell particles embedded into a uniform matrix of permittivity . The characteristic function of the entire region occupied by the substance of permittivity in this dispersion is given by the coefficient function in front of the corresponding . It is readily verified that the characteristic functions of regions with different permittivities are mutually orthogonal.

Further, we limit ourselves to the case of particles with hard cores. Then where is the Kronecker delta, and reduces to

 Π1(r)=N∑a=1χ(1)a(r), (30)

with the additional restriction on the locations of any two balls.

For a macroscopically homogeneous and isotropic system

where is the volume concentration of the hard cores. In view of this fact and the mutual orthogonality of the characteristic functions of regions with different permittivities, the moments of the function (28) can be represented in the limit as

 ⟨[δ^ε(r)]s⟩ = c(Δ^ε1)s+(ϕ(c,δ)−c)(Δ^ε2)s (31) +(1−ϕ(c,δ))(Δ^ε0)s,

where

 ϕ(c,δ) ≡ ⟨Π2(r)⟩=⟨[1−N∏a=1(1−χ(2)a(r))]⟩ = ⟨∑1≤a≤Nχ(2)a(r)−∑1≤a

is the effective volume concentration of hard-core–penetrable-shell particles Torquato2013 (). Besides , it depends on the relative thickness of the shell ; in particular, . The averaged values of the sums in Eq. (LABEL:eq:effectiveconcentration) are calculated using the partial distribution functions for the system under consideration.

For hard-core–hard-shell particles Eq. (LABEL:eq:effectiveconcentration) gives

 ϕ(c,δ)=c(1+δ)3. (33)

## Vi The Case of Inhomogeneous Isotropic Shells

To extend the results of Sec. V to dispersions of spherically symmetrical particles with hard cores and adjacent inhomogeneous penetrable shells, we begin with the situation where each shell consists of concentric spherical layers with outer radii (grouped in the order of increasing magnitude) and constant dielectric permittivities , . Next, we suggest that the local permittivity distribution within the dispersion is given by this law, generalizing Eq. (27):

 ^ε(r)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩^ε1iflR2,M. (34)

Let be the characteristic functions of balls centered at point , having radii , and allowed to overlap. Then the characteristic functions of the collections of balls with radii are

 Π2,j(r)=1−N∏a=1(1−χ(2,j)a). (35)

Repeating almost literally the reasoning in Sec. V, we can represent the distribution (34) in form (10) with

 δ^ε(r) = [1−Π2,M(r)]Δ^ε0+Π1(r)Δ^ε1 (36) + [Π2,1(r)−Π1(r)]Δ^ε2,1 + M∑j=2[Π2,j(r)−Π2,j−1(r)]Δ^ε2,j,

where . Correspondingly,

 ⟨[δ^ε(r)]s⟩ = [1−ϕ(c,δM)](Δ^ε0)s+c(Δ^ε1)s (37) + M∑j=1[ϕ(c,δj)−ϕ(c,δj−1)](Δ^ε2)s,

where , is given by Eq. (LABEL:eq:effectiveconcentration) at , and we denoted . Finally, passing to the limits , () and assuming to be differentiable with respect to , for a dispersion of particles with a piecewise-continuous complex permittivity profile of the shells we obtain

 ⟨[δ^ε(r)]s⟩ = [1−ϕ(c,δM)](Δ^ε0)s+c(Δ^ε1)s (38) + δM∫0∂ϕ(c,u)∂u[Δ^ε2(u)]sdu,

where is the deviation as a function of and corresponds to the outermost edge of the shell.

For a uniform shell (), Eq. (38) immediately reduces to Eq. (31) with .

## Vii Equations for Effective Conductivity

In the case of uniform shells, where the moments are given by Eq. (31), the sums involved in Eq. (26) take the form . For , they reduce to infinite geometric series, so

 ∞∑s=1(−Δ^εi3^εeff)s=−^εi−^εeff2^εeff+^εi. (39)

For , the left-hand side in Eq. (39) can be treated, as was shown in Sec. III, as an asymptotic series of the right-hand side, so that the restriction can be omitted. The resulting equation for is

 [1−ϕ(c,δ)]^ε0−^εeff2^εeff+^ε0+c^ε1−^εeff2^εeff+^ε1 +[ϕ(c,δ)−c]^ε2−^εeff2^εeff+^ε2=0. (40)

To extract the equation for the quasistatic , we pass in Eq. (40) to the limit and assume that

 2σeff+σi≫ω4π(2εeff+εi),|σi−σeff|≫ω4π|εi−εeff|. (41)

Then, retaining the first term in the formal perturbation series in for the left-hand side of Eq. (40), we obtain

 [1−ϕ(c,δ)]σ0−σeff2σeff+σ0+cσ1−σeff2σeff+σ1 +[ϕ(c,δ)−c]σ2−σeff2σeff+σ2=0. (42)

In view of Eq. (41), the sufficient condition for the validity of Eq. (42) can be represented as

 |σeff−σi|≫ω4π(2εeff+εi),i=1,2,3. (43)

The generalizations of Eqs. (40) and (42) to dispersions of particles with inhomogeneous isotropic shells are evident [see Eq. (38)]:

 [1−ϕ(c,δM)]^ε0−^εeff2^εeff+^ε0+c^ε1−^εeff2^εeff+^ε1 +δM∫0∂ϕ(c,u)∂u^ε2(u)−^εeff2^εeff+^ε2(u)du=0, (44)
 [1−ϕ(c,δM)]σ0−σeff2σeff+σ0+cσ1−σeff2σeff+σ1 +δM∫0∂ϕ(c,u)∂uσ2(u)−σeff2σeff+σ2(u)du=0. (45)

Based upon the volume averaging procedure, Eqs. (40) and (42) were first proposed in Sushko2013 (), and Eqs. (44) and (45) in Sushko2018polymer (). Here, Eqs. (42) and (45) are finally substantiated with a statistical mechanics formalism and an internally closed homogenization procedure.

## Viii Comparison with analytical results and numerical simulations

Before we proceed to processing experimental data with the above theory, it is of crucial importance to test its validity by contrasting it with other authors’ analytical and computer simulation results.

For 3D systems of particles with hard cores and fully penetrable shells calculations Rikvold85 (), done within the scaled-particle approximation Reiss1959 () for hard-sphere fluids, give

 ϕ(c,δ)=1−(1−c)exp[−((1+δ)3−1)c1−c] ×exp{−3(1+δ)3c22(1−c)3[2−31+δ+1(1+δ)3 −(31+δ−6(1+δ)2+3(1+δ)3)c]}. (46)

This result is in very good agreement with Monte Carlo simulations Lee1988 () (see also Rottereau2003 ()). Simulation results for (and ) of 3D systems of monosized particles are available in Siekierski2005 (); Siekierski2007 (); Siekierski2006 (), where was calculated using the RRN approach Kirkpatrick1971 (); Kirkpatrick1973 ().

In simulations Siekierski2005 (); Siekierski2007 (); Siekierski2006 (), the virtual RRN samples, representing the morphology and phase structure of simulated dispersions, were built from a matrix with cubic cells by placing spherical grains randomly into the matrix and using a special algorithm for avoiding conflicts of spatial restrictions on the grain locations. Each cell was marked as belonging to a grain if the center of the cell lay inside it. The procedure was repeated until the assumed volume fraction of the filler was attained. The residual part of the virtual sample (not belonging to the grains) was attributed to either the shells, with prescribed thickness , or the matrix (all cells belonging to neither the grains nor the shell). Finally, each cell was represented as a parallel combination of a resistor and a capacitor, with their parameters taken according to the assumed material parameters of all the phases; the electrical parameters used are summarized in Table 1. Replacing each pair of neighboring cells by an equivalent electrical circuit with the corresponding impedance and the nodes at the centers of the cells, the virtual samples were analyzed as 3D networks of such impedances.

We use results Rikvold85 (); Siekierski2005 (); Siekierski2007 (); Siekierski2006 () to test the validity of Eqs. (42) and (45) in the following five steps.

### viii.1 Mapping the geometrical parameters in models Rikvold85 (); Siekierski2007 () onto each other

This step is to contrast results Rikvold85 (); Siekierski2007 () for in order to find the relations between the parameters and for a dispersion of spherical core-shell particles Rikvold85 () and their counterparts and in simulations Siekierski2005 (); Siekierski2007 (); Siekierski2006 ().

Evidently, for a given and under the condition that the assumed volume fraction of the filler is attained, , the simulation procedure Siekierski2007 () leads to values of different from those of . Indeed, suppose that spherical grains of radius , where is the edge of the cubic cell, are used to generate a virtual sample of volume . Then and , for only one cell can belong to each grain. To achieve this filler concentration in the same-volume dispersion of spherical particles with shell thickness , their core radius must be equal to . For such particles, the relative shell thickness . Correspondingly,

 δ=Kδ′, (47)

where, in our example, . Considering as a fitting parameter, one can generalize Eq. (47) to the situations where each grain contains a large number of cells. The greater this number is, the closer to unity is expected to be. However, for , used in Siekierski2005 (); Siekierski2007 (); Siekierski2006 (), the deviation of from unity should be noticeable.

The results of applying Eqs. (46) and (47) to simulation data Siekierski2007 () for the composition of dispersions of spherical hard-core–penetrable-shell particles at different values of are presented in Figs. 2 and 3. They clearly demonstrate that under the proper choice of , these equations describe data Siekierski2007 () very well; the found values of turn out to be close to our above estimates.

### viii.2 Verifying functional relationship (42) between σeff and ϕ

With this object in view, the dispersion composition is assumed to be known for different values of at fixed and . Taking the corresponding values of from simulations Siekierski2007 () (Figs. 2 and 3), we then use Eq. (42) to calculate as a function of for given and without referring to Eqs. (46) and (47).

The results so obtained are shown in Fig. 4, together with conductivity simulation results Siekierski2007 (). If , the agreement between both theories is good for all three sets of data (, 7, and at fixed ). At lower values of , our theory predicts the percolation-type behavior of (see also Fig. 5 and the inset), with the threshold concentration that can be estimated from the relation Sushko2013 (). For the indicated sets of data, the estimations with Eqs. (42) and (47) give, respectively, (, ), (, ), and 0.046 (, ). Contrastingly, the simulated values of seem to increase gradually even at the lowest values of , and the percolation thresholds, if any, are hard to detect. This situation is typical of conductivity simulations for finite-size systems where the percolation threshold is a random non-Gaussian variable Berlyand1997 ().

### viii.3 Testing our model for the case of uniform penetrable shells

This step consists in fitting conductivity data Siekierski2005 (); Siekierski2007 () using Eq. (42) with given by Eq. (46) and given by Eq. (47). As Fig. 5 demonstrates, the value of , determined by fitting the composition data Siekierski2007 () (step A), is also appropriate to reproduce conductivity data Siekierski2007 () (this step). Similarly, Figs. 6 and 7 clearly indicate that the parameter alone, with a reasonable fitting value for each series, is sufficient to reproduce all ten series of simulation data Siekierski2007 () for of dispersions of particles with uniform penetrable shells. This fact is a strong argument in favor of the model expressed by Eqs. (42) and (46).

### viii.4 Scrutinizing the conductivity maximum positions

Under the condition , typical of simulations Siekierski2005 (); Siekierski2007 (), Eq. (42) can be greatly simplified by passing to the limit where it takes the form

 4σ3eff−2[(2−3ϕ)σ0−(1+3c−3ϕ)σ2]σ2eff −(2−3c)σ0σ2σeff=0. (48)

A nontrivial physically meaningful solution to Eq. (48) is

 σeff=34(A+√B+A2), (49)

where

 A≡(23−ϕ)σ0+(ϕ−c−13)σ2, (49a) B≡43(23−c)σ0σ2. (49b)

For the data series in Figs. 6 and 7, the versus plots given by Eqs. (49) and Eq. (42) are indistinguishable.

The concentrations where the conductivity maxima occur are found from the conditions and . Since near these maxima , it follows from Eq. (48) and the first condition that

 ∂ϕ(c,δ)/∂c|c=cmax=1, (50)

and that the derivatives and have the same sign at . According to Eq. (46), for . So, for such a , the second condition is fulfilled, and has a local maximum at indeed. Its value is given by Eqs. (49) at