1 Introduction

Accessible and applicable characterizations of the irregular geometry and connectivity of disordered media are currently enjoying a surge of attention [1,2,3,4,5,6,7,8,9,10,11]. Disordered systems as diverse as nano-structured composites [12,13,14], randomly doped crystals [15], metallic alloys [16], sandstones [17], sedimentary rocks [18, 19], chalk [20], wood [21], soil [22], zeolites [23], polymer blends [24], porous silicon [25], or ceramics [26] all require quantitative characterizations of their connectivity as a prerequisite for predicting the transport of charge, volume, mass, momentum or energy in these materials.

Much attention was paid in recent years [1, 11] to the study of certain additive quantities [27] known as intrinsic volumes, quermaßintegrals [28] (p.240), valuations [29] (p.173) or Minkowski functionals [30, 31]. Additivity means that adding the values for two samples equals the sum of the value for the union plus that for the intersection [27], Eq. (42). In particular the additive Euler characteristic [27] (p.31), [28] (p.237) is often believed to determine transport in disordered media [32, 33] or to “account for connectivity in a macroscale way” [1](p.325). On the other hand it is widely appreciated, that transport in random systems depends upon the existence of percolating paths from inlet to outlet [34,35,36,37,38]. Rigorous relations between the Euler characteristic and percolation were established in Eqs. (48)–(52) in Ref. [39], p. 386. Except for certain special geometries these relations show that in general the Euler characteristic cannot possibly suffice to characterize connectivity or transport. Most importantly, the Euler characteristic is additive, while connectivity, connectedness or “percolativity” are not.

Disordered systems often exhibit anisotropic geometry and/or connectivity. Examples include metamaterials [40], composites [41,42,43] metallic filters [44], shales and clay [45, 46], wood [21, 47], layered porous media [48, 49], fibrous media [50, 51], piezoelectric polymer foams, photonic crystals [52] or liquid crystals. It is therefore important to generalize theories and geometric characterizations to the case of anisotropic media.

Geometric quantities for anisotropic media, especially vectorial [53] and tensorial [54, 55] extensions of additive functionals are currently under active investigation [56]. Let us emphasize that higher order Minkowski functionals as well as their tensorial extensions are studied primarily as descriptors of geometry and have, to the best of our knowledge, not been used as predictors for physical transport. Our intention in this work is to introduce novel effective medium approximations to predict percolation and conduction in disordered materials. Recall that classical effective medium theories [35, 57,58,59,60] contain only volume, more precisely volume density i.e. volume fraction, as a geometrical descriptor. In self-consistent effective medium theories for mixtures with infinite material contrast a percolation transition emerges at volume fraction 1/3. Anisotropic variants of effective medium theory based on volume fraction [41, 61, 62] exhibit the same percolation transition and they continue to attract much interest [42, 52, 63]. Minkowski functionals other than volume, however, have, as far as we know, not been used to generalize effective medium theories to anisotropic media, and the same applies to pair connectivity [38], the pair connectedness [64], or n-connectedness [65].

The self-consistent effective medium approximation for the effective dc-conductivity of homogeneous and isotropic media with porosity \({\overline{\phi }}\) reads [35, 57, 58]

$$\begin{aligned} {\overline{\phi }}\;\frac{(\upsigma _\mathbbm {p}^\prime -{\overline{\upsigma }}^\prime )}{\upsigma _\mathbbm {p}^\prime +2{\overline{\upsigma }}^\prime }+ (1-{\overline{\phi }})\frac{(\upsigma _\mathbbm {m}^\prime -{\overline{\upsigma }}^\prime )}{\upsigma _\mathbbm {m}^\prime +2{\overline{\upsigma }}^\prime }= 0 \end{aligned}$$
(1)

where \(\upsigma _\mathbbm {p}^\prime \) and \(\upsigma _\mathbbm {m}^\prime \) are the dc-conductivities of two isotropic materials filling pores \(\mathbbm {p}\) and matrix \(\mathbbm {m}\) (they will be introduced more precisely in Eq. (7) below). For infinite contrast \(\upsigma _\mathbbm {p}^\prime /\upsigma _\mathbbm {m}^\prime =\infty \) this equation has a singularity at \({\overline{\phi }}={\overline{\phi }}_c=1/3\) which marks the percolation transition. Nonzero solutions for \({\overline{\upsigma }}^\prime \) exist only for \({\overline{\phi }}>{\overline{\phi }}_c\), while for \({{\overline{\phi }}\le 1/3}\) the effective conductivity \({\overline{\upsigma }}^\prime =0\) vanishes.

Equation (1) has been generalized to anisotropic media in the special case where \(\mathbbm {p}\) is a homogeneous mixture of aligned oblate ellipsoids \((x/b)^2+(y/b)^2+(z/c)^2\le 1\) with half axes \(b\ge c>0\). In this case the effective conductivity \({\overline{\upsigma }}^\prime \) becomes tensorial. Its eigenvalues \({\overline{\upsigma }}^\prime _x,{\overline{\upsigma }}^\prime _y,{\overline{\upsigma }}^\prime _z\) obey \({\overline{\upsigma }}^\prime _x={\overline{\upsigma }}^\prime _y\) and Eq. (1) becomes essentially two coupled equations [41], Eq. (8)

$$\begin{aligned} 0=&\frac{{\overline{\phi }}(\upsigma _\mathbbm {p}^\prime -{\overline{\upsigma }}^\prime _x)}{{\overline{\upsigma }}^\prime _x+N_a(\upsigma _\mathbbm {p}^\prime -{\overline{\upsigma }}^\prime _x)}+ \frac{(1-{\overline{\phi }})(\upsigma _\mathbbm {m}^\prime -{\overline{\upsigma }}^\prime _x)}{{\overline{\upsigma }}^\prime _x+N_a(\upsigma _\mathbbm {m}^\prime -{\overline{\upsigma }}^\prime _x)} \end{aligned}$$
(2a)
$$\begin{aligned} 0=&\frac{{\overline{\phi }}(\upsigma _\mathbbm {p}^\prime -{\overline{\upsigma }}^\prime _z)}{{\overline{\upsigma }}^\prime _z+N_c(\upsigma _\mathbbm {p}^\prime -{\overline{\upsigma }}^\prime _z)}+ \frac{(1-{\overline{\phi }})(\upsigma _\mathbbm {m}^\prime -{\overline{\upsigma }}^\prime _z)}{{\overline{\upsigma }}^\prime _z+N_c(\upsigma _\mathbbm {m}^\prime -{\overline{\upsigma }}^\prime _z)} \end{aligned}$$
(2b)

for \({\overline{\upsigma }}^\prime _x={\overline{\upsigma }}^\prime _y\) and \({\overline{\upsigma }}^\prime _z\). The depolarization factors \(N_a,N_c\) depend on a combined ratio R of the two unknowns \({\overline{\upsigma }}^\prime _x,{\overline{\upsigma }}^\prime _z\) and the two half axes bc through the relations

$$\begin{aligned} 1=&2N_a+N_c \end{aligned}$$
(2c)
$$\begin{aligned} N_c =&R^{-3}(1+R^2)(1-\arctan R) \end{aligned}$$
(2d)
$$\begin{aligned} R=&\sqrt{\frac{{\overline{\upsigma }}^\prime _z}{{\overline{\upsigma }}^\prime _x}\frac{b^2}{c^2}-1} \end{aligned}$$
(2e)

that couple the two Eqs. (2a) and (2b). As in the isotropic case, these two coupled equations exhibit a percolation singularity at \({\overline{\phi }}={\overline{\phi }}_c=1/3\).

Local porosity theory [36] is a generalized self-consistent effective medium theory that contains local porosity distributions and local percolation probabilities above and beyond volume fractions to characterize the geometry and connectedness of a disordered medium. An important consequence are percolation thresholds \({\overline{\phi }}_c\ne 1/3\). Predictions from local porosity theory agree with observations such as the empirically known Archie’s law [36]. Local porosity distributions and local percolation probabilities have become standard tools for image analysis of porous media [66, 67]. They were applied to clay rocks [68], fractured coal [69], Opalinus clay [70, 71], chalk [72], reticulate porous ceramics [73], macroporous reticulated silicon oxycarbide [74], dielectric characterization of wood [21], highly porous polymeric matrices [75], metallic filters [44], beryllium pebbles and arctic firns [64], and recently even to Ibuprofen tablets in pharmaceutical research [76]. Numerous publications [17, 32, 77,78,79,80,81] have demonstrated the predictivity of percolation probabilities and the total fraction of percolating cells for physical transport.

2 Problem and objective

A shortcoming of local porosity theory for electrical transport in heterogeneous media (introduced in Ref. [36] and reviewed in Refs. [37, 82]) is that the local percolation probability \(\uplambda \), which characterizes the geometrical connectivity of the medium, is a scalar but not a tensorial quantity. Disordered or heterogeneous media, however, generally have tensorial transport properties even if the constituent materials are isotropic.

An attempt to extend LPT  to anisotropic systems was made in Ref. [83]. The anisotropy is described using eigenvalues of the “orientation matrix”. The orientation matrix is constructed from the position vectors of those voxels within a measurement cell, that belong to the largest percolating cluster in the cell. The position vectors are computed relative to the “cluster center” defined as an arithmetic average. Two eigenvalue ratios of the orientation matrix are collected into local anisotropy distribution. The relation of the local anisotropy distribution to local percolation probabilities is not discussed in Ref. [83]. As a consequence the generalized effective medium equation is unchanged and transport quantities such as conductivities or permeabilities remain scalar quantities.

The objective of this work is to extend local porosity theory and effective medium theories to anisotropic media in such a way, that the transport quantities become tensorial. For brevity and simplicity we consider only transport of charge, i.e. the case of electrical transport.

3 Model

Let \(\mathbbm {p},\mathbbm {m}\subset \mathbbm {R}^d\) denote the domains (pore space, matrix) occupied by the two components of a heterogeneous or porous sample \(\mathbbm {S}=\mathbbm {p}\cup \mathbbm {m}\subset \mathbbm {R}^d\). The Maxwell equations [84] read

(3a)
(3b)
(3c)
(3d)

where \({\varvec{\mathrm {\mathrm {D}}}}({\varvec{\mathrm {r}}},t)\) is the electric displacement, \({\varvec{\mathrm {\mathrm {B}}}}({\varvec{\mathrm {r}}},t)\) is the magnetic induction, \({\varvec{\mathrm {\mathrm {E}}}}({\varvec{\mathrm {r}}},t)\) is the electric field, and \({\varvec{\mathrm {\mathrm {H}}}}({\varvec{\mathrm {r}}},t)\) is the magnetic field. Taking the divergence of Eq. (3d) gives the equation of continuity

(3e)

It relates the macroscopic charge density \(Q({\varvec{\mathrm {r}}},t)\) and the macroscopic current density \({\varvec{\mathrm {\mathrm {J}}}}({\varvec{\mathrm {r}}},t)\).

Causality and locality in space restrict the constitutive equations to convolutions in time written as

$$\begin{aligned}&{\varvec{\mathrm {\mathrm {J}}}}_c({\varvec{\mathrm {r}}},t)=\int \limits _{-\infty }^t {\varvec{\mathsf {\upsigma }}}^\prime ({\varvec{\mathrm {r}}},t-t^\prime ){\varvec{\mathrm {\mathrm {E}}}}({\varvec{\mathrm {r}}},t^\prime )\;{\mathrm d}t^\prime \end{aligned}$$
(4a)
$$\begin{aligned}&{\varvec{\mathrm {\mathrm {D}}}}({\varvec{\mathrm {r}}},t)=\upepsilon _0\int \limits _{-\infty }^t {\varvec{\mathsf {\upepsilon }}}_r({\varvec{\mathrm {r}}},t-t^\prime ){\varvec{\mathrm {\mathrm {E}}}}({\varvec{\mathrm {r}}},t^\prime )\;{\mathrm d}t^\prime \end{aligned}$$
(4b)
$$\begin{aligned}&{\varvec{\mathrm {\mathrm {H}}}}({\varvec{\mathrm {r}}},t)=\frac{1}{\upmu _0}\int \limits _{-\infty }^t \left[ {\varvec{\mathsf {\upmu }}}_r({\varvec{\mathrm {r}}},t-t^\prime )\right] ^{-1} {\varvec{\mathrm {\mathrm {B}}}}({\varvec{\mathrm {r}}},t^\prime ) \;{\mathrm d}t^\prime \end{aligned}$$
(4c)

for \({\varvec{\mathrm {r}}}\in \mathbbm {p}\cup \mathbbm {m}\), where \({\varvec{\mathsf {\upsigma }}}^\prime \) is the electrical conductivity tensor. The symbols \(\upepsilon _0,{\varvec{\mathsf {\upepsilon }}}_r\) denote the of dielectric permittivity of the vacuum resp. the permittivity tensor of the medium, and \(\upmu _0,{\varvec{\mathsf {\upmu }}}_r\) are the scalar resp. tensorial magnetic permeabilities of the vacuum resp. the medium with \(\upmu _0=4\pi \times 10^{-7}\)N/A\(^2\) and \(\upepsilon _0=1/(\upmu _0 c^2)\approx 8.8542\times 10^{-12}\)F/m. The conduction current due to free charges \({\varvec{\mathrm {\mathrm {J}}}}_c\) appearing in Ohm’s law, Eq. (4a), is related to the total current via the current of bound charges defined as \({\varvec{\mathrm {\mathrm {J}}}}_b({\varvec{\mathrm {r}}},t)={\varvec{\mathrm {\mathrm {J}}}}({\varvec{\mathrm {r}}},t)-{\varvec{\mathrm {\mathrm {J}}}}_c({\varvec{\mathrm {r}}},t)\). Fourier transformation with respect to t gives

(5a)
(5b)
(5c)

where \(\omega =2\pi \nu \) is the angular frequency and \(\nu \) denotes frequency. One has

$$\begin{aligned} {\varvec{\mathsf {\upepsilon }}}({\varvec{\mathrm {r}}},\omega )=\,&{\varvec{\mathsf {\upepsilon }}}^\prime ({\varvec{\mathrm {r}}},\omega ) +\mathrm {i}{\varvec{\mathsf {\upepsilon }}}^{\prime \prime }({\varvec{\mathrm {r}}},\omega )\nonumber \\ =\,&{\varvec{\mathsf {\upepsilon }}}_r({\varvec{\mathrm {r}}},\omega )+\mathrm {i}\frac{{\varvec{\mathsf {\upsigma }}}^{\prime }({\varvec{\mathrm {r}}},\omega )}{\omega } \end{aligned}$$
(6a)
$$\begin{aligned} {\varvec{\mathsf {\upsigma }}}({\varvec{\mathrm {r}}},\omega )=\,&{\varvec{\mathsf {\upsigma }}}^\prime ({\varvec{\mathrm {r}}},\omega ) +\mathrm {i}{\varvec{\mathsf {\upsigma }}}^{\prime \prime }({\varvec{\mathrm {r}}},\omega )\nonumber \\ =\,&{\varvec{\mathsf {\upsigma }}}^\prime ({\varvec{\mathrm {r}}},\omega )+ \mathrm {i}\omega \left( {\varvec{{\mathsf {1}}}}-{\varvec{\mathsf {\upepsilon }}}_r({\varvec{\mathrm {r}}},\omega )\right) \end{aligned}$$
(6b)

for the complex conductivity and dielectric tensors.

The material parameter functions for a composite material consisting of two homogeneous but anisotropic components are assumed to have the form

(7a)
(7b)
(7c)

where

$$\begin{aligned} \chi _{_\mathbbm {X}}({\varvec{y}})= {\left\{ \begin{array}{ll} 1 ,&{} {\varvec{y}}\in \mathbbm {X}\\ 0 ,&{} {\varvec{y}}\notin \mathbbm {X} \end{array}\right. } \end{aligned}$$
(8)

is the characteristic (or indicator) function of a subset \(\mathbbm {X}\subset \mathbbm {R}^d\). The relaxation frequencies, defined as

(9)

are assumed to exist. They can be used to make \(\omega \) dimensionless.

The discontinuity of the material parameters at the interface requires one to specify boundary conditions at \(\partial \mathbbm {p}=\partial \mathbbm {m}\). Mathematically the equations are interpreted in a weak sense as equations for distributions. Depending on the boundary conditions a suitable domain could be a Sobolev space for the potentials resp. a space of potential fields for the electromagnetic fields [85, 86]. The boundary conditions at the interface are

(10a)
(10b)
(10c)
(10d)

for \({\varvec{\mathrm {r}}}\in \partial \mathbbm {p}=\partial \mathbbm {m}\), where \(Q_{\partial }({\varvec{\mathrm {r}}},\omega )\) (resp. \({\varvec{\mathrm {\mathrm {J}}}}_{\partial }({\varvec{\mathrm {r}}},\omega )\)) are the Fourier transforms of (possibly time dependent) surface charge (resp. surface current) densities with support in \(\partial \mathbbm {p}=\partial \mathbbm {m}\). The notation \({\varvec{\mathrm {\mathrm {B}}}}_{\partial \mathbbm {p}}({\varvec{\mathrm {r}}},\omega )\) (resp. \({\varvec{\mathrm {\mathrm {B}}}}_{\partial \mathbbm {m}}({\varvec{\mathrm {r}}},\omega )\)) is the limiting value of the vector field as the point \({\varvec{\mathrm {r}}}\in \partial \mathbbm {p}=\partial \mathbbm {m}\) is approached from within \(\mathbbm {p}\) (resp. from within \(\mathbbm {m}\)).

It is assumed that there are neither volume nor surface charges or currents inside the medium so that \(Q=0\), \(Q_{\partial }=0\), \({\varvec{\mathrm {\mathrm {J}}}}=0\) and \({\varvec{\mathrm {\mathrm {J}}}}_{\partial }=0\) from now on. This assumption is commonly made to exclude phenomena and effects arising from electrical double layers and other physico-chemical processes in the interfacial region.

When the wavelength is large compared to the scale of heterogeneities the time derivatives in in Faraday’s law, Eq. (3), are small compared to spatial derivatives. Setting them to zero results in the quasistatic approximation valid for small frequencies

$$\begin{aligned} \omega \ll \frac{1}{\ell \sqrt{\upepsilon _0\upmu _0}} =\frac{c}{\ell }=\frac{299792458}{\ell }\mathrm {Hz} \end{aligned}$$
(11)

where \(\ell \) denotes the scale of heterogeneities and c is the speed of light. For heterogeneities with \(\ell \approx 10...100\mu \)m one finds \(\omega \ll 3...30\)THz approaching the infrared. The quasistatic approximation decouples the full equations (3). Fourier transformation of eqs. (3) with respect to t results in the quasistatic equations

(12a)
(12b)
(12c)
(12d)
(12e)

where \({\varvec{\mathrm {r}}}\in \mathbbm {p}\cup \mathbbm {m}\) in Eqs. (12a)–(12c), \({\varvec{\mathrm {r}}}\in \partial \mathbbm {p}=\partial \mathbbm {m}\) in Eqs. (12d)–(12e), and where \(\upepsilon ({\varvec{\mathrm {r}}},\omega )\) denotes the complex frequency dependent local dielectric function.

The effective macroscopic dielectric tensor \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) of self-consistent effective medium approximations [58] is usually defined by averaging Eq. (12c)

$$\begin{aligned} \langle {\varvec{\mathrm {\mathrm {D}}}}({\varvec{\mathrm {r}}},\omega )\rangle =\upepsilon _0\,{\overline{{\varvec{\mathsf {\upepsilon }}}}}(\omega ;{\varvec{\mathrm {\Pi }}})\,\langle {\varvec{\mathrm {\mathrm {E}}}}({\varvec{\mathrm {r}}},\omega )\rangle \end{aligned}$$
(13)

where the angular brackets \(\langle \cdot \rangle \) denote (spatial or ensemble) averaging and \({\varvec{\mathrm {\Pi }}}\) are parameters of physical importance (such as porosity or connectivity) that arise from averaging the smaller scale heterogeneities.

4 Local geometry

4.1 Disordered media

A disordered (or porous) medium \({\mathbbm {m}}\) (or \({\mathbbm {p}}\)) is assumed to be a random set. A random set is a set-valued random variable [87]. The probability space of random disordered geometries is denoted \((\Omega ,{\mathcal {A}},P)\). The set \(\Omega \) is the set of all closed subsets of \({\mathbbm {R}^d}\) including \(\emptyset \). The \(\sigma \)-algebra is generated by the system of closed sets having nonempty intersection with compact sets and the random geometry distribution P is the image measure of the measure on the probability space underlying the set-valued random variable. For mathematical details on random sets the interested reader is referred to Ref. [82].

4.2 Local porosity

Let \({\mathbbm {K}\subset \mathbbm {R}^d}\) be a convex and compact set with centroid at the origin \({{\varvec{\mathrm {0}}}\in \mathbbm {R}^d}\). The centroid is the geometrical center of \({\mathbbm {K}}\). Then

$$\begin{aligned} \mathbbm {K}({\varvec{\mathrm {r}}})={\varvec{\mathrm {r}}}+\mathbbm {K}=\{{\varvec{\mathrm {r}}}+{{\varvec{\mathrm {q}}}}\in \mathbbm {R}^d:{{\varvec{\mathrm {q}}}}\in \mathbbm {K}\} \end{aligned}$$
(14)

denotes its translate by a vector \({{\varvec{\mathrm {r}}}\in \mathbbm {R}^d}\). The local porosity in a measurement cell \({\mathbbm {K}({\varvec{\mathrm {r}}})}\) placed at position \({\varvec{\mathrm {r}}}\) is defined as

$$\begin{aligned} \phi (\mathbbm {K}({\varvec{\mathrm {r}}}))=\frac{|\mathbbm {p}\cap \mathbbm {K}({\varvec{\mathrm {r}}})|}{|\mathbbm {K}({\varvec{\mathrm {r}}})|} \end{aligned}$$
(15)

where

$$\begin{aligned} |\mathbbm {K}|=\int \limits _\mathbbm {K}\;{\mathrm d}^d{{\varvec{\mathrm {q}}}}= \int \limits _{\mathbbm {R}^d}\chi _{\mathbbm {K}}({{\varvec{\mathrm {q}}}})\;{\mathrm d}^d{{\varvec{\mathrm {q}}}} \end{aligned}$$
(16)

is the volume of a set \({\mathbbm {K}}\). The local matrix volume fraction is \({\phi _{\mathbbm {m}}=1-\phi }\).

4.3 Anisotropic local matrix fraction

Consider from now on \(d=3\) and let \({\mathbbm {m}_i}\) denote the local matrix space clusters defined as (path-)connectedness components of \({\mathbbm {m}\cap \mathbbm {K}({{\varvec{\mathrm {r}}}})}\). Then

$$\begin{aligned} \mathbbm {m}\cap \mathbbm {K}({{\varvec{\mathrm {r}}}})=\bigcup _{i=1}^{n_\mathbbm {m}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\mathbbm {m}_i \end{aligned}$$
(17)

where \(\mathbbm {m}_i\cap \mathbbm {m}_j=\emptyset \) for \(i\ne j\) and \(n_\mathbbm {m}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))\) is the total number of local clusters within \({\mathbbm {K}({{\varvec{\mathrm {r}}}})}\). The barycenters of these local clusters are

$$\begin{aligned} {{\varvec{\mathrm {r}}}_\mathbbm {m}}_i&= \frac{1}{|\mathbbm {m}_i|}\int \limits _{\mathbbm {m}_i}{{\varvec{\mathrm {q}}}}\;{\mathrm d}^3{{\varvec{\mathrm {q}}}} \end{aligned}$$
(18)

and

$$\begin{aligned} {\phi _{\mathbbm {m}}}_i(\mathbbm {K}({{\varvec{\mathrm {r}}}}))=\frac{|\mathbbm {m}_i|}{|\mathbbm {K}({{\varvec{\mathrm {r}}}})|} \end{aligned}$$
(19)

is their local volume fraction. The matrix is assumed to have a constant density taken to be unity. Averaging over the local clusters the tensor

(20a)

is introduced, where \({\varvec{{\mathsf {X}}}}({{\varvec{\mathrm {r}}}})={{\varvec{\mathrm {r}}}}\otimes {{\varvec{\mathrm {r}}}}\), is the tensor product. The tensor is real, symmetric and non-negative definite. In the following the tensor is normalized as

$$\begin{aligned} {{\varvec{{\mathsf {J}}}}_\mathbbm {m}}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))= \frac{{\varvec{{\mathsf {I}}}}_\mathbbm {m}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}{\mathrm {Tr}\,{\varvec{{\mathsf {I}}}}_\mathbbm {m}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))} \end{aligned}$$
(20b)

by dividing with its trace. The non-negative eigenvalues \({J_{\mathbbm {m}a}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\), \({J_{\mathbbm {m}b}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\), \({J_{\mathbbm {m}c}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\) obey \(J_{\mathbbm {m}a}(\mathbbm {K}({{\varvec{\mathrm {r}}}})) +J_{\mathbbm {m}b}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))+J_{\mathbbm {m}c}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))=1\), and are assumed to be ordered as \({1\ge J_{\mathbbm {m}a}\ge J_{\mathbbm {m}b}\ge J_{\mathbbm {m}c}\ge 0}\).

4.4 Local percolation

Let \({\partial \mathbbm {K}({{\varvec{\mathrm {r}}}})}\) denote the boundary of a measurement cell with centroid \({{\varvec{\mathrm {r}}}}\). A measurement cell \({\mathbbm {K}({{\varvec{\mathrm {r}}}})={{\varvec{\mathrm {r}}}}+\mathbbm {K}}\) is called percolating if there exists a continuous path

$$\begin{aligned} {\varvec{\mathrm {p}}}:[0,1]&\rightarrow \mathbbm {R}^3\nonumber \\ t&\mapsto {\varvec{\mathrm {p}}}(t)={\varvec{\mathrm {p}}}_t \end{aligned}$$
(21a)

with

(21b)
(21c)
(21d)

where \({\mathbbm {B}({\varvec{\mathrm {x}}},a)=\{{{\varvec{\mathrm {r}}}}\in \mathbbm {R}^3\!:\!|{{\varvec{\mathrm {r}}}}-{\varvec{\mathrm {x}}}|\le a\}}\) denotes a ball of radius a centered at \({\varvec{\mathrm {x}}}\). Equation (21) is illustrated in Figure 1.

Fig. 1
figure 1

Schematic illustration from Ref.[88] of the local percolation criterion (21). The grey hatched circle centered at \({\varvec{\mathrm {p}}}_0\) is the forbidden ball \({\mathbbm {B}\left( {\varvec{\mathrm {p}}}_0,|{{\varvec{\mathrm {r}}}}-{\varvec{\mathrm {p}}}_0|\right) }\). Two paths are shown that start at \({\varvec{\mathrm {p}}}_0\) on the left boundary. The dashed path is percolating according to the percolation criterion (21), because it ends outside the forbidden ball. The dash-dotted path is not percolating, because it ends inside the forbidden ball

The set of local percolating directions for a measurement cell \({\mathbbm {K}({{\varvec{\mathrm {r}}}})}\) centered at \({{\varvec{\mathrm {r}}}}\) is defined as

$$\begin{aligned} \mathbbm {L}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))=&\left\{ \!\frac{{\varvec{\mathrm {p}}}_0\!-{\varvec{\mathrm {p}}}_1}{|{\varvec{\mathrm {p}}}_0\!-{\varvec{\mathrm {p}}}_1|}\!\in \!\mathbbm {S}^2\!:\! {\varvec{\mathrm {p}}}_0,{\varvec{\mathrm {p}}}_1\!\in \!\partial \mathbbm {K}({{\varvec{\mathrm {r}}}}),\right. \nonumber \\ {}&\left. {\varvec{\mathrm {p}}}_0\!\rightsquigarrow \!{\varvec{\mathrm {p}}}_1\!\right\} \end{aligned}$$
(22)

where \({\mathbbm {S}^2\!=\!\partial \mathbbm {B}({\varvec{\mathrm {0}}},1)}\) is the unit sphere and \({\varvec{\mathrm {p}}}_0\rightsquigarrow {\varvec{\mathrm {p}}}_1\) means that \({\varvec{\mathrm {p}}}_0\) is the starting point and \({\varvec{\mathrm {p}}}_1\) is the end point of a percolation path in the sense of Eq. (21). The local percolation indicators, defined as

$$\begin{aligned} {\Lambda (\mathbbm {K}({\varvec{\mathrm {r}}}))= {\left\{ \begin{array}{ll} 0 &{} \text { if } \mathbbm {L}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))=\emptyset \\ 1 &{} \text {otherwise} , \end{array}\right. }} \end{aligned}$$
(23)

are a family of \(\{0,1\}\)-valued random variables

$$\begin{aligned} \Lambda (\;\cdot \;;\mathbbm {K}):(\Omega ,{\mathcal {A}},P)&\rightarrow \{0,1\}\nonumber \\ \mathbbm {p}&\mapsto \Lambda (\mathbbm {p};\mathbbm {K}) \end{aligned}$$
(24)

indexed by \({\mathbbm {K}}\). The image of the measure P under the map

$$\begin{aligned} \Lambda (P)(\;\cdot \;;\mathbbm {K})&:\{\emptyset ,\{0\},\{1\},\{0,1\}\}\rightarrow [0,1]\nonumber \\ \emptyset&\mapsto 0, \qquad \qquad \quad \forall \,\mathbbm {K} \nonumber \\ \{0\}&\mapsto 1-p(\mathbbm {K}), \nonumber \\ \{1\}&\mapsto p(\mathbbm {K}), \nonumber \\ \{0,1\}&\mapsto 1 \qquad \qquad \qquad \forall \,\mathbbm {K} \end{aligned}$$
(25)

defines the total percolation probability \({p(\mathbbm {K})}\) of the measurement cell \({\mathbbm {K}}\).

4.5 Anisotropic local percolation

The set \({\mathbbm {L}}\) of local percolating directions defines a tensor

(26a)

with respect to the cell center \({{\varvec{\mathrm {r}}}}\). It is again normalized

(26b)

by dividing with the trace. The non-negative eigenvalues \({J_{\mathbbm {p}a}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\), \({J_{\mathbbm {p}b}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\), \({J_{\mathbbm {p}c}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\) obey \({J_{\mathbbm {p}a}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\) \({+J_{\mathbbm {p}b}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))+J_{\mathbbm {p}c}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))=1}\), and are assumed to be ordered as \({1\ge J_{\mathbbm {p}a}\ge J_{\mathbbm {p}b}\ge J_{\mathbbm {p}c}\ge 0}\).

4.6 Local geometry distribution

The local geometry of a random set \({\mathbbm {p}}\) inside a measurement cell \({\mathbbm {K}}\) is characterized by the geometric parameters \({g(\mathbbm {p};\mathbbm {K})=(\Lambda ,\phi ,d,e)}\) where

$$\begin{aligned} d&=\Lambda d_\mathbbm {p}+(1-\Lambda )d_\mathbbm {m} \end{aligned}$$
(27a)
$$\begin{aligned} e&=\Lambda e_\mathbbm {p}+(1-\Lambda )e_\mathbbm {m} \end{aligned}$$
(27b)

and

(27c)
(27d)

parametrize the anisotropy of the measurement cell \({\mathbbm {K}}\) with local porosity \(\phi \). These quantities form a family of \({\mathbbm {G}}\)-valued random variables

(28)

indexed by \({\mathbbm {K}}\) where \({\mathbbm {G}=\{0,1\}\times [0,1]^3 }\) and \({\mathcal {B}}\) is the standard Borel \(\sigma \)-algebra of the unit hypercube in \({\mathbbm {R}^4}\). The family of image measures

$$\begin{aligned} \mu (\;\cdot \;;\mathbbm {K}):{\mathcal {B}}\rightarrow & {} [0,\infty ]\nonumber \\ B\mapsto & {} P(g^{-1}(B;\mathbbm {K})) \end{aligned}$$
(29)

of P under g is the family of local geometry distributions with family index \({\mathbbm {K}}\).

In practical applications the distribution P of the random geometries is unknown. Instead one is given a sample assumed to be a realization of a random geometry. The local geometry distribution is then estimated as the joint empirical probability measure

$$\begin{aligned} {\mu (g;\mathbbm {K})= \lim _{N\rightarrow \infty }\frac{1}{N}\sum ^N_{i=1} \delta \left( g-g(\mathbbm {K}({\varvec{\mathrm {r}}}_i))\right) } \end{aligned}$$

where \({\varvec{\mathrm {r}}}_i\) is a sequence of cell centers such that all \({\mathbbm {K}({\varvec{\mathrm {r}}}_i)}\) are mutually disjoint, and N denotes the total number of different placements of \({\mathbbm {K}}\). Mathematically this empirical density is an estimator for a Young measure (see Refs. [39, 82, 89] for details).

4.7 Limiting local geometry distributions

The dependence of \({\mu (g;\mathbbm {K})}\) on \({\mathbbm {K}}\) disappears in the limit where \({\mathbbm {K}}\) shrinks to a point or expands to become \({\mathbbm {R}^3}\). In the limit \({|\mathbbm {K}|\rightarrow 0}\) of small pointlike measurement cells one has with Eq. (27)

(30a)
(30b)

where \({{\overline{\phi }}=\phi (\mathbbm {S})}\) is the total porosity of the sample.

In the opposite \({|\mathbbm {K}|\rightarrow \infty }\) limit the limiting local geometry distributions are also expected to become independent [89] of \({\mathbbm {K}}\). Then Eq. (27) gives

(31a)
(31b)
(31c)

where \({\mu _1(\phi ,d_\mathbbm {p},e_\mathbbm {p})}\) is the conditional probability density that a cell is percolating in the scaling limit, has porosity \(\phi \) and percolation anisotropy parameters \({d_\mathbbm {p},e_\mathbbm {p}}\), while \({\mu _0(\phi ,d_\mathbbm {m},e_\mathbbm {m})}\) is the conditional probability density that a cell is non-percolating in the scaling limit, has porosity \(\phi \) and percolation anisotropy parameters \({d_\mathbbm {m},e_\mathbbm {m}}\). The quantity \(\uplambda (\phi )\) is the conditional probability that a very large measurement cell is percolating given that it has local porosity \(\phi \). It may be called limiting local percolation probability (cf. Ref. [36]) in the sense that

(31d)

is expected to hold from Eq. (25).

The possible macroscopic limits have been classified in Ref. [89], Sec.VII.b, into four cases. In the following the macroscopically heterogeneous cases will be neglected and the porous medium is assumed to be macroscopically homogeneous and anisotropic. If it is macroscopically non-random one has

(32a)
(32b)

where \({\overline{\phi }}\) is the macroscopic porosity and \({{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}}\) are the macroscopic anisotropy parameters for pore and matrix. The simple choice

$$\begin{aligned} \uplambda (\phi )&= {\left\{ \begin{array}{ll} 0, &{} \text { for }0\le \phi<{\overline{\phi }}\\ {\overline{\uplambda }}, &{} \text { for }\phi ={\overline{\phi }}\\ 1, &{} \text { for }{\overline{\phi }}<\phi \le 1 \end{array}\right. } \end{aligned}$$
(32c)

will be used below. Expectation values of measurable functions \({f:\mathbbm {G}\rightarrow \mathbbm {R}}\) are computed as

(33)

5 Coated ellipsoidal inclusion

5.1 Formulation of the local problem

Let \(d=3\) and let \({{\varvec{\mathrm {r}}}^{\!\mathrm {T}}=(x,y,z)\in \mathbbm {R}^3}\) denote a vector in cartesian coordinates. The potential problem for the unknown electrostatic potentials \({U_\mathbbm {c},U_\mathbbm {s},U_\mathbbm {u}}\) in the domains \({\mathbbm {c},\mathbbm {s},\mathbbm {u}\subset \mathbbm {R}^3}\) reads

(34a)
(34b)
(34c)
(34d)
(34e)
(34f)
(34g)
(34h)

where \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}}\) are the dielectric tensors of the homogeneous but anisotropic dielectric materials filling the domains \({\mathbbm {c},\mathbbm {s},\mathbbm {u}}\) and \({\varvec{\mathrm {\mathrm {E}}}}_{\mathrm {ext}}\) denotes a constant applied field. The union of \({\mathbbm {c}}\) and \({\mathbbm {s}}\) forms an ellipsoidal inclusion

(35a)

consisting of an ellipsoidal core

(35b)

surrounded by a shell

$$\begin{aligned} {\mathbbm {s}=\mathbbm {e}\setminus \mathbbm {c}} \end{aligned}$$
(35c)

that is embedded into an environment (or background)

$$\begin{aligned} {\mathbbm {u}=\mathbbm {R}^3\setminus \mathbbm {e}} \end{aligned}$$
(35d)

defined as the complement of \({\mathbbm {e}}\). It is assumed that the quadratic forms \({{\varvec{{\mathsf {Q}}}}_\mathbbm {e},{\varvec{{\mathsf {Q}}}}_\mathbbm {c}}\) are non-degenerate, and that \({\mathbbm {c}\subset \mathbbm {s}}\) and \({\mathbbm {s}\setminus \mathbbm {c}\ne \emptyset }\) holds throughout. The finite two-dimensional boundaries are denoted

(35e)
(35f)
(35g)
(35h)

and the vector fields \({{\varvec{\mathrm {\mathrm {n}}}}_{\partial \mathbbm {e}},{\varvec{\mathrm {\mathrm {n}}}}_{\partial \mathbbm {c}}}\) in Eq. (34) are unit normal vector fields on \({\mathbbm {e}}\) and \({\mathbbm {c}}\). The boundary conditions (34d), (34e) at \({\partial \mathbbm {u}\!\cap \!\partial \mathbbm {s}\!=\!\partial \mathbbm {e}}\) and (34f), (34g) at \({\partial \mathbbm {s}\cap \partial \mathbbm {c}\!=\!\partial \mathbbm {c}}\) express continuity of the potentials and the normal component of the dielectric displacement. The boundary condition (34h) demands that the potential vanishes at infinity \(|{{\varvec{\mathrm {r}}}}|\rightarrow \infty \).

The non-degenerate quadratic form matrices \({{\varvec{{\mathsf {Q}}}}_\mathbbm {e},{\varvec{{\mathsf {Q}}}}_\mathbbm {c}}\) (with \({\det {\varvec{{\mathsf {Q}}}}_\mathbbm {e}\ne 0, \det {\varvec{{\mathsf {Q}}}}_\mathbbm {c}\ne 0}\)) depend on the material properties \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}}\) of the shell as follows. Assume that \({\det {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}\ne 0}\) and define the coordinate transformation \({{{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}}\) by the formula

(36a)
(36b)

as the square root of \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}}\). The linear, non-orthogonal and non-singular coordinate transformation (36b) relates the vector \({{\varvec{\mathrm {r}}}}=(x,y,z)^{\!\mathrm {T}}\) in the original coordinate system with the vector \({{{{\varvec{\mathrm {r}}}}_\mathbbm {s}}=({x_\mathbbm {s}},{y_\mathbbm {s}},{z_\mathbbm {s}})^{\!\mathrm {T}}}\) in the transformed coordinates. To specify the transformation uniquely it is assumed that

(37a)
(37b)

with

(37c)
(37d)

holds true. Here \({a_{\mathbbm {e},\mathbbm {s}}, b_{\mathbbm {e},\mathbbm {s}}, c_{\mathbbm {e},\mathbbm {s}}}\) denote the semiaxes of the outer shell ellipsoid \({\mathbbm {e}}\) and \({a_{\mathbbm {c},\mathbbm {s}}, b_{\mathbbm {c},\mathbbm {s}}, c_{\mathbbm {c},\mathbbm {s}}}\) for the core ellipsoid \({\mathbbm {c}}\) in \({({x_\mathbbm {s}},{y_\mathbbm {s}},{z_\mathbbm {s}})}\)-coordinates.

Finally, it is assumed that the ellipsoids \({\mathbbm {e}_{\mathbbm {s}}={{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{-1}\mathbbm {e}}\) and \({\mathbbm {c}_{\mathbbm {s}}={{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{-1}\mathbbm {c}}\) are confocal, i.e. it is assumed there exists a number \({0<C<c_{\mathbbm {e},\mathbbm {s}}^2}\) such that

(38a)
(38b)
(38c)

holds. Then \({a_{\mathbbm {c},\mathbbm {s}}<a_{\mathbbm {e},\mathbbm {s}}}\), \({b_{\mathbbm {c},\mathbbm {s}}<b_{\mathbbm {e},\mathbbm {s}}}\) and \({c_{\mathbbm {c},\mathbbm {s}}<c_{\mathbbm {e},\mathbbm {s}}}\) and the coated ellipsoid is parameterized by four parameters. If the half axes \({a_{\mathbbm {e},\mathbbm {s}},b_{\mathbbm {e},\mathbbm {s}},c_{\mathbbm {e},\mathbbm {s}}}\) are normalized by division with their sum \({a_{\mathbbm {e},\mathbbm {s}}+b_{\mathbbm {e},\mathbbm {s}}+c_{\mathbbm {e},\mathbbm {s}}}\), then three parameters suffice to characterize the coated ellipsoid. They could be chosen as the volume fraction \({\phi _\mathbbm {c}}\), and two normalized half axes, or as \({\phi _\mathbbm {c}}\) and two of the axis ratios

(39a)
(39b)
(39c)

for shell and core. Here the axis ratio \({\varphi _{\mathbbm {c}}}\) is a core to shell length fraction for the intermediate axes. The confocality condition (38) implies the relations

(40a)
(40b)

between these axis ratios. The core \({\mathbbm {c}}\) within the coated ellipsoidal inclusion \({\mathbbm {e}}\) has a local volume fraction

(41)

that can be used to eliminate \({\varphi _{\mathbbm {c}}}\) in exchange for \({\phi _\mathbbm {c}}\). Then with \({\phi _\mathbbm {c}\rightarrow 1}\) for \(C\rightarrow 0\), while \({\phi _\mathbbm {c}\rightarrow 0}\) for \({C\rightarrow c_{\mathbbm {e},\mathbbm {s}}^2}\). Given \({d_\mathbbm {s},e_\mathbbm {s}}\) and \({\phi _\mathbbm {c}}\) these equations are solved for \({\varphi _{\mathbbm {c}}}\) to determine

(42a)
(42b)

as functions of \({d_\mathbbm {s},e_\mathbbm {s}}\) and \({\phi _\mathbbm {c}}\). The functions \(D_\mathbbm {c},E_\mathbbm {c}:[0,1]^3\rightarrow [0,1]\) show clearly, that \({d_\mathbbm {c},e_\mathbbm {c}}\) are not independent parameters. Instead, they are fixed as soon as \({\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s}}\) are given.

5.2 Solution of the local problem

To solve Eq. (34), one first transforms Eq. (34c) into the Laplace equation in vacuum using the transformation

(43a)
(43b)

where it is assumed that \({\det {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}\ne 0}\). To specify the linear, non-orthogonal and non-singular coordinate transformation \({{{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {u}}}\) uniquely it is assumed that

(44)

holds true in \(({x_\mathbbm {u}},{y_\mathbbm {u}},{z_\mathbbm {u}})\)-coordinates with \(a_{\mathbbm {e},\mathbbm {u}} \ge b_{\mathbbm {e},\mathbbm {u}} \ge c_{\mathbbm {e},\mathbbm {u}} > 0\).

Following Refs. [90] and [91], Eq. (44), the solution of the boundary value problem (34) can be written as

(45a)
(45b)
(45c)

where \({{\xi _\mathbbm {u}}={\xi _\mathbbm {u}}({{\varvec{\mathrm {r}}}})}\) (resp. \({{\xi _\mathbbm {s}}={\xi _\mathbbm {s}}({{\varvec{\mathrm {r}}}})}\)) is the first ellipsoidal coordinate of the vector \({{\varvec{\mathrm {r}}}}\) in \({({x_\mathbbm {u}},{y_\mathbbm {u}},{z_\mathbbm {u}})}\)-coordinates (resp. \({({x_\mathbbm {s}},{y_\mathbbm {s}},{z_\mathbbm {s}})}\)-coordinates) and \(C>0\) is the confocal constant in Eq. (38). The ellipsoidal coordinates \({({\xi _\mathbbm {u}},{\eta _\mathbbm {u}},{\zeta _\mathbbm {u}})}\) of a point \({({x_\mathbbm {u}},{y_\mathbbm {u}},{z_\mathbbm {u}})}\) are defined as the three roots \({\xi _\mathbbm {u}},{\eta _\mathbbm {u}},{\zeta _\mathbbm {u}}\) of the equation [92] (p. 19)

$$\begin{aligned}&\frac{{x_\mathbbm {u}}^2}{a_{\mathbbm {e},\mathbbm {u}}^2+ u} +\frac{{y_\mathbbm {u}}^2}{b_{\mathbbm {e},\mathbbm {u}}^2+ u} +\frac{{z_\mathbbm {u}}^2}{c_{\mathbbm {e},\mathbbm {u}}^2+ u} = 1 \end{aligned}$$
(46a)

which obey

$$\begin{aligned} -&a_{\mathbbm {e},\mathbbm {u}}^2 \!\le \! {\zeta _\mathbbm {u}} \!\le \! -b_{\mathbbm {e},\mathbbm {u}}^2 \!\le \! {\eta _\mathbbm {u}} \!\le \! -c_{\mathbbm {e},\mathbbm {u}}^2 \!\le \! {\xi _\mathbbm {u}} < \infty . \end{aligned}$$
(46b)

Similarly, the ellipsoidal coordinates \(({\xi _\mathbbm {s}},{\eta _\mathbbm {s}},{\zeta _\mathbbm {s}})\) of a point \(({x_\mathbbm {s}},{y_\mathbbm {s}},{z_\mathbbm {s}})\) are defined as the three roots of the equation

$$\begin{aligned}&\frac{{x_\mathbbm {s}}^2}{a_{\mathbbm {c},\mathbbm {s}}^2+u} +\frac{{y_\mathbbm {s}}^2}{b_{\mathbbm {c},\mathbbm {s}}^2+u} +\frac{{z_\mathbbm {s}}^2}{c_{\mathbbm {c},\mathbbm {s}}^2+u} = 1 \end{aligned}$$
(47a)

which obey

$$\begin{aligned} -&a_{\mathbbm {c},\mathbbm {s}}^2\!\le \!{\zeta _\mathbbm {s}} \!\le \! -b_{\mathbbm {c},\mathbbm {s}}^2 \!\le \! {\eta _\mathbbm {s}} \le -c_{\mathbbm {c},\mathbbm {s}}^2 \!\le \! {\xi _\mathbbm {s}}\!<\!\infty \! . \end{aligned}$$
(47b)

To formulate expressions for \({\varvec{{\mathsf {G}}}}_\mathbbm {c},{\varvec{{\mathsf {G}}}}_\mathbbm {s}({\xi _\mathbbm {s}}),{\varvec{{\mathsf {G}}}}_\mathbbm {u}({\xi _\mathbbm {u}})\) three diagonal tensors, \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}}({\xi _\mathbbm {u}})\), \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}({\xi _\mathbbm {s}})\) and \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}({\xi _\mathbbm {s}})\) of depolarizing factors are needed. All of them have unit trace, i.e. are normalized. The diagonal entries \(N_{\mathbbm {e},\mathbbm {u}}^{({x_\mathbbm {u}})}\), \(N_{\mathbbm {e},\mathbbm {u}}^{({y_\mathbbm {u}})}\), and \(N_{\mathbbm {e},\mathbbm {u}}^{({z_\mathbbm {u}})}\) of \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}}({\xi _\mathbbm {u}})\) are

$$\begin{aligned} N_{\mathbbm {e},\mathbbm {u}}^{(i)}({\xi _\mathbbm {u}})&=\frac{3|\mathbbm {e}|}{8\pi } \int \limits _{\xi _\mathbbm {u}}^\infty \frac{1}{(u+w^2)R_{\mathbbm {e},\mathbbm {u}}(u)}{\mathrm d} u \end{aligned}$$
(48a)
$$\begin{aligned} R_{\mathbbm {e},\mathbbm {u}}(u)&\!\!=\!\!\sqrt{\!(\!u+a_{\mathbbm {e},\!\mathbbm {u}}^2)\!(u\!+b_{\mathbbm {e},\mathbbm {u}}^2)\!(\!u\!+c_{\mathbbm {e},\!\mathbbm {u}}^2)} \end{aligned}$$
(48b)

where \(w=a_{\mathbbm {e},\mathbbm {u}},b_{\mathbbm {e},\mathbbm {u}},c_{\mathbbm {e},\mathbbm {u}}\) for \(i={x_\mathbbm {u}},{y_\mathbbm {u}},{z_\mathbbm {u}}\) and \(|\mathbbm {e}|\) is the volume of \(\mathbbm {e}\). The diagonal entries \(N_{\mathbbm {c},\mathbbm {s}}^{({x_\mathbbm {s}})},N_{\mathbbm {c},\mathbbm {s}}^{({y_\mathbbm {s}})},N_{\mathbbm {c},\mathbbm {s}}^{({z_\mathbbm {s}})}\) of \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}({\xi _\mathbbm {s}})\) are

$$\begin{aligned} N_{\mathbbm {c},\mathbbm {s}}^{(i)}({\xi _\mathbbm {s}})&=\frac{3|\mathbbm {c}|}{8\pi } \int \limits _{\xi _\mathbbm {s}}^\infty \frac{1}{(u+w^2)R_{\mathbbm {c},\mathbbm {s}}(u)}{\mathrm d} u \end{aligned}$$
(49a)
$$\begin{aligned} R_{\mathbbm {c},\mathbbm {s}}(u)&\!\!=\!\!\sqrt{\!(u+a_{\mathbbm {c},\!\mathbbm {s}}^2)\!(u+b_{\mathbbm {c},\mathbbm {s}}^2)\!(u+c_{\mathbbm {c},\!\mathbbm {s}}^2)} \end{aligned}$$
(49b)

where \(|\mathbbm {c}|\) is the volume of \(\mathbbm {c}\) and \(w=a_{\mathbbm {c},\mathbbm {s}},b_{\mathbbm {c},\mathbbm {s}},c_{\mathbbm {c},\mathbbm {s}}\) for \(i={x_\mathbbm {s}}_\mathbbm {u},{y_\mathbbm {s}}_\mathbbm {u},{z_\mathbbm {s}}_\mathbbm {u}\). The diagonal entries \(N_{\mathbbm {e},\mathbbm {s}}^{({x_\mathbbm {s}})},N_{\mathbbm {e},\mathbbm {s}}^{({y_\mathbbm {s}})},N_{\mathbbm {e},\mathbbm {s}}^{({z_\mathbbm {s}})}\) of \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}({\xi _\mathbbm {s}})\) are

$$\begin{aligned} N_{\mathbbm {e},\mathbbm {s}}^{(i)}({\xi _\mathbbm {s}})&=\frac{3|\mathbbm {e}_{\mathbbm {s}}|}{8\pi } \int \limits _{\xi _\mathbbm {s}}^\infty \frac{1}{(u+w^2)R_{\mathbbm {e},\mathbbm {s}}(u)}{\mathrm d} u \end{aligned}$$
(50a)
$$\begin{aligned} R_{\mathbbm {e},\mathbbm {s}}(u)&\!\!=\!\!\sqrt{\!(u\!+a_{\mathbbm {e},\!\mathbbm {s}}^2)\!(u+b_{\mathbbm {e},\!\mathbbm {s}}^2)\!(u+c_{\mathbbm {e},\!\mathbbm {s}}^2)} \end{aligned}$$
(50b)

where \(w=a_{\mathbbm {e},\mathbbm {s}},b_{\mathbbm {e},\mathbbm {s}},c_{\mathbbm {e},\mathbbm {s}}\) for \(i={x_\mathbbm {s}}_\mathbbm {u},{y_\mathbbm {s}}_\mathbbm {u},{z_\mathbbm {s}}_\mathbbm {u}\) and \(|\mathbbm {e}_\mathbbm {s}|=|\mathbbm {e}|\) is the volume of \(\mathbbm {e}\). Then

$$\begin{aligned} {{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({\xi _\mathbbm {u}})&=({{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {u}}^{-1})^{\!\mathrm {T}}{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}}({\xi _\mathbbm {u}}){{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {u}}^{-1} \end{aligned}$$
(51a)
$$\begin{aligned} {{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}({\xi _\mathbbm {s}})&=({{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{-1})^{\!\mathrm {T}}{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}({\xi _\mathbbm {s}}){{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{-1} \end{aligned}$$
(51b)
$$\begin{aligned} {{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}}({\xi _\mathbbm {s}})&=({{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{-1})^{\!\mathrm {T}}{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}({\xi _\mathbbm {s}}){{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{-1} \end{aligned}$$
(51c)

are the depolarization tensors in (xyz)-coordinates. With the abbreviations \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}}={{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}}(0)\), \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}={{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}(0)\), \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}={{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}(0)\), \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}={{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}(0)\) etc., one has

$$\begin{aligned} {\varvec{{\mathsf {G}}}}_\mathbbm {c}&=\left\{ \phi _\mathbbm {c}{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}{\varvec{\mathsf {\Delta }}} +\left[ {\varvec{{\mathsf {1}}}}+{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}-{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\right] \right. \nonumber \\&\qquad \left. \left[ {\varvec{{\mathsf {1}}}}+({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}-\phi _\mathbbm {c}{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}}) {\varvec{\mathsf {\Delta }}}\right] \right\} ^{-1} \end{aligned}$$
(52a)
$$\begin{aligned} {\varvec{{\mathsf {G}}}}_\mathbbm {s}({\xi _\mathbbm {s}})&= -\left[ {\varvec{{\mathsf {1}}}}+{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}{\varvec{\mathsf {\Delta }}}\right] {\varvec{{\mathsf {G}}}}_\mathbbm {c} +{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}({\xi _\mathbbm {s}}){\varvec{\mathsf {\Delta }}}{\varvec{{\mathsf {G}}}}_\mathbbm {c} \end{aligned}$$
(52b)
$$\begin{aligned} {\varvec{{\mathsf {G}}}}_\mathbbm {u}({\xi _\mathbbm {u}})&= -{\varvec{{\mathsf {1}}}} + {{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({\xi _\mathbbm {u}}) \left\{ \phi _\mathbbm {c}{\varvec{\mathsf {\Delta }}}+({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}-{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\right. \nonumber \\&\qquad \left. \left[ {\varvec{{\mathsf {1}}}}+({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}-\phi _\mathbbm {c}{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}}) {\varvec{\mathsf {\Delta }}}\right] \right\} {\varvec{{\mathsf {G}}}}_\mathbbm {c} \end{aligned}$$
(52c)

where \(\phi _\mathbbm {c}=|\mathbbm {c}|/|\mathbbm {e}|\) is the local core volume fraction and

$$\begin{aligned} {\varvec{\mathsf {\Delta }}}={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}}-{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}} \end{aligned}$$
(52d)

is the dielectric difference between the core and the shell material.

6 Self-consistent approximation

The potentials \(U_\mathbbm {c}({{\varvec{\mathrm {r}}}};{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\), \(U_\mathbbm {s}({{\varvec{\mathrm {r}}}};{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\), and \(U_\mathbbm {u}({{\varvec{\mathrm {r}}}};{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\) in Eq. (45) depend on the data \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}\). The dielectric permittivity \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) of the coated ellipsoidal inclusion is determined as the solution of the equation

$$\begin{aligned} U_\mathbbm {u}({{\varvec{\mathrm {r}}}};{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})= U_\mathbbm {u}({{\varvec{\mathrm {r}}}};{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}) \end{aligned}$$
(53)

for all \({{\varvec{\mathrm {r}}}}\in \mathbbm {u}\) by demanding that the potential outside of the inclusion is identical to that of the full solution. Setting \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}}={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) in Eq. (52) gives for the potential

$$\begin{aligned} {\varvec{{\mathsf {G}}}}_\mathbbm {e}&= \left[ {\varvec{{\mathsf {1}}}}+{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}-{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\right] ^{-1} \end{aligned}$$
(54a)
$$\begin{aligned} {\varvec{{\mathsf {G}}}}_\mathbbm {u}\!&=\! -{\varvec{{\mathsf {1}}}}\!+\!{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({\xi _\mathbbm {u}})({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\!-\!{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}){\varvec{{\mathsf {G}}}}_\mathbbm {e} \end{aligned}$$
(54b)

and

$$\begin{aligned} {\varvec{\mathrm {\mathrm {E}}}}_\mathbbm {e}&= \left[ {\varvec{{\mathsf {1}}}}+{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}-{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\right] ^{-1}{\varvec{\mathrm {\mathrm {E}}}}_{\mathrm {ext}} \end{aligned}$$
(55a)
$$\begin{aligned} {\varvec{\mathrm {\mathrm {D}}}}_\mathbbm {e}&= {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}{\varvec{\mathrm {\mathrm {E}}}}_\mathbbm {e} \end{aligned}$$
(55b)

for the fields. Equating Eq. (54b) with Eq. (52c) yields

$$\begin{aligned} ({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}-{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})&\left[ {\varvec{{\mathsf {1}}}}+{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\!-\!{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\right] ^{-1} =\nonumber \\&[\phi _\mathbbm {c}{\varvec{\mathsf {\Delta }}}+({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}\!-\!{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}){\varvec{{\mathsf {B}}}}]{\varvec{{\mathsf {G}}}}_\mathbbm {c} \end{aligned}$$
(56)

where

$$\begin{aligned} {\varvec{{\mathsf {B}}}} = {\varvec{{\mathsf {1}}}}\!+\!({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}\!\!-\!\phi _\mathbbm {c}{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}}){\varvec{\mathsf {\Delta }}} . \end{aligned}$$
(57)

and

$$\begin{aligned} {\varvec{{\mathsf {G}}}}_\mathbbm {c}\!=\!\left\{ \phi _\mathbbm {c}{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}{\varvec{\mathsf {\Delta }}} \!+\!\left[ {\varvec{{\mathsf {1}}}}+{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}-{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\right] {\varvec{{\mathsf {B}}}}\right\} ^{-1} . \end{aligned}$$
(58)

Solving Eq. (56) for \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) gives

$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}} {=} {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}} {+} \left[ \left( [\phi _\mathbbm {c}{\varvec{\mathsf {\Delta }}}{+}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}\!-\!{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}){\varvec{{\mathsf {B}}}}]{\varvec{{\mathsf {G}}}}_\mathbbm {c}\right) ^{-1} {-}{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}\right] ^{-1} \end{aligned}$$
(59)

which simplifies to become

$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}};\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s})= {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}+\phi _\mathbbm {c}{\varvec{\mathsf {\Delta }}}{\varvec{{\mathsf {B}}}}^{-1} \end{aligned}$$
(60)

for the dielectric permittivity of the ellipsoidal inclusion in terms of its geometric and material parameters. The geometric parameters are restricted to \(\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s}\) because \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}\) follows from \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}}\) and the latter depends only on \(d_\mathbbm {s},e_\mathbbm {s}\).

Following Ref. [58] and Eq. (13), the effective dielectric permittivity \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) is defined by averaging Eq. (55) over the disorder to get

$$\begin{aligned}&\left\langle {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}{\varvec{\mathrm {\mathrm {E}}}}_\mathbbm {e}\right\rangle = {\overline{{\varvec{\mathsf {\upepsilon }}}}}\left\langle {\varvec{\mathrm {\mathrm {E}}}}_\mathbbm {e}\right\rangle \end{aligned}$$
(61a)
$$\begin{aligned} \Longrightarrow \quad&\left\langle ({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}-{\overline{{\varvec{\mathsf {\upepsilon }}}}}) \!\left[ {\varvec{{\mathsf {1}}}}\!+\!{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\!\mathbbm {u}}}\!(\!{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\!-\!{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}\!)\!\right] ^{-1}\right\rangle \!= 0 \end{aligned}$$
(61b)

where \(\left\langle \dots \right\rangle \) is the disorder average. The implicit dependencies of \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) are clear from Eq. (60). The implicit dependencies of \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}\) emerge from noting that the quadratic form matrix \({\varvec{{\mathsf {Q}}}}_\mathbbm {e}\) of the ellipsoidal inclusion from Eq. (35a) can be expressed equivalently as

$$\begin{aligned} {\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}&= {{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {u}} \begin{pmatrix} a_{\mathbbm {e},\mathbbm {u}}^{2} &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad b_{\mathbbm {e},\mathbbm {u}}^{2} &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad c_{\mathbbm {e},\mathbbm {u}}^{2} \end{pmatrix} {{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {u}}^{\!\mathrm {T}} \end{aligned}$$
(62a)
$$\begin{aligned}&= {{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}} \begin{pmatrix} a_{\mathbbm {e},\mathbbm {s}}^{2} &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad b_{\mathbbm {e},\mathbbm {s}}^{2} &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad c_{\mathbbm {e},\mathbbm {s}}^{2} \end{pmatrix} {{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{\!\mathrm {T}} \end{aligned}$$
(62b)

in \({{{\varvec{\mathrm {r}}}}_\mathbbm {u}}\)-coordinates or \({{{\varvec{\mathrm {r}}}}_\mathbbm {s}}\)-coordinates. Noting that Eq. (48) can be written as

$$\begin{aligned} {{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}} =\frac{\sqrt{\det ({\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1})}}{2}\int \limits _0^\infty \frac{[{\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}+u{\varvec{{\mathsf {1}}}}]^{-1}}{\sqrt{\det ({\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}+u{\varvec{{\mathsf {1}}}})}}{\mathrm d} u \end{aligned}$$
(63)

allows to rewrite \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}\) as

$$\begin{aligned} {{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}&= \frac{1}{2\sqrt{\det ({\varvec{{\mathsf {Q}}}}_\mathbbm {e})}}\int \limits _0^\infty \frac{[{\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}+s{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}]^{-1}}{\sqrt{\det ({\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}+s{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})}}{\mathrm d} s \end{aligned}$$
(64)

and this shows that \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}={{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}})\) depends on \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}\) and \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}\).

Self-consistency determines \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) and \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) such that the background equals the effective medium. Setting

$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}&={\overline{{\varvec{\mathsf {\upepsilon }}}}} \end{aligned}$$
(65a)
$$\begin{aligned} {{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}&={{\overline{{{\varvec{{\mathsf {L}}}}}}}} \end{aligned}$$
(65b)

in Eq. (61b) yields

$$\begin{aligned} \left\langle ({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}-{\overline{{\varvec{\mathsf {\upepsilon }}}}}) \left\{ {\varvec{{\mathsf {1}}}}+ {{\overline{{{\varvec{{\mathsf {L}}}}}}}}\left[ {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}-{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] \right\} ^{-1}\right\rangle = 0 , \end{aligned}$$
(66)

the self-consistent effective medium approximation for \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\). The effective depolarisation tensor \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) defined as

$$\begin{aligned} {{\overline{{{\varvec{{\mathsf {L}}}}}}}}&= \frac{1}{2\sqrt{\det ({\varvec{{\mathsf {Q}}}}_\mathbbm {e})}}\int \limits _0^\infty \frac{[{\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}+s{\overline{{\varvec{\mathsf {\upepsilon }}}}}]^{-1}}{\sqrt{\det ({\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}+s{\overline{{\varvec{\mathsf {\upepsilon }}}}})}}{\mathrm d} s \end{aligned}$$
(67)

with \({\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}\) given in Eq. (62b) depends on \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\). In addition to \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) also the shape of the ellipsoidal inclusion \(\mathbbm {e}\) has to be determined self-consistently. This differs from conventional Bruggemann theory.

7 Local porosity theory

The self-consistent equation of local porosity theory is obtained by expressing \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) and \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) in terms of local geometric quantities and replacing the disorder average with an average over the local geometry distribution \(\mu \). This gives

$$\begin{aligned} \left\langle \left\{ ({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}-{\overline{{\varvec{\mathsf {\upepsilon }}}}})^{-1} +{{\overline{{{\varvec{{\mathsf {L}}}}}}}}\;\right\} ^{-1}\right\rangle _\mu = 0 \end{aligned}$$
(68a)

with \(\langle \dots \rangle _\mu \) given in Eq. (33). The dependencies

$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}&={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}};\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s}) \end{aligned}$$
(68b)
$$\begin{aligned} {{\overline{{{\varvec{{\mathsf {L}}}}}}}}&={{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}};d_\mathbbm {s},e_\mathbbm {s}) \end{aligned}$$
(68c)

have to be related to the materials filling \(\mathbbm {p}\) and \(\mathbbm {m}\). The cases of \(\Lambda =0\) (nonpercolating local geometry) and \({\Lambda =1}\) (percolating local geometry) need to be distinguished, because for percolating local geometries the shell \(\mathbbm {s}\) should contain \(\mathbbm {p}\)-material, while it should contain \(\mathbbm {m}\)-material in the nonpercolating case. The core and shell properties are specified as

$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}}&=(1-\Lambda ){{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}}+\Lambda {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}} \end{aligned}$$
(69a)
$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}&=(1-\Lambda ){{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}}+\Lambda {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}} \end{aligned}$$
(69b)
$$\begin{aligned} \phi _\mathbbm {c}&=(1-\Lambda )\phi +\Lambda (1-\phi ) \end{aligned}$$
(69c)
$$\begin{aligned} d_\mathbbm {s}&=(1-\Lambda )d_\mathbbm {m}+\Lambda d_\mathbbm {p} \end{aligned}$$
(69d)
$$\begin{aligned} e_\mathbbm {s}&=(1-\Lambda )e_\mathbbm {m}+\Lambda e_\mathbbm {p} . \end{aligned}$$
(69e)

Inserting Eq. (69) into Eq. (68) and using (33) gives

$$\begin{aligned}&\!\!\!\!\!\!0=\int \limits _0^1\int \limits _{0}^{1}\!\!\int \limits _0^{1} \uplambda (\phi ) \left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}};1-\phi ,d_\mathbbm {p},e_\mathbbm {p})-{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] ^{-1}\right. \nonumber \\&\!\!\!+\left. {{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}};d_\mathbbm {p},e_\mathbbm {p})\right\} ^{-1} \mu _1(\phi ,d_\mathbbm {p},e_\mathbbm {p}){\mathrm d}e_\mathbbm {p}{\mathrm d}d_\mathbbm {p}{\mathrm d}\phi \nonumber \\&\!\!\!+\int \limits _0^1\int \limits _{0}^{1}\!\!\int \limits _0^{1} (1-\uplambda (\phi )) \left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}};\phi ,d_\mathbbm {m},e_\mathbbm {m})-{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] ^{-1}\right. \nonumber \\&\!\!\!+\left. {{\overline{{{\varvec{{\mathsf {L}}}}}}}} ( {\overline{{\varvec{\mathsf {\upepsilon }}}}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}};d_\mathbbm {m},e_\mathbbm {m} )\right\} ^{-1} \mu _0(\phi ,d_\mathbbm {m},e_\mathbbm {m})\;{\mathrm d}e_\mathbbm {m}{\mathrm d}d_\mathbbm {m}{\mathrm d}\phi \end{aligned}$$
(70)

as the mixing law of local porosity theory. The unknown quantity in Eq. (70) is \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\). The functional \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) is given in Eq. (60) with Eqs. (57) and (52d), while \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) is from eq.(67).

8 Analytical results

8.1 Small measurement cells

In the limit \(|\mathbbm {K}|\rightarrow 0\) of pointlike measurement cells the mixing law of local porosity theory becomes

$$\begin{aligned} 0&= {\overline{\phi }}\left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}};\!0,\!1,\!1) \!-\!{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] ^{-1} \!\!+{{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}};\! 1,\! 1) \right\} ^{-1} \nonumber \\&+(\!1\!-\!{\overline{\phi }} ) \left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}}\!,{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}}\!;\!1,\!1,\!1\!)\! -\!{\overline{{\varvec{\mathsf {\upepsilon }}}}}\!\right] ^{-1} \!\!+\!\!{{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}}\!,\!{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}}\!;1\!,1\!) \right\} ^{-1}\nonumber \\{} \end{aligned}$$
(71)

after inserting Eq. (30) into Eq. (70). For isotropic materials \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}}=\upepsilon _\mathbbm {p}{\varvec{{\mathsf {1}}}}}, {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}}=\upepsilon _\mathbbm {m}{\varvec{{\mathsf {1}}}}\) the classic Bruggemann effective medium approximation (EMA) equation

$$\begin{aligned} 0={\overline{\phi }}\frac{\upepsilon _\mathbbm {p}-{\overline{\upepsilon }}}{\upepsilon _\mathbbm {p}+2{\overline{\upepsilon }}}+ (1-{\overline{\phi }})\frac{\upepsilon _\mathbbm {m}-{\overline{\upepsilon }}}{\upepsilon _\mathbbm {m}+2{\overline{\upepsilon }}} \end{aligned}$$
(72)

ensues, which should be compared to Eq. (1). The known anisotropic effective medium formula (2) is recovered setting \(\mu (1,\phi ,d_\mathbbm {p},e_\mathbbm {p})\,{=}\, {\overline{\phi }}\delta \Big [(\phi ,d_\mathbbm {p}, e_\mathbbm {p})\!-\!(1,\!1,\!c/b)\Big ]\) and \(\mu (0,\phi ,d_\mathbbm {m},e_\mathbbm {m})=(1\!-\!{\overline{\phi }}) \delta [(\phi ,d_\mathbbm {m},e_\mathbbm {m}) -\!(0,\!1,\!c/b)]\) in Eq. (30), where bc are the parameters in Eq. (2). The limit of small measurement cells is the limit in which the basic assumption of statistical independence underlying effective medium approximations is violated. It is also worth to remark that, in a limiting sense, Eq. (2) requires measurement cells to remain anisotropic when they shrink to a point.

8.2 Large measurement cells

For large measurement cells \(|\mathbbm {K}|\rightarrow \infty \) Eq. (32) gives

$$\begin{aligned} 0=&\uplambda ({\overline{\phi }}) \left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}};1-{\overline{\phi }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})-{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] ^{-1}\right. \nonumber \\&\qquad \qquad \qquad +\left. {{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}};{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\right\} ^{-1} \nonumber \\&+ (1-\uplambda ({\overline{\phi }})) \left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}};{\overline{\phi }},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})-{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] ^{-1}\right. \nonumber \\&\qquad \qquad \qquad +\left. {{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}};{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\right\} ^{-1} . \end{aligned}$$
(73)

For geometrical isotropy \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}={{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}={\varvec{{\mathsf {1}}}}/3\) and material isotropy \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}}=\upepsilon _\mathbbm {p}{\varvec{{\mathsf {1}}}}}\), \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}}=\upepsilon _\mathbbm {m}{\varvec{{\mathsf {1}}}}\) this becomes

$$\begin{aligned} 0=\uplambda ({\overline{\phi }})\frac{{\upepsilon _{C}}-{\overline{\upepsilon }}}{{\upepsilon _{C}}+2{\overline{\upepsilon }}}+ (1-\uplambda ({\overline{\phi }}))\frac{{\upepsilon _{B}}-{\overline{\upepsilon }}}{{\upepsilon _{B}}+2{\overline{\upepsilon }}} \end{aligned}$$
(74)

the generalized effective medium approximation from Ref. [36], Sec.V.A, where

$$\begin{aligned} {\upepsilon _{C}}&={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}(\upepsilon _\mathbbm {m},\upepsilon _\mathbbm {p};1-{\overline{\phi }},1,1) \end{aligned}$$
(75a)
$$\begin{aligned} {\upepsilon _{B}}&={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}(\upepsilon _\mathbbm {p},\upepsilon _\mathbbm {m};{\overline{\phi }},1,1) \end{aligned}$$
(75b)

are the local dielectric functions.

8.3 Geometrical isotropy

Geometrical isotropy is the special case where the local geometry is spherical instead of ellipsoidal. Computing the integrals gives

$$\begin{aligned} {{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}={{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}=\frac{1}{3}\;{\varvec{{\mathsf {1}}}} \end{aligned}$$
(76)

and hence

$$\begin{aligned} {{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}={{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}}=\left( 3{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}\right) ^{-1} . \end{aligned}$$
(77)

Inserting this into Eq. (57) and (60) gives

$$\begin{aligned}&{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}};\phi _\mathbbm {c}) = {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}};\phi _\mathbbm {c},1,1) = \nonumber \\&\upepsilon _\mathbbm {s}+ \phi _\mathbbm {c}\left\{ ({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}}-{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}})^{-1}+ (1-\phi _\mathbbm {c})\left( 3{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}\right) ^{-1} \right\} ^{-1} \end{aligned}$$
(78)

for coated spheres. In this case Eq. (70) becomes

$$\begin{aligned}&0=\int \limits _0^1\!\left[ \uplambda (\phi ) \left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}};1\!-\!\phi )\!-\!{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] ^{-1} \!\!+{{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}})\right\} ^{-1}\right. \nonumber \\&\quad \left. + (1\!-\!\uplambda (\phi )) \left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}};\phi )\!-\!{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] ^{-1} \!\!+{{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}})\right\} ^{-1} \right] \nonumber \\&\quad \mu (\phi )\,{\mathrm d}\phi \end{aligned}$$
(79)

where the shorthand notation \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},\,\cdot \,)={{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},\;\cdot \;;1,1)\), was used and \(\mu (\phi )=\mu _1(\phi )=\mu _0(\phi )\) is the local porosity distribution. This generalizes local porosity theory to anisotropic materials.

For \(|\mathbbm {K}|\rightarrow 0\) (pointlike measurement cells) Eqs. (30), (57), (60) and (77) for spheres give again Eq. (71) as the generalized effective medium equation for \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) in the case of anisotropic media.

9 Material isotropy

A situation of great practical and applied interest is the case of material isotropy but geometrical anisotropy. In this case the effective material parameters the effective dielectric properties become anisotropic only due to geometric properties of the mixture, while the materials themselves are isotropic. General results for macroscopically random systems will be investigated in the first Sect. 9.1 of this section. The non-random case will be discussed in more detail in the remaining subsections.

9.1 Macroscopically random systems

Isotropic materials occupying \(\mathbbm {p}\) and \(\mathbbm {m}\) have dielectric tensors proportional to \({\varvec{{\mathsf {1}}}}\)

$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}}&=\upepsilon _\mathbbm {p}{\varvec{{\mathsf {1}}}},\qquad \upepsilon _\mathbbm {p}\in \mathbbm {C}, \end{aligned}$$
(80a)
$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}}&=\upepsilon _\mathbbm {m}{\varvec{{\mathsf {1}}}},\qquad \upepsilon _\mathbbm {m}\in \mathbbm {C}, \end{aligned}$$
(80b)

where the scalar prefactors

$$\begin{aligned} \upepsilon _\mathbbm {p}(\omega )&=\upepsilon _\mathbbm {p}^\prime (\omega )+\mathrm {i}\frac{\upsigma _\mathbbm {p}^\prime (\omega )}{\omega } \end{aligned}$$
(80c)
$$\begin{aligned} \upepsilon _\mathbbm {m}(\omega )&=\upepsilon _\mathbbm {m}^\prime (\omega )+\mathrm {i}\frac{\upsigma _\mathbbm {m}^\prime (\omega )}{\omega } \end{aligned}$$
(80d)

are in general frequency dependent. Then \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}}=\upepsilon _\mathbbm {c}{\varvec{{\mathsf {1}}}}\), \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}=\upepsilon _\mathbbm {s}{\varvec{{\mathsf {1}}}}\) so that

$$\begin{aligned} {{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}} = \frac{1}{\upepsilon _\mathbbm {s}}{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}} \end{aligned}$$
(81a)
$$\begin{aligned} {{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}} = \frac{1}{\upepsilon _\mathbbm {c}}{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}} \end{aligned}$$
(81b)

where \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}},{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}\) and hence \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}},{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}\) are all diagonal. It follows now from Eq. (62b) that \({\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}\) is diagonal and hence from Eq. (62a) that \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}\) is diagonal, and thus that \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}\) is diagonal. As a result

$$\begin{aligned} {\overline{{\varvec{\mathsf {\upepsilon }}}}}=\text {diag}({\overline{\upepsilon }}_x,{\overline{\upepsilon }}_y,{\overline{\upepsilon }}_z) \end{aligned}$$
(82)

and \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) are diagonal. Equation (60) becomes

$$\begin{aligned}&{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}= \upepsilon _\mathbbm {s}\left[ {\varvec{{\mathsf {1}}}}{+}\phi _\mathbbm {c} \left\{ \frac{\upepsilon _\mathbbm {s}}{\upepsilon _\mathbbm {c}{-}\upepsilon _\mathbbm {s}}{\varvec{{\mathsf {1}}}} {+}{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}{-}\phi _\mathbbm {c}{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}} \right\} ^{-1}\right] \end{aligned}$$
(83)

where the diagonal matrices \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}},{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}\), given in Eqs. (49) and (50), can be written as

$$\begin{aligned} N_{\mathbbm {e},\mathbbm {s}}^{({x_\mathbbm {s}})}= & {} N_a(d_\mathbbm {s},e_\mathbbm {s}), \end{aligned}$$
(84a)
$$\begin{aligned} N_{\mathbbm {e},\mathbbm {s}}^{({y_\mathbbm {s}})}= & {} N_b(d_\mathbbm {s},e_\mathbbm {s}), \end{aligned}$$
(84b)
$$\begin{aligned} N_{\mathbbm {e},\mathbbm {s}}^{({z_\mathbbm {s}})}= & {} N_c(d_\mathbbm {s},e_\mathbbm {s}), \end{aligned}$$
(84c)
$$\begin{aligned} N_{\mathbbm {c},\mathbbm {s}}^{({x_\mathbbm {s}})}= & {} N_a\left( D_\mathbbm {c}(\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s}),E_\mathbbm {c}(\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s})\right) , \end{aligned}$$
(84d)
$$\begin{aligned} N_{\mathbbm {c},\mathbbm {s}}^{({y_\mathbbm {s}})}= & {} N_b\left( D_\mathbbm {c}(\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s}),E_\mathbbm {c}(\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s})\right) , \end{aligned}$$
(84e)
$$\begin{aligned} N_{\mathbbm {c},\mathbbm {s}}^{({z_\mathbbm {s}})}= & {} N_c\left( D_\mathbbm {c}(\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s}),E_\mathbbm {c}(\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s})\right) \end{aligned}$$
(84f)

with

$$\begin{aligned} N_a(d,e)&\!=\! \frac{e}{2d}\!\! \int \limits _0^\infty \!\! \frac{{\mathrm d}u}{\left[ u\!+\!\frac{1}{d^2}\right] ^{\!\frac{3}{2}}\!\left[ u\!+\!1\right] ^{\!\frac{1}{2}} \left[ u\!+\!e^2\right] ^{\!\frac{1}{2}} } \end{aligned}$$
(85a)
$$\begin{aligned} N_b(d,e)&\!=\! \frac{e}{2d}\!\! \int \limits _0^\infty \!\! \frac{{\mathrm d}u}{\left[ u\!+\!\frac{1}{d^2}\right] ^{\!\frac{1}{2}}\left[ u\!+\!1\right] ^{\!\frac{3}{2}} \left[ u\!+\!e^2\right] ^{\!\frac{1}{2}} } \end{aligned}$$
(85b)
$$\begin{aligned} N_c(d,e)&\!=\! \frac{e}{2d}\!\! \int \limits _0^\infty \!\! \frac{{\mathrm d}u}{\left[ u\!+\!\frac{1}{d^2}\right] ^{\!\frac{1}{2}}\left[ u\!+\!1\right] ^{\!\frac{1}{2}} \left[ u\!+\!e^2\right] ^{\!\frac{3}{2}} } \end{aligned}$$
(85c)

and \(D_\mathbbm {c},E_\mathbbm {c}\) from Eq. (42). Inserting Eq. (82) into Eq. (67), using

$$\begin{aligned} {\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1} = b_{\mathbbm {e},\mathbbm {s}}^2{{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}} \begin{pmatrix} d_\mathbbm {s}^{-2} &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 \\ 0 &{} 0 &{}e_\mathbbm {s}^{2} \end{pmatrix} {{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{\!\mathrm {T}} \end{aligned}$$
(86)

and substituting \(s = (b_{\mathbbm {e},\mathbbm {s}}^2 \upepsilon _\mathbbm {s}u)/{\overline{\upepsilon }}_y\) in (67) gives

$$\begin{aligned} {{\overline{{{\varvec{{\mathsf {L}}}}}}}}={\overline{{\varvec{\mathsf {\upepsilon }}}}}^{-1}{{\overline{{{\varvec{{\mathsf {N}}}}}}}}={{\overline{{{\varvec{{\mathsf {N}}}}}}}}{\overline{{\varvec{\mathsf {\upepsilon }}}}}^{-1} \end{aligned}$$
(87)

with

$$\begin{aligned} {{\overline{{{\varvec{{\mathsf {N}}}}}}}}&{=} \text {diag}( N_a(d_{{\overline{\upepsilon }}\mathbbm {s}},e_{{\overline{\upepsilon }}\mathbbm {s}}), N_b(d_{{\overline{\upepsilon }}\mathbbm {s}},e_{{\overline{\upepsilon }}\mathbbm {s}}), N_c(d_{{\overline{\upepsilon }}\mathbbm {s}},e_{{\overline{\upepsilon }}\mathbbm {s}}))\qquad \end{aligned}$$
(88a)

where

$$\begin{aligned} d_{{\overline{\upepsilon }}\mathbbm {s}}=d_\mathbbm {s}\sqrt{\frac{{\overline{\upepsilon }}_x}{{\overline{\upepsilon }}_y}},\qquad e_{{\overline{\upepsilon }}\mathbbm {s}}=e_\mathbbm {s}\sqrt{\frac{{\overline{\upepsilon }}_y}{{\overline{\upepsilon }}_z}} \end{aligned}$$
(88b)

are self-consistently determined parameters that depend upon the unknown \({\overline{\upepsilon }}\). Equation (68) becomes diagonal

$$\begin{aligned} \left\langle \left\{ ({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}-{\overline{{\varvec{\mathsf {\upepsilon }}}}})^{-1} +{{\overline{{{\varvec{{\mathsf {N}}}}}}}}\;{\overline{{\varvec{\mathsf {\upepsilon }}}}}^{-1}\;\right\} ^{-1}\right\rangle _\mu = 0 \end{aligned}$$
(89)

with \({{\overline{{{\varvec{{\mathsf {N}}}}}}}}={{\overline{{{\varvec{{\mathsf {N}}}}}}}}(d_{{\overline{\upepsilon }}\mathbbm {s}},e_{{\overline{\upepsilon }}\mathbbm {s}})\) depending on \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\). This coupled system of three equations needs to be solved for \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) with material data from (69).

Substituting \(u=e^2s\) in Eq. (85) gives a form that is useful when \(d=0\) or \(e=0\)

$$\begin{aligned} N_a(d,e)= & {} \frac{1}{2}\int \limits _0^\infty \frac{d^2e^2\,{\mathrm d} s}{(d^2e^2s+1)^{\frac{3}{2}}(e^2s+1)^{\frac{1}{2}}(s+1)^{\frac{1}{2}}} \end{aligned}$$
(90a)
$$\begin{aligned} N_b(d,e)= & {} \frac{1}{2}\int \limits _0^\infty \frac{e^2\,{\mathrm d} s}{(d^2e^2s+1)^{\frac{1}{2}}(e^2s+1)^{\frac{3}{2}}(s+1)^{\frac{1}{2}}} \end{aligned}$$
(90b)
$$\begin{aligned} N_c(d,e)= & {} \frac{1}{2}\int \limits _0^\infty \frac{{\mathrm d} s}{(d^2e^2s+1)^{\frac{1}{2}}(e^2s+1)^{\frac{1}{2}}(s+1)^{\frac{3}{2}}} . \end{aligned}$$
(90c)

This form is convenient, because it allows one to read off the limiting cases

$$\begin{aligned} N_a(1,1)&= N_b(1,1) =N_c(1,1) = \frac{1}{3} \end{aligned}$$
(91a)
$$\begin{aligned} N_a(0,1)&= 0,\qquad N_b(0,1) =N_c(0,1) = \frac{1}{2} \end{aligned}$$
(91b)
$$\begin{aligned} N_a(1,0)&= N_b(1,0) = 0, \qquad N_c(1,0) = 1 \end{aligned}$$
(91c)
$$\begin{aligned} N_a(0,0)&= N_b(0,0) = 0, \qquad N_c(0,0) = 1 \end{aligned}$$
(91d)

and

$$\begin{aligned} N_a(d,0)= & {} N_b(d,0) = 0, \qquad N_c(d,0) = 1 \end{aligned}$$
(91e)
$$\begin{aligned} N_a(0,e)= & {} 0,~N_b(0,e)\!=\!\frac{e}{1+e},~ N_c(0,e)\!=\!\frac{1}{1+e} \end{aligned}$$
(91f)
$$\begin{aligned} N_a(1,e)= & {} N_b(1,e), \end{aligned}$$
(91g)
$$\begin{aligned} N_b(d,1)= & {} N_c(d,1), \end{aligned}$$
(91h)

for all \(0\le d,e\le 1\).

Compared to the elliptic integrals in Eqs. (48), (49), (50) the parameters \(d_{{\overline{\upepsilon }}\mathbbm {s}},e_{{\overline{\upepsilon }}\mathbbm {s}}\) in Eq. (88b) are neither from the unit interval nor real, but in general complex numbers. As a consequence the well known formulae for real numbers given in Ref. [93] do not apply. The results for \(N_a,N_b,N_c\) with complex parameters de have been computed below in Appendix A.

9.2 Macroscopically non-random systems

For macroscopically homogeneous but non-random systems Eq. (32) yields the diagonal equation

$$\begin{aligned} 0= & {} \uplambda ({\overline{\phi }})\left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_{C}}-{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] ^{-1} +{{\overline{{{\varvec{{\mathsf {N}}}}}}}}(d_{{\overline{\upepsilon }}\mathbbm {p}},e_{{\overline{\upepsilon }}\mathbbm {p}})\;{\overline{{\varvec{\mathsf {\upepsilon }}}}}^{-1} \right\} ^{-1}\\&+(1-\uplambda ({\overline{\phi }})) \left\{ \left[ {{\varvec{\mathsf {\upepsilon }}}_{B}}-{\overline{{\varvec{\mathsf {\upepsilon }}}}}\right] ^{-1} +{{\overline{{{\varvec{{\mathsf {N}}}}}}}}(d_{{\overline{\upepsilon }}\mathbbm {m}},e_{{\overline{\upepsilon }}\mathbbm {m}}){\overline{{\varvec{\mathsf {\upepsilon }}}}}^{-1} \right\} ^{-1}\nonumber \end{aligned}$$
(92)

for \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) where the notation

$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_{C}}&={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}(\upepsilon _\mathbbm {m},\upepsilon _\mathbbm {p};1-{\overline{\phi }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) \end{aligned}$$
(93a)
$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_{B}}&={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}(\upepsilon _\mathbbm {p},\upepsilon _\mathbbm {m};{\overline{\phi }},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}) \end{aligned}$$
(93b)

is used again for local functions. Written in components the diagonal Eq. (92) becomes three coupled equations (\(j=x,y,z\))

$$\begin{aligned} 0&={\overline{\uplambda }}\left[ \frac{1}{{\upepsilon _{C}}_j-{\overline{\upepsilon }}_j}+\frac{{{\overline{N}}}_{\mathbbm {p} j}}{{\overline{\upepsilon }}_j} \right] ^{-1} \nonumber \\&\quad + (1-{\overline{\uplambda }})\left[ \frac{1}{{\upepsilon _{B}}_j\!-\!{\overline{\upepsilon }}_j}+\frac{{{\overline{N}}}_{\mathbbm {m} j}}{{\overline{\upepsilon }}_j} \right] ^{-1} \end{aligned}$$
(94)

where the abbreviations

$$\begin{aligned} {{\overline{N}}}_{\mathbbm {p} j}&\!=\![{{\overline{{{\varvec{{\mathsf {N}}}}}}}}(d_{{\overline{\upepsilon }}\mathbbm {p}},e_{{\overline{\upepsilon }}\mathbbm {p}})]_j\!=\! \left[ {{\overline{{{\varvec{{\mathsf {N}}}}}}}}\left( {{\overline{d}}}_\mathbbm {p}\sqrt{\frac{{\overline{\upepsilon }}_x}{{\overline{\upepsilon }}_y}}, {{\overline{e}}}_\mathbbm {p}\sqrt{\frac{{\overline{\upepsilon }}_y}{{\overline{\upepsilon }}_z}}\right) \right] _j \end{aligned}$$
(95a)
$$\begin{aligned} {{\overline{N}}}_{\mathbbm {m} j}&\!=\![{{\overline{{{\varvec{{\mathsf {N}}}}}}}}(d_{{\overline{\upepsilon }}\mathbbm {m}},e_{{\overline{\upepsilon }}\mathbbm {m}})]_j\!=\!\! \left[ {{\overline{{{\varvec{{\mathsf {N}}}}}}}}\!\left( \!{{\overline{d}}}_\mathbbm {m}\sqrt{\frac{{\overline{\upepsilon }}_x}{{\overline{\upepsilon }}_y}}, {{\overline{e}}}_\mathbbm {m}\sqrt{\frac{{\overline{\upepsilon }}_y}{{\overline{\upepsilon }}_z}}\right) \!\right] _j \end{aligned}$$
(95b)

were introduced. The local functions are obtained from Eqs. (68b), (93), (83) as

$$\begin{aligned} {\upepsilon _{C}}_j&\!=\!\upepsilon _\mathbbm {p} \frac{\upepsilon _\mathbbm {p}\!+\! (\upepsilon _\mathbbm {m}\!-\!\upepsilon _\mathbbm {p})({N_C}_j+1-{\overline{\phi }})}{\upepsilon _\mathbbm {p}+(\upepsilon _\mathbbm {m}-\upepsilon _\mathbbm {p})\;{N_C}_j} \end{aligned}$$
(96a)
$$\begin{aligned} {{\varvec{\mathsf {\upepsilon }}}_{B}}_j&\!=\! \upepsilon _\mathbbm {m} \frac{\upepsilon _\mathbbm {m}\!+\! (\upepsilon _\mathbbm {p}\!-\!\upepsilon _\mathbbm {m})({N_B}_j+{\overline{\phi }})}{\upepsilon _\mathbbm {m}+(\upepsilon _\mathbbm {p}-\upepsilon _\mathbbm {m})\;{N_B}_j} \end{aligned}$$
(96b)

where

$$\begin{aligned}&{N_C}_x({\overline{\phi }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) \!=\! N_a\!\left( D_\mathbbm {p},E_\mathbbm {p}\right) \!-\!(1\!-\!{\overline{\phi }}) N_a({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) \end{aligned}$$
(96c)
$$\begin{aligned}&{N_C}_y({\overline{\phi }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) \!=\! N_b\!\left( D_\mathbbm {p},E_\mathbbm {p}\right) -(1\!-\!{\overline{\phi }}) N_b({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\qquad \end{aligned}$$
(96d)
$$\begin{aligned}&{N_C}_z({\overline{\phi }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) \!=\! N_c\!\left( D_\mathbbm {p},E_\mathbbm {p}\right) -(1\!-\!{\overline{\phi }}) N_c({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\quad \end{aligned}$$
(96e)
$$\begin{aligned}&{N_B}_x({\overline{\phi }},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}) \!=\! N_a\!\left( D_\mathbbm {m},E_\mathbbm {m}\right) -{\overline{\phi }}N_a({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}) \end{aligned}$$
(96f)
$$\begin{aligned}&{N_B}_y({\overline{\phi }},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}) \!=\! N_b\!\left( D_\mathbbm {m},E_\mathbbm {m}\right) -{\overline{\phi }} N_b({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}) \end{aligned}$$
(96g)
$$\begin{aligned}&{N_B}_z({\overline{\phi }},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}) \!=\! N_c\!\left( D_\mathbbm {m},E_\mathbbm {m}\right) -{\overline{\phi }}N_c({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}) \end{aligned}$$
(96h)

are the diagonal elements in \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}\), \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}\) from Eq. (84), and \(\upepsilon _\mathbbm {p}\), \(\upepsilon _\mathbbm {m}\) are from Eq. (80). Here

$$\begin{aligned} D_\mathbbm {p}&= D_\mathbbm {c}(1\!-\!{\overline{\phi }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) \end{aligned}$$
(96i)
$$\begin{aligned} E_\mathbbm {p}&= E_\mathbbm {c}(1\!-\!{\overline{\phi }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) \end{aligned}$$
(96j)
$$\begin{aligned} D_\mathbbm {m}&= D_\mathbbm {c}({\overline{\phi }},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}) \end{aligned}$$
(96k)
$$\begin{aligned} E_\mathbbm {m}&= E_\mathbbm {c}({\overline{\phi }},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}) \end{aligned}$$
(96l)

from Eq. (42). The material parameters \(\upepsilon _\mathbbm {p}^\prime ,\upepsilon _\mathbbm {m}^\prime ,\upsigma _\mathbbm {p}^\prime ,\upsigma _\mathbbm {m}^\prime \) in \(\upepsilon _\mathbbm {p},\upepsilon _\mathbbm {m}\) and the geometrical parameters \({\overline{\phi }},{\overline{\uplambda }},{{\overline{d}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {m}\) in these equations are measurable or known from experiment.

The complex effective depolarization factors for \(\mathbbm {p}\) and \(\mathbbm {m}\) depend on \({\overline{\upepsilon }}\). They couple the three equations (94) via Eq. (88b). Clearly, \({\overline{\upepsilon }}_j=0\) solves Eq. (94) for any j. For \({\overline{\upepsilon }}_j\ne 0\) one can rewrite the system of equations as

$$\begin{aligned} \frac{{\overline{\uplambda }}{\widehat{\upepsilon }}_j}{M_j-{\widehat{\upepsilon }}_j} +\frac{(1-{\overline{\uplambda }}){\widehat{\upepsilon }}_j}{P_j-{\widehat{\upepsilon }}_j} = R_j({\widehat{\upepsilon }}) \end{aligned}$$
(97a)

where the normalized effective dielectric permittivity

$$\begin{aligned} {\widehat{\upepsilon }}_j&=\frac{{\overline{\upepsilon }}_j}{\upepsilon _\mathbbm {p}} \end{aligned}$$
(97b)

is the unknown,

$$\begin{aligned} R_j({\widehat{\upepsilon }})={\overline{\uplambda }} \left[ {{\overline{N}}}_{\mathbbm {p} j}({\widehat{\upepsilon }})-{{\overline{N}}}_{\mathbbm {m} j}({\widehat{\upepsilon }})\right] -{{\overline{N}}}_{\mathbbm {p} j}({\widehat{\upepsilon }}) \end{aligned}$$
(97c)

is the right hand side, and

$$\begin{aligned} P_j&\!=\! \frac{{\upepsilon _{C}}_j}{\upepsilon _\mathbbm {p}} \!=\! \frac{1+K-K({N_C}_ j + 1-{\overline{\phi }})}{1+K-K{N_C}_{j}} , \end{aligned}$$
(97d)
$$\begin{aligned} M_j&\!=\!\frac{{\upepsilon _{B}}_j}{\upepsilon _\mathbbm {p}}\!=\! \left( \!\frac{1}{K\!+\!1}\!\right) \frac{1+K({N_B}_{j}+{\overline{\phi }})}{1+K\,{N_B}_{j}} , \end{aligned}$$
(97e)
$$\begin{aligned} K&= \frac{\upepsilon _\mathbbm {p}}{\upepsilon _\mathbbm {m}}-1 \end{aligned}$$
(97f)

with \(j=x,y,z\) are parameters. Equation (97) can also be formulated as

$$\begin{aligned} {\widehat{\upepsilon }}_j&=A_j\left( 1\pm \sqrt{1-\frac{B_j}{A_j^2}}\right) \end{aligned}$$
(98a)
$$\begin{aligned} A_j&= \frac{R_j({\widehat{\upepsilon }})(P_j+M_j)+{\overline{\uplambda }} P_j+(1-{\overline{\uplambda }})M_j}{2(1+R_j({\widehat{\upepsilon }}))} \end{aligned}$$
(98b)
$$\begin{aligned} B_j&= \frac{R_j({\widehat{\upepsilon }})P_jM_j}{1+R_j({\widehat{\upepsilon }})} . \end{aligned}$$
(98c)

This is not a solution of Eqs. (94) but a reformulation, because \(A_j,B_j\) depend through \(R_j\) on the two ratios \({\widehat{\upepsilon }}_x/{\widehat{\upepsilon }}_y\) and \({\widehat{\upepsilon }}_y/{\widehat{\upepsilon }}_z\) .

To obtain analytical information about solutions it is useful to consider limiting and special cases. The identity

$$\begin{aligned} \sum _{j=x,y,z} \frac{{\overline{\uplambda }}{\widehat{\upepsilon }}_j}{M_j-{\widehat{\upepsilon }}_j}+\frac{(1-{\overline{\uplambda }}){\widehat{\upepsilon }}_j}{P_j-{\widehat{\upepsilon }}_j} = -1 \end{aligned}$$
(99)

holds generally. It is obtained by summation of Eq. (97a) with the help of Eq. (A13).

9.3 Spherical (Isotropic) Limit

The spherically isotropic limit is the limit \({{\overline{d}}}_\mathbbm {p}\rightarrow 1\), and \({{\overline{d}}}_\mathbbm {m}\!\rightarrow 1\), and \({{\overline{e}}}_\mathbbm {p}\rightarrow 1\), and \({{\overline{e}}}_\mathbbm {m}\!\rightarrow 1\). In this limit \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}={{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}={\varvec{{\mathsf {1}}}}/3\) and \({{\overline{N}}}_{\mathbbm {p} j}\!=\!{{\overline{N}}}_{\mathbbm {m} j}\). Inserting this into Eq. (97) renders \(P_j\), \(M_j\) independent of j and \(R_j\!=\!-{{\overline{N}}}_{\mathbbm {p} j}\) follows from Eq. (95a) and (95b). Subtracting Eq. (97a) with \(j=x\) from (97a) with \(j=y\) gives

(100)

where Eq. (97c) was used. The integral on the right hand side depends on \({\widehat{\upepsilon }}_z\), while the left hand side does not. Thus the equation can only be fulfilled for \({\widehat{\upepsilon }}_x={\widehat{\upepsilon }}_y\). Repeating the same subtraction with \(j=y\) and \(j=z\) gives \({\widehat{\upepsilon }}_y={\widehat{\upepsilon }}_z\), now invoking that the integral depends on \({\widehat{\upepsilon }}_x\). In this way the equality

$$\begin{aligned} {\widehat{\upepsilon }}_x={\widehat{\upepsilon }}_y={\widehat{\upepsilon }}_z \end{aligned}$$
(101)

and the generalized effective medium equation from Ref. [36] is recovered.

9.4 Planar (flat) limits

The planar or flat limits arise when \({{\overline{e}}}_\mathbbm {p}\rightarrow 0\) or \({{\overline{e}}}_\mathbbm {m}\rightarrow 0\). Then Eq. (41) implies \(E_\mathbbm {p}\rightarrow 0\) and \(E_\mathbbm {m}\rightarrow 0\) and Eq. (40) gives \(\varphi _{\mathbbm {c}}=1\), \({{\overline{d}}}_\mathbbm {p}^2=D_\mathbbm {p}^2\) and \({{\overline{d}}}_\mathbbm {m}^2=D_\mathbbm {m}^2\). From Eq. (96) follows \({N_C}_x,{N_C}_y,{N_B}_x,{N_B}_y\rightarrow 0\), \({N_C}_z\!\rightarrow \!{\overline{\phi }}\), \({N_B}_z\!\rightarrow \! 1-{\overline{\phi }}\). With this (97d) and (97e) give

$$\begin{aligned} P_x=P_y&= \frac{1+K{\overline{\phi }}}{1+K} = M_x = M_y \end{aligned}$$
(102a)
$$\begin{aligned} P_z&=\frac{1}{1+K(1-{\overline{\phi }})} = M_z \end{aligned}$$
(102b)

and eqs. (97a) simplify to become

$$\begin{aligned} {\widehat{\upepsilon }}_j = R_j({\widehat{\upepsilon }})(P_j-{\widehat{\upepsilon }}_j) \end{aligned}$$
(103)

for \(j=x,y,z\). Note that these equations do not decouple. For \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}\) Eq. (103) is independent of \({\overline{\uplambda }}\) and hence the solution \({\widehat{\upepsilon }}\) becomes independent of \({\overline{\uplambda }}\).

9.5 Percolation limit

For sufficiently small frequencies \(\omega \ll \omega _\mathbbm {p}\), \(\omega \ll \omega _\mathbbm {m}\) the effective dielectric function \({\widehat{\upepsilon }}_j(\omega )\rightarrow {\overline{\upsigma }}^\prime (0)\) approaches the effective d.c. conductivity and the same holds for the material parameter functions \(\upepsilon _\mathbbm {p}(\omega )\rightarrow \upsigma _\mathbbm {p}^\prime (0)\) and \(\upepsilon _\mathbbm {m}(\omega )\rightarrow \upsigma _\mathbbm {m}^\prime (0)\). The complex quantity K approaches a non-negative real value

$$\begin{aligned} \kappa =\lim _{\omega \rightarrow 0}K(\omega )= \frac{\upsigma _\mathbbm {p}^\prime (0)}{\upsigma _\mathbbm {m}^\prime (0)}-1 \end{aligned}$$
(104)

which will be called conductivity contrast. The d.c. conductivities obey \(\upsigma _\mathbbm {p}^\prime (0),\upsigma _\mathbbm {m}^\prime (0)\ge 0\). Without loss of generality it can be assumed that

$$\begin{aligned} \upsigma _\mathbbm {m}^\prime (0)\le \upsigma _\mathbbm {p}^\prime (0) \quad \iff \quad \kappa \ge 0 \end{aligned}$$
(105)

holds true, i.e. that the material of higher conductivity occupies the pore space \(\mathbbm {p}\). Vanishing contrast \(\kappa =0\) means \(\upsigma _\mathbbm {p}^\prime =\upsigma _\mathbbm {m}^\prime \).

In the limit \(\kappa \rightarrow \infty \) of infinite conductivity contrast one has for \(j=x,y,z\)

$$\begin{aligned} \lim _{\kappa \rightarrow \infty }P_j&= \frac{{\overline{\phi }}-{N_C}_{j}}{1-{N_C}_{j}} =: P_j^\infty \end{aligned}$$
(106a)
$$\begin{aligned} \lim _{\kappa \rightarrow \infty }M_j&= 0 \end{aligned}$$
(106b)

in Eq. (97d) and (97e). Then Eq. (99) becomes

$$\begin{aligned} \sum _{j=x,y,z} \frac{{\widehat{\upsigma }}^\prime _j}{P^\infty _j-{\widehat{\upsigma }}^\prime _j} = \frac{3{\overline{\uplambda }} -1}{1-{\overline{\uplambda }}} \end{aligned}$$
(107)

where

$$\begin{aligned} 0\le {\widehat{\upsigma }}^\prime _j=\frac{{\overline{\upsigma }}^\prime _j(0)}{\upsigma _\mathbbm {p}^\prime (0)}\le 1 \end{aligned}$$
(108)

is the normalized d.c.-conductivity. Because \({\widehat{\upsigma }}^\prime _j\le P_j^\infty \), the left hand side is non-negative. For \({\overline{\uplambda }}\le 1/3\) the right hand side becomes negative, and hence \({\overline{\uplambda }}=1/3\) is a percolation threshold.

Because \(M_j\!=\! 0\), Eq. (98c) gives \(B_j=0\). Thus \({\widehat{\upsigma }}^\prime _j\!=\! 2A_j\) and

$$\begin{aligned} \frac{{\widehat{\upsigma }}^\prime _j}{P_j^\infty } = \frac{{\overline{\uplambda }} + R_j({\widehat{\upsigma }}^\prime )}{1+R_j({\widehat{\upsigma }}^\prime )}\, \end{aligned}$$
(109)

in the percolation limit. These are three coupled equations, that can be reduced to two coupled equations

$$\begin{aligned} u_{xy}&:=\frac{{\widehat{\upsigma }}^\prime _x}{{\widehat{\upsigma }}^\prime _y}= \frac{P_x^\infty }{P_y^\infty } \frac{({\overline{\uplambda }}+R_x)(1+R_y)}{(1+R_x)({\overline{\uplambda }}+R_y)} \end{aligned}$$
(110a)
$$\begin{aligned} u_{yz}&:=\frac{{\widehat{\upsigma }}^\prime _y}{{\widehat{\upsigma }}^\prime _z}= \frac{P_y^\infty }{P_z^\infty } \frac{({\overline{\uplambda }}+R_y)(1+R_z)}{(1+R_y)({\overline{\uplambda }}+R_z)} \end{aligned}$$
(110b)

by taking ratios. They have to be solved for the two unknown ratios \(u_{xy}\), \(u_{yz}\). The functions \(R_x,R_y,R_z\) depend only on these two ratios. For \({\overline{\uplambda }}=1\) and for \({\overline{\uplambda }}=0\) one finds

$$\begin{aligned} {\widehat{\upsigma }}^\prime _j={\left\{ \begin{array}{ll} P_j^\infty &{}\text { for }{\overline{\uplambda }}=1\\ 0 &{}\text { for } {\overline{\uplambda }}=0 \end{array}\right. } \end{aligned}$$
(111)

because the expression \(P_j^\infty {{\overline{N}}}_{\mathbbm {p} j}/({{\overline{N}}}_{\mathbbm {p} j}-1\)) resulting from Eq. (109) becomes negative.

9.6 Degenerate cases of Eq. (97)

Neither equation (97) nor Eq. (98) can be solved for \({\widehat{\upepsilon }}\) in full generality. To simplify the equations consider the special case \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{d}}}\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}\), with \(0<{{\overline{d}}}\), \({{\overline{e}}}\le 1\). Then \({{\overline{N}}}_{\mathbbm {p} j}={{\overline{N}}}_{\mathbbm {m} j}={{\overline{N}}}_j\) and \(R_j=-{{\overline{N}}}_j\) ensues.

9.6.1 Spherical LPT

The spherical case \({{\overline{d}}}=1={{\overline{e}}}\) implies \({\overline{\upsigma }}^\prime _x={\overline{\upsigma }}^\prime _y = {\overline{\upsigma }}^\prime _z\) and \(P_x^\infty =P_y^\infty =P_z^\infty \) from Section 9.3. Then Eq. (107) recovers

$$\begin{aligned}&{\widehat{\upsigma }}^\prime _x={\widehat{\upsigma }}^\prime _y = {\widehat{\upsigma }}^\prime _z= {\left\{ \begin{array}{ll} \frac{\displaystyle (3{\overline{\uplambda }}-1){\overline{\phi }}}{\displaystyle {\overline{\uplambda }}(3-{\overline{\phi }})} &{} \text { for } \frac{1}{3}\le {\overline{\uplambda }}\le 1\\ 0 &{} \text { for } 0\le {\overline{\uplambda }}\le \frac{1}{3}\qquad \end{array}\right. } \end{aligned}$$
(112)

the known results from Ref. [36], eqs.(5.1),(6.3) and Section 9.3 above.

Note, that \({\widehat{\upsigma }}^\prime \) depends on two geometric parameters, namely porosity \({\overline{\phi }}\) and connectivity \({\overline{\uplambda }}\) of the medium. Here \({\overline{\uplambda }}=\uplambda ({\overline{\phi }})\) was assumed to be constant. In general \({\overline{\uplambda }}({\overline{\phi }})\) is a nonlinear or singular function of \({\overline{\phi }}\) (see Ref. [36], Sec. VI.C).

9.6.2 Planar (flat) LPT

In the planar case with \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}\rightarrow 0\) the ellipsoids degenerate into flat ellipses. As discussed in Section 9.4 one has \({{\overline{d}}}_\mathbbm {p}^2=D_\mathbbm {p}^2\) in this case. Inserting Eq. (C7) into Eq. (96) and (106a) yields

$$\begin{aligned} P_x^\infty =P_y^\infty =&{\overline{\phi }} \end{aligned}$$
(113a)
$$\begin{aligned} P_z^\infty =&0 \end{aligned}$$
(113b)

for the local functions. Then Eq. (109) implies that \({\widehat{\upsigma }}^\prime _z =0\) for all \({\overline{\uplambda }}\). Furthermore, \({\widehat{\upsigma }}^\prime _x={\widehat{\upsigma }}^\prime _y={\overline{\phi }}\) must hold for \({\overline{\uplambda }}=1\), because of Eq. (111).

9.6.3 Oblate LPT

Here \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}=1\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}\) with \({0<{{\overline{e}}}<1}\). Evaluating the integrals in Eq. (90) for \(d=1\) and using Eqs. (B3), (B4) yields Eq. (C3). Inserting Eq. (C3a) into Eq. (97) gives \(P_x=P_y\) and \(M_x=M_y\). Subtracting Eq. (97a) with \(j=y\) from (97a) with \(j=x\) yields

$$\begin{aligned} {\overline{\upsigma }}^\prime _x = {\overline{\upsigma }}^\prime _y \end{aligned}$$
(114)

by an argument analogous to that in Section 9.3 for the spherically isotropic limit.

With this result Eq. (109) reduces to a system of two coupled equations

$$\begin{aligned} \frac{{\widehat{\upsigma }}^\prime _y}{P_y^\infty }&= \frac{{\overline{\uplambda }} - N_b(1,{{\overline{e}}} u)}{1-N_b(1,{{\overline{e}}} u)}, \end{aligned}$$
(115a)
$$\begin{aligned} \frac{{\widehat{\upsigma }}^\prime _z}{P_z^\infty }&= \frac{{\overline{\uplambda }} - N_c(1,{{\overline{e}}} u)}{1-N_c(1,{{\overline{e}}} u)} \end{aligned}$$
(115b)

where \(u=\sqrt{{\widehat{\upsigma }}^\prime _y/{\widehat{\upsigma }}^\prime _z}\) generates the coupling, and the functions \(N_b(1,e)\), \(N_c(1,e)\) are given in Eq. (C3). They obey \(2N_b+N_c=1\). From Eq. (40) one finds

$$\begin{aligned} 1-{\overline{\phi }}=\phi _\mathbbm {c}&= \frac{E_\mathbbm {p}}{{{\overline{e}}}} \left( \frac{1-{{\overline{e}}}^2}{1-E_\mathbbm {p}^2}\right) ^{\frac{3}{2}} \end{aligned}$$
(116)

and Eqs. (106) and (96) give

$$\begin{aligned}&P_y^\infty = \frac{{\overline{\phi }}-N_b(1,E_\mathbbm {p})+(1-{\overline{\phi }})N_b(1,{{\overline{e}}})}{1-N_b(1,E_\mathbbm {p},)+(1-{\overline{\phi }})N_b(1,{{\overline{e}}})}, \end{aligned}$$
(117a)
$$\begin{aligned}&P_z^\infty = \frac{{\overline{\phi }}-N_c(1,E_\mathbbm {p})+(1-{\overline{\phi }})N_c(1,{{\overline{e}}})}{1-N_c(1,E_\mathbbm {p})+(1-{\overline{\phi }})N_c(1,{{\overline{e}}})}. \end{aligned}$$
(117b)

9.6.4 Prolate LPT

In this case \({{\overline{e}}}=1\), i.e. one has \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}=1\) and \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{d}}}\) with \(0<{{\overline{d}}}<1\).

Evaluating the integrals in Eq. (90) for \(e=1\) using Eq. (B3) and Eq. (B4) yields Eq. (C1). Inserting Eq. (C1c) into Eq. (97) shows \(P_y=P_z\) and \(M_y=M_z\). Subtracting Eq. (97a) with \(j=z\) from (97a) with \(j=y\) yields

$$\begin{aligned} {\overline{\upsigma }}^\prime _y = {\overline{\upsigma }}^\prime _z \end{aligned}$$
(118)

by an argument analogous to that in Section 9.3 for the spherically isotropic limit.

With this result Eq. (109) reduces to a system of two coupled equations

$$\begin{aligned} \frac{{\widehat{\upsigma }}^\prime _x}{P_x^\infty }&= \frac{{\overline{\uplambda }} - N_a({{\overline{d}}} u,1)}{1-N_a({{\overline{d}}} u,1)}, \end{aligned}$$
(119a)
$$\begin{aligned} \frac{{\widehat{\upsigma }}^\prime _y}{P_y^\infty }&= \frac{{\overline{\uplambda }} - N_b({{\overline{d}}} u,1)}{1-N_b({{\overline{d}}} u,1)} \end{aligned}$$
(119b)

where \(u=\sqrt{{\widehat{\upsigma }}^\prime _x/{\widehat{\upsigma }}^\prime _y}\) generates the coupling, and the functions \(N_a(d,1)\), \(N_b(d,1)\) are given in Eq. (C1). They obey \(N_a+2N_b=1\). From Eq. (40) one finds

$$\begin{aligned} 1-{\overline{\phi }}=\phi _\mathbbm {c}&= \frac{D_\mathbbm {p}^2}{{{\overline{d}}}^2} \left( \frac{1-{{\overline{d}}}^2}{1-D_\mathbbm {p}^2}\right) ^{\frac{3}{2}} \end{aligned}$$
(120)

and Eqs. (106) and (96) give

$$\begin{aligned} P_x^\infty= & {} \frac{{\overline{\phi }}-N_a(D_\mathbbm {p},1)+(1-{\overline{\phi }})N_a({{\overline{d}}},1)}{1-N_a(D_\mathbbm {p},1)+(1-{\overline{\phi }})N_a({{\overline{d}}},1)}, \end{aligned}$$
(121a)
$$\begin{aligned} P_y^\infty= & {} \frac{{\overline{\phi }}-N_b(D_\mathbbm {p},1)+(1-{\overline{\phi }})N_b({{\overline{d}}},1)}{1-N_b(D_\mathbbm {p},1)+(1-{\overline{\phi }})N_b({{\overline{d}}},1)}. \end{aligned}$$
(121b)

For \({{\overline{d}}}\approx 1\) an approximate solution of Eq. (119) can be obtained for \(1/3\le {\overline{\uplambda }}\le 1\) by setting

$$\begin{aligned} u\approx \sqrt{\frac{P_x^\infty [{\overline{\uplambda }}-N_a({{\overline{d}}},1)][1-N_b({{\overline{d}}},1)]}{P_y^\infty [{\overline{\uplambda }}-N_b({{\overline{d}}},1)][1-N_a({{\overline{d}}},1)]}} \end{aligned}$$
(122)

with \(N_a,N_b\) from Eq. (C1). This decouples Eqs. (119) and yields \({\widehat{\upsigma }}^\prime \) as a function of \({\overline{\phi }},{\overline{\uplambda }},{{\overline{d}}}\).

For \({{\overline{d}}}\approx 0\) the prolate ellipsoids become strongly elongated and Eq. (C1) can be approximated by Eq. (C2) to lowest order in d. Approximating \(N_a({{\overline{d}}} u,1)\approx 0\) and \(N_b({{\overline{d}}} u,1)\approx 1/2\) in Eq. (119) by their limiting values gives

$$\begin{aligned} {\widehat{\upsigma }}^\prime _x&\approx {\overline{\uplambda }} P_x^\infty \end{aligned}$$
(123a)
$$\begin{aligned} {\widehat{\upsigma }}^\prime _y&\approx (2{\overline{\uplambda }}-1)P_y^\infty \end{aligned}$$
(123b)

for \(1/2\le {\overline{\uplambda }}\le 1\). Then \({\widehat{\upsigma }}^\prime _j\ge 0\) implies \({\widehat{\upsigma }}^\prime _y={\widehat{\upsigma }}^\prime _z=0\) for \({\overline{\uplambda }}\le 1/2\). An anisotropic percolation threshold arises at \({\overline{\uplambda }}=1/2\). Inserting \({\widehat{\upsigma }}^\prime _y={\widehat{\upsigma }}^\prime _z=0\) into Eq. (107) implies

$$\begin{aligned} {\widehat{\upsigma }}^\prime _x \approx \frac{3{\overline{\uplambda }}-1}{2{\overline{\uplambda }}}\;P_x^\infty \end{aligned}$$
(124)

for \(1/3\le {\overline{\uplambda }}\le 1/2\). In this case, i.e. for \({{\overline{d}}}\approx 0\) an approximate solution of Eq. (119) can be obtained for \(1/3\le {\overline{\uplambda }}\le 1\) by using u from Eq. (122) with \(N_a,N_b\) from Eq. (C2).

9.6.5 Circocylindrical LPT

Here \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{d}}}=0\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}=1\). This case resembles the spherical case. Inserting Eq. (91b) decouples Eq. (119) and yields

$$\begin{aligned} {\widehat{\upsigma }}^\prime _x= & {} {\left\{ \begin{array}{ll} {\overline{\uplambda }}{\overline{\phi }} &{} \text { for } \frac{1}{2}\le {\overline{\uplambda }}\le 1\\ \frac{\displaystyle (3{\overline{\uplambda }}-1){\overline{\phi }}}{\displaystyle 2{\overline{\uplambda }}} &{} \text { for } \frac{1}{3}\le {\overline{\uplambda }}\le \frac{1}{2}\\ 0 &{} \text { for } 0\le {\overline{\uplambda }}\le \frac{1}{3} \end{array}\right. } \end{aligned}$$
(125a)
$$\begin{aligned} {\widehat{\upsigma }}^\prime _y = {\widehat{\upsigma }}^\prime _z= & {} {\left\{ \begin{array}{ll} \frac{\displaystyle (2{\overline{\uplambda }}-1){\overline{\phi }}}{\displaystyle (2-{\overline{\phi }})} &{} \text { for } \frac{1}{2}\le {\overline{\uplambda }}\le 1\\ 0 &{} \text { for } 0\le {\overline{\uplambda }}\le \frac{1}{2} \end{array}\right. } \end{aligned}$$
(125b)

as a function of \({\overline{\phi }},{\overline{\uplambda }}\). Here the case \(1/3\le {\overline{\uplambda }}\le 1/2\) in Eq. (125a) follows from inserting \({\widehat{\upsigma }}^\prime _y = {\widehat{\upsigma }}^\prime _z=0\) in Eq. (107).

9.6.6 Ellipsocylindrical LPT

In this case \({{\overline{d}}}=0\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}\) and \(0<{{\overline{e}}}<1\).

Evaluating the integrals in Eq. (90) for \(d=0\) gives eqs. (91f) resp. (C5). From Eq. (40) one finds

$$\begin{aligned} 1-{\overline{\phi }}=\phi _\mathbbm {c}&=\frac{E_\mathbbm {p}}{{{\overline{e}}}}\frac{{{\overline{e}}}^2-1}{{E_\mathbbm {p}}^2-1} \end{aligned}$$
(126)

or

$$\begin{aligned} E_\mathbbm {p}&= \frac{{{\overline{e}}}^2-1}{2{{\overline{e}}}(1-{\overline{\phi }})} +\sqrt{1+\left( \frac{{{\overline{e}}}^2-1}{2{{\overline{e}}}(1-{\overline{\phi }})}\right) ^2} \end{aligned}$$
(127)

so that Eqs. (106) and (96) yield

$$\begin{aligned} P_y^\infty&=\frac{{{\overline{e}}}-E_\mathbbm {p}}{{{\overline{e}}}(1-{{\overline{e}}}E_\mathbbm {p})} \end{aligned}$$
(128a)
$$\begin{aligned} P_z^\infty&=\frac{{{\overline{e}}}({{\overline{e}}}-E_\mathbbm {p})}{1-{{\overline{e}}}E_\mathbbm {p}} \end{aligned}$$
(128b)

for the local functions, and

$$\begin{aligned} \frac{P_y^\infty }{P_z^\infty }=\frac{1}{{{\overline{e}}}^2} \end{aligned}$$
(129)

holds true for their ratio. Dividing Eq. (109) for \(j=y\) by Eq. (109) for \(j=z\) results in the simple expression

$$\begin{aligned} \sqrt{\frac{{\overline{\upsigma }}^\prime _y}{{\overline{\upsigma }}^\prime _z}}=\frac{1}{{{\overline{e}}}} \end{aligned}$$
(130)

for the ratio of the unknown effective conductivity components. Inserting this into Eq. (95a) gives the surprising result

$$\begin{aligned} {{\overline{N}}}_{\mathbbm {p} x}=0,\qquad {{\overline{N}}}_{\mathbbm {p} y}= {{\overline{N}}}_{\mathbbm {p} z}=\frac{1}{2} \end{aligned}$$
(131)

for all \({{\overline{e}}}\), i.e. independent of \({\overline{\upsigma }}^\prime _j\) or \({{\overline{e}}}\). The result decouples eqs. (109) to give with the help of Eq. (106)

$$\begin{aligned} {\widehat{\upsigma }}^\prime _x&= {\overline{\uplambda }}\,{\overline{\phi }} \end{aligned}$$
(132a)
$$\begin{aligned} {\widehat{\upsigma }}^\prime _y&= \frac{(2{\overline{\uplambda }}-1)({{\overline{e}}}-E_\mathbbm {p})}{{{\overline{e}}}(1-{{\overline{e}}}E_\mathbbm {p})} \end{aligned}$$
(132b)
$$\begin{aligned} {\widehat{\upsigma }}^\prime _z&= \frac{(2{\overline{\uplambda }}-1)\,{{\overline{e}}}\,({{\overline{e}}}-E_\mathbbm {p})}{1-{{\overline{e}}}E_\mathbbm {p}} \end{aligned}$$
(132c)

for \(1/2\le {\overline{\uplambda }}\le 1\) where \(E_\mathbbm {p}\) is from Eq. (127). For \({\overline{\uplambda }}\le 1/2\) it follows that \({\widehat{\upsigma }}^\prime _y={\widehat{\upsigma }}^\prime _z=0\) and inserting this into Eq. (107) gives

$$\begin{aligned} {\widehat{\upsigma }}^\prime _x =\frac{(3{\overline{\uplambda }}-1)\,{\overline{\phi }}}{2{\overline{\uplambda }}} \end{aligned}$$
(133)

for \(1/3\le {\overline{\uplambda }}\le 1/2\).

10 Numerical results

10.1 Percolation transitions

The Bruggemann effective medium approximation (1) exhibits a percolation transition at porosity \({\overline{\phi }}_c=1/3\) at which the dc-conductivity \({\widehat{\upsigma }}^\prime (0)\) vanishes. While qualitatively correct, the percolation threshold is generally too high for direct application to porous and disordered media. They often become percolating at much lower volume fractions. Indeed, Archie’s empirical law [36], which relates electrical conductivity to porosity, suggests that there is no lower limit for the percolation threshold.

Fig. 2
figure 2

Comparison of effective medium approximations for \({\widehat{\upsigma }}^\prime \) as functions of porosity \({\overline{\phi }}\). The solution of Eq. (1) is shown as a solid line, the solution of Eq. (2) with \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})=(1.0,0.5)\) is shown as two dash-dotted lines corresponding to \({\widehat{\upsigma }}^\prime _{x,\mathrm {Sch}}={\widehat{\upsigma }}^\prime _{y,\mathrm {Sch}}\) (upper and long dash-dotted) and \({\widehat{\upsigma }}^\prime _{z,\mathrm {Sch}}\) (lower and short dash-dotted). In both cases the percolation threshold is \({\overline{\phi }}_c=1/3\). Two dashed lines are seen for LPT corresponding to \({\widehat{\upsigma }}^\prime _x = {\widehat{\upsigma }}^\prime _y\) (long dashed) and \({\widehat{\upsigma }}^\prime _z\) (short dashed) with \(\uplambda ({\overline{\phi }})={\overline{\phi }}^{1/3}\). This implies a percolation threshold at \({\overline{\phi }}_c=1/27\approx 3.7\%\). The axis ratios are \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) = (1.0,0.5)\) and \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})=(1.0,0.01)\)

The Bruggemann result (1) has been generalized from spheres to oblate ellipsoids in Ref. [41]. It is governed by the generalized Eq. (2) which reduces to Eq. (1) for \(b=c\). The numerical solution of the anisotropic or ellipsoidal Bruggemann equation (2) is shown as dash-dotted lines in Fig. 2 as a function of \({\overline{\phi }}\). The isotropic or spherical Bruggemann result is plotted as the solid line Fig. 2. The percolation threshold is \({\overline{\phi }}_c=1/3\) in both cases.

The numerical solution of Eq. (109) for a comparable oblate case with \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) = (1.0,0.5)\) is shown as two dashed lines in Fig. 2. The long dashed line is \({\widehat{\upsigma }}^\prime _x = {\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) is short dashed. Here, \(\uplambda ({\overline{\phi }}) = {\overline{\phi }}^{1/3}\) has been assumed for the local percolation probability, and the remaining two axis ratios are \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})=(1.0,0.01)\). The percolation threshold \({\overline{\phi }}_c\) is then found as the solution of the equation \(\uplambda ({\overline{\phi }}) = 1/3\), which gives \({\overline{\phi }}_c = 1/27 \approx 3.7\%\). This is confirmed in Fig. 2. It is much more realistic than \({\overline{\phi }}_c = 1/3\) for natural media. The value of \({\overline{\phi }}_c\) depends on the local percolation probability \(\uplambda ({\overline{\phi }})\), thereby highlighting the importance of this geometrical quantity.

The solutions of Eq. (2) and Eq. (109) both display anisotropic \({\widehat{\upsigma }}^\prime \) conductivity. In both cases \({\widehat{\upsigma }}^\prime _x = {\widehat{\upsigma }}^\prime _y\), because the anisotropy is of oblate type. However, several differences exist between Eq. (2) and Eq. (109).

While \({\widehat{\upsigma }}^\prime _{x,\mathrm {Sch}}\) and \({\widehat{\upsigma }}^\prime _{z,\mathrm {Sch}}\) in Fig. 2 do not differ much, the difference between \({\widehat{\upsigma }}^\prime _{x}\) and \({\widehat{\upsigma }}^\prime _{z}\) is significant for \({\overline{\phi }} < 1\). This is in part due to the small value of \({{\overline{e}}}_\mathbbm {m}\). In addition \({\widehat{\upsigma }}^\prime _{z}\) does not approach unity for \({\overline{\phi }}\rightarrow 1\). This behaviour is caused by the local functions and stems from the confocality condition. It appears in the limit \({\overline{\phi }}\rightarrow 1\) whenever \(e_\mathbbm {p}\ne 1\) and is absent in the prolate case when \(e_\mathbbm {p}= 1\). Of course, the two additional axis ratios \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) in Eq. (109) allow for a greater variety of solutions than for Eq. (2). Recall also that Eq. (109) is only one of the many special cases from the generalized LPT in Eq. (70).

10.2 dc-Conductivity at infinite contrast

This section gives numerical solutions of the anisotropic effective medium equations (109) for the unknown normalized dc-conductivity \({\widehat{\upsigma }}^\prime (0)\). Equations (109) hold for material isotropy and geometric anisotropy in the percolation limit, i.e. at frequency \(\omega =0\) and infinite conductivity contrast \(\kappa =\infty \). The solutions for the unknown real part of the normalized effective conductivity \({\widehat{\upsigma }}^\prime \) are displayed in twelve subfigures of Fig. 3.They depend on six parameters \(({\overline{\phi }},{\overline{\uplambda }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) with values in the six dimensional unit cube \([0,1]^6\). The porosity is fixed at \({\overline{\phi }}=0.1\) in all cases.

The central subfigure in Fig. 3 visualizes the four dimensional parameter space \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\in [0,1]^4\). The anisotropy parameters \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) for \(\mathbbm {p}\) are indicated by open symbols, the parameters \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) for \(\mathbbm {m}\) by filled symbols in the central unit square. Every parameter quadruple \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\in [0,1]^4\) is represented by a pair consisting of an open and a closed symbol. The two symbols of such a pair are connected by a dashed line.

Twelve subfigures in Fig. 3 are arranged around the central subfigure. They show \({\widehat{\upsigma }}^\prime =({\widehat{\upsigma }}^\prime _x,{\widehat{\upsigma }}^\prime _y,{\widehat{\upsigma }}^\prime _z)\) as functions of \({\overline{\uplambda }}\in [0,1]\). Because \({\overline{\phi }}=0.1\) the range of values is bounded by \({\widehat{\upsigma }}^\prime _j\le 0.1\) for \(j=x,y,z\). Note also that \({\widehat{\upsigma }}^\prime _x({\overline{\uplambda }})\ge {\widehat{\upsigma }}^\prime _y({\overline{\uplambda }})\ge {\widehat{\upsigma }}^\prime _z({\overline{\uplambda }})\) holds true for all \({\overline{\uplambda }}\) in all subfigures.

The twelve subfigures are labeled from (a) through (l) as follows

  1. (a)

    spherical - - - spherical

  2. (b)

    oblate - - - prolate

  3. (c)

    prolate - - - oblate

  4. (d)

    circular - - - circular

  5. (e)

    (nearly) prolate - - - flat

  6. (f)

    flat - - - (nearly) prolate

  7. (g)

    flat - - - flat

  8. (h)

    flat - - - circocylindrical

  9. (i)

    cylindrical - - - flat

  10. (j)

    circocylindrical - - - circocylindrical

  11. (k)

    cylindrical - - - none

  12. (l)

    none - - - cylindrical

The first word in e.g. “oblate - - - prolate” characterizes the position of the pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) in the unit square, while the second gives the position of the pair \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\).

Fig. 3
figure 3

Effective medium approximations for the real part \({\widehat{\upsigma }}^\prime \) of the normalized effective conductivity from Eq. (108) at infinite conductivity contrast \(\kappa =\infty \). The parameters are \(({\overline{\phi }},{\overline{\uplambda }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\in [0,1]^6\). The porosity is \({\overline{\phi }}=0.1\) in all cases. The central subfigure shows the parameter quadruple \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\in [0,1]^4\) as \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\)- - - - -\(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) pairs where \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) is represented by open symbols and \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) by corresponding filled symbols inside the unit square with \({{\overline{d}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m}\) on the abscissa and \({{\overline{e}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {m}\) on the ordinate. To facilitate finding them, parameter pairs belonging together are connected with a dashed line as a guide to the eye. The twelve subfigures surrounding the central subfigure each show the normalized effective conductivity components \({\widehat{\upsigma }}^\prime _x\) (solid lines), \({\widehat{\upsigma }}^\prime _y\) (dashed lines), and \({\widehat{\upsigma }}^\prime _z\) (dash-dotted lines) as functions of \({\overline{\uplambda }}\). The symbol displayed in each subfigure indicates the values of the corresponding parameter pair shown in the central subfigure

10.2.1 \({{\overline{d}}}_\mathbbm {p}\approx {{\overline{d}}}_\mathbbm {m}\) and \({{\overline{e}}}_\mathbbm {p}\approx {{\overline{e}}}_\mathbbm {m}\)

The subfigures (a), (d), (g) and (j) located on the four corners of Fig. 3 illustrate cases with \({{\overline{d}}}_\mathbbm {p}\approx {{\overline{d}}}_\mathbbm {m}\) and \({{\overline{e}}}_\mathbbm {p}\approx {{\overline{e}}}_\mathbbm {m}\). In the isotropic case (a) the three components \({\widehat{\upsigma }}^\prime _x,{\widehat{\upsigma }}^\prime _y,{\widehat{\upsigma }}^\prime _z\) are nearly identical as expected from Eq. (101). In the circular - - - circular case (d) one has \({\widehat{\upsigma }}^\prime _x\approx {\widehat{\upsigma }}^\prime _y\) but \({\widehat{\upsigma }}^\prime _z \approx 0\). In both these cases the percolation transition occurs at \({\overline{\uplambda }}=1/3\). In the flat - - - flat case (g) a second percolation transition starts to appear in \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\) at \({\overline{\uplambda }}=1/2\) while \({\widehat{\upsigma }}^\prime _z\approx 0\) remains small. Note also that \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\) are nearly linear above \({\overline{\uplambda }}\approx 1/2\). In the circocylindrical - - - circocylindrical case (j) the percolation transition in \({\widehat{\upsigma }}^\prime _x\) at \({\overline{\uplambda }}=1/3\) and at \({\overline{\uplambda }}=1/2\) in \({\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) remain clearly visible. Now all components are nearly linear above \({\overline{\uplambda }}\approx 1/2\). The main difference to the flat - - - flat case (g) is that \({\widehat{\upsigma }}^\prime _z\ne 0\) and \({\widehat{\upsigma }}^\prime _x\ne {\widehat{\upsigma }}^\prime _y\) at \({\overline{\uplambda }}=1\). Other cases with \({{\overline{d}}}_\mathbbm {p}\approx {{\overline{d}}}_\mathbbm {m}\) and \({{\overline{e}}}_\mathbbm {p}\approx {{\overline{e}}}_\mathbbm {m}\) interpolate between these four limiting cases.

10.2.2 Flat pores \({{\overline{e}}}_\mathbbm {p}\approx 0\)

In the case \({{\overline{e}}}_\mathbbm {p}\approx 0\) the pores approach flat ellipses and the open symbol falls near the lower boundary of the central square. The corresponding subfigures are (d), (f), (g) and (h). As discussed in Section 9.6.2 one expects \({\widehat{\upsigma }}^\prime _z\approx 0\) for all \({\overline{\uplambda }}\), because \(P_z^\infty \rightarrow 0\) for \({{\overline{e}}}_\mathbbm {p}\approx 0\). Secondly one expects \({\widehat{\upsigma }}^\prime _x={\widehat{\upsigma }}^\prime _y=0.1\) at \({\overline{\uplambda }}=1\), because \(P_x^\infty =P_y^\infty ={\overline{\phi }}\) for \({{\overline{e}}}_\mathbbm {p}\approx 0\). Indeed, these expectations are confirmed in all cases.

10.2.3 Cylindrical pores \(d_\mathbbm {p}\approx 0\)

For \({{\overline{d}}}_\mathbbm {p}\approx 0\) the pores tend to become cylindrical and the open symbol \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) falls near the left boundary of the central square. This case is seen in subfigures (g), (i), (j) and (k). In all these cases the percolation transition in \({\widehat{\upsigma }}^\prime _x\) occurs at threshold 1/3, while \({\widehat{\upsigma }}^\prime _y\) approaches a percolation transition with threshold \({\overline{\uplambda }}=1/2\) as \({{\overline{d}}}_\mathbbm {p}\rightarrow 0\).

10.2.4 Exchanging \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) with \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\)

The parameter pairs \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) and \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) have very different impact on the result. This is illustrated in an approximate way by exchanging them.

Subfigures (b) and (c) illustrate such an exchange for parameters near the prolate and oblate boundary. Note that the curvature of \({\widehat{\upsigma }}^\prime _y\) changes.

Subfigures (e) and (f) illustrate the exchange between positions close to the flat boundary and positions somewhat removed from the upper prolate boundary. Again the curvature of \({\widehat{\upsigma }}^\prime _y\) changes, but in addition the change in \({\widehat{\upsigma }}^\prime _z\) is significant. In the case (f) the component \({\widehat{\upsigma }}^\prime _z\) is nearly zero for all \({\overline{\uplambda }}\), while it becomes strongly curved upwards in subfigure (e). This indicates a percolation transition at \({\overline{\uplambda }} = 1\) in the limit \({{\overline{e}}}_\mathbbm {m} \rightarrow 0\). Keep in mind, that the condition \({\overline{\uplambda }} = 1\) can be fulfilled for any \(0\le {\overline{\phi }}\le 1\).

The subfigures (h) and (i) illustrate the exchange between points near the flat and cylindrical boundary. Subfigure (h) shows that \({\widehat{\upsigma }}^\prime _z \approx 0\) holds for all \({\overline{\uplambda }}\). There is no percolation transition in \({\widehat{\upsigma }}^\prime _y\) at \({\overline{\uplambda }} = 1/2\) in (h). In (i) a percolation transition in \({\widehat{\upsigma }}^\prime _y\) occurs at \({\overline{\uplambda }} = 1/2\), while the strong curvature of \({\widehat{\upsigma }}^\prime _z\) indicates again a percolation transition at \({\overline{\uplambda }} = 1\) for \({{\overline{e}}}_\mathbbm {m} \rightarrow 0\).

The cylindrical cases with \({{\overline{d}}}_\mathbbm {p} \approx 0\) or \({{\overline{d}}}_\mathbbm {m} \approx 0\) are shown in subfigures (k) and (l). In (k) the percolation transition for \({\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) at \({\overline{\uplambda }} = 1/2\) is clearly visible. This transition does not appear in subfigure (l).

The eight subfigures illustrate the difference between the parameter pairs \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) and \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\). It is mainly the pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) that determines whether or not an additional percolation transition emerges at \({\overline{\uplambda }}=1/2\). It is also responsible for the smallness of \({\widehat{\upsigma }}^\prime _z\) in the limit \({{\overline{d}}}_\mathbbm {p} \rightarrow 0\) or \({{\overline{e}}}_\mathbbm {p} \rightarrow 0\).

10.2.5 Interpolation

The twelve results for \({\widehat{\upsigma }}^\prime \) shown in Fig. 3 depend in many cases continuously on the parameters \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m}, {{\overline{e}}}_\mathbbm {m})\). The continuous parameter dependence can be used to interpolate between different cases and to envisage results for other cases. This will be illustrated with two examples.

Consider first the sequence (d) \(\rightarrow \) (f) \(\rightarrow \) (h) \(\rightarrow \) (j). In the first transition the pair \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) moves from the lower boundary (flat) to the upper boundary (prolate) while \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) moves slightly away from the oblate boundary. The change in \({\widehat{\upsigma }}^\prime \) is small and mainly seen in \({\widehat{\upsigma }}^\prime _y\). In the transition (f) \(\rightarrow \) (h) both parameter pairs move towards the cylindrical boundary. The result is that the curvature in \({\widehat{\upsigma }}^\prime _y\) becomes much more pronounced and also \({\widehat{\upsigma }}^\prime _x\) becomes curved. The curvature in \({\widehat{\upsigma }}^\prime _x\) is due to the percolation transition that emerges at \({\overline{\uplambda }}=1/2\), because the open symbol has moved to the left. This is confirmed in the transition (h) \(\rightarrow \) (j) as the parameter pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) moves even closer to the cylindrical boundary. The change in \({\widehat{\upsigma }}^\prime _y\) on the other hand is caused by the closed symbol being close to the cylindrical boundary.

Next, consider the sequence (c) \(\rightarrow \) (e) \(\rightarrow \) (i) \(\rightarrow \) (l). In the first transition (c) \(\rightarrow \) (e) the pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) remains essentially fixed, while \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) moves from prolate to flat. As a result \({\widehat{\upsigma }}^\prime _z\) becomes strongly curved upwards, indicating a percolation transition at \({\overline{\uplambda }} = 1\). Shifting \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) close to the cylindrical boundary and \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) only slightly keeps the upward curvature in \({\widehat{\upsigma }}^\prime _z\) and exposes the emerging percolation transition at \({\overline{\uplambda }}=1/2\) in \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\). The transition (i) \(\rightarrow \) (l) then undoes most of these changes and suggests to close the circle with an additional transition (l) \(\rightarrow \) (c).

10.3 ac-Conductivity at infinite contrast

Figure 4 shows numerical solutions of Eq. (94) for the normalized frequency dependent conductivity tensor \({\widehat{\upsigma }}^\prime (\omega )\) at infinite conductivity contrast \(\kappa = \infty \). The line styles for \({\widehat{\upsigma }}^\prime _x,{\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) in Fig. 4 are the same as in Fig. 3. Each of the subfigures in Fig. 4 has a fixed parameter quadruple \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) and is labeled (a), (d), (g), (j), (e), and (i), to indicate the analogous pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) - - - - \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) labelled identically in the central subfigure of Fig. 3. The six subfigures in Fig. 4 are further subdivided into four subfigures. The two columns are distinguished by \({\overline{\uplambda }} = 0.35\) (left column) and \({\overline{\uplambda }} = 0.55\) (right column). The values are chosen in this way, because they are slightly above the two percolation thresholds at \({\overline{\uplambda }}=1/3\) and \({\overline{\uplambda }}=1/2\). In each subplot the upper two plots show \(\log _{10}({\widehat{\upsigma }}^\prime )\) against \(\log _{10}(\omega )\), the lower row displays the behaviour of the derivative \({\mathrm d} \log _{10}({\widehat{\upsigma }}^\prime (\omega ))/{\mathrm d}\log _{10}(\omega )\) against \(\log _{10}(\omega )\). The numerical value of the derivatives may be viewed as a local scaling exponent. The abscissa in all plots is \(\log _{10}(\omega )\). The frequency is dimensionless and given in units of the relaxation frequency \(\omega _\mathbbm {p}\) defined in Eq. (9).

Fig. 4
figure 4

Analysis of the real part of the frequency dependent normalized conductivity tensor with eigenvalues \({\widehat{\upsigma }}^\prime _j(\omega )\) (\(j=x,y,z\)). In all subfigures frequency is dimensionless and given in units of the relaxation frequency \(\omega _\mathbbm {p}\) defined in Eq. (9). The six subfigures are labeled (a), (d), (g), (j), (e), and (i) to indicate the corresponding parameter quadruple \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) shown in the central subfigure of Fig. 3. Each subfigure consists of four subfigures. The two columns show plots for \({\overline{\uplambda }} = 0.35\) (left) resp. \({\overline{\uplambda }} =0.55\) (right) as indicated. In each subplot the upper two plots show \(\log _{10}({\widehat{\upsigma }}^\prime )\) against \(\log _{10}(\omega )\), the lower two show the local exponent \({\mathrm d}\log _{10}({\widehat{\upsigma }}^\prime )/{\mathrm d} \log _{10}(\omega )\) versus \(\log _{10}(\omega )\), i.e. the derivative of the upper plot. The line styles for \({\widehat{\upsigma }}^\prime _x\) (solid), \({\widehat{\upsigma }}^\prime _y\) (dashed), and \({\widehat{\upsigma }}^\prime _y\) (dash-dotted) are the same as in Fig. 3. The porosity is \({\overline{\phi }} = 0.1\) in all cases

The porosity is \({\overline{\phi }} = 0.1\) in all cases.

10.3.1 General observations

Two different sources can be identified which cause dispersion of \({\widehat{\upsigma }}^\prime \). The first source is the dispersion of the components \(P_j\) and \(M_j\), defined in Eq. (97d) and Eq. (97e). It is the dispersion in the response of a single local geometry. The second source causing dispersion arises from disorder and the percolation limit.

The spherical case (a) is shown for reference in the upper left corner of Fig. 4. In this case the dispersion is weak and it has a limited range in frequency.

The ordering \({\widehat{\upsigma }}^\prime _x(\omega ) \ge {\widehat{\upsigma }}^\prime _y(\omega ) \ge {\widehat{\upsigma }}^\prime _z(\omega )\) holds in all subfigures as for dc-conductivities.

In general, the dispersion decreases for increasing values of \({\overline{\uplambda }}\). Contrary to the limit \(\omega = 0\), the values of \({\widehat{\upsigma }}^\prime \) at high frequencies are barely affected by \({\overline{\uplambda }}\).

10.3.2 Flat pores \({{\overline{e}}}_\mathbbm {p} \approx 0\)

The case (d) from Fig. 3 is shown in the upper right corner of Fig. 4. Clearly, \({\widehat{\upsigma }}^\prime _z \ll {\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _z \ll {\widehat{\upsigma }}^\prime _y\) holds for all frequencies. As in the case \(\omega = 0\), this results from \(P_z \rightarrow 0\) for \({{\overline{e}}}_\mathbbm {p} \approx 0\). For \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\), there is one dispersive frequency range, which arises from the percolation transition at \({\overline{\uplambda }} = 1/3\). For \({\widehat{\upsigma }}^\prime _z\), there is an additional dispersion at higher frequencies. It arises from the dispersion of \(P_z\).

10.3.3 Cylindrical pores \({{\overline{d}}}_\mathbbm {p} \approx 0\)

The case (j) from Fig. 3 is placed on the right side of the middle row in Fig. 4. As discussed in section 10.2.3, \({\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) approach a percolation transition at \({\overline{\uplambda }} = 1/2\). This yields two effects for the dispersion of these components at \({\overline{\uplambda }} = 0.35\). First, the dispersion is much larger than the dispersion of \({\widehat{\upsigma }}^\prime _x\), and second, it occurs in a much broader frequency range. Both observations can be explained by the fact that the percolation limits at \({\overline{\uplambda }} = 1/3\) and \({\overline{\uplambda }} = 1/2\) lead to dispersion in two different frequency ranges which overlap. This effect cannot be observed at \({\overline{\uplambda }} = 0.55 > 1/2\).

10.3.4 Completely flat case

The case (g) from Fig. 3 is shown on the middle left side of Fig. 4. The dispersion of \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\) is similar to case (j), because \({\widehat{\upsigma }}^\prime _x\) approaches the percolation limit at \({\overline{\uplambda }} = 1/3\) and \({\widehat{\upsigma }}^\prime _y\) at \({\overline{\uplambda }} = 1/2\). The behavior of \({\widehat{\upsigma }}^\prime _z\) can be qualitatively understood as a combination of the cylindrical (j) and the circular case (d). Compared to case (j), there is an additional dispersion at high frequencies which arises from the dispersion of \(P_z\) for \({{\overline{e}}}_\mathbbm {p} \approx 0\).

10.3.5 Flat matrix \({{\overline{e}}}_\mathbbm {m} \approx 0\)

The case (e) from Fig. 3 is shown in the lower left corner of Fig. 4. At first sight the dispersion looks similar to the flat case (g). However, the additional dispersion at high frequencies for \({\widehat{\upsigma }}^\prime _z\) results from the percolation limit at \({\overline{\uplambda }} = 1\), and not from the dispersion of \(P_z\) for \({{\overline{e}}}_\mathbbm {p} \approx 0\). This can be concluded from the fact that this dispersion barely changes between \({\overline{\uplambda }} = 0.35\) and \({\overline{\uplambda }} = 0.55\).

10.3.6 Mixed case \({{\overline{d}}}_\mathbbm {p} \approx 0\), \({{\overline{e}}}_\mathbbm {m} \approx 0\)

The case (i) from Fig. 3 is shown in the lower right corner of Fig. 4. Because \({{\overline{d}}}_\mathbbm {p} \approx 0\), the conductivity eigenvalue \({\widehat{\upsigma }}^\prime _y\) approaches the percolation limit at \({\overline{\uplambda }} = 1/2\), and \({\widehat{\upsigma }}^\prime _z\) approaches the percolation limit at \({\overline{\uplambda }} = 1\) because \({{\overline{e}}}_\mathbbm {m} \approx 0\) holds. This means that the three different eigenvalues of \({\widehat{\upsigma }}^\prime \) approach three different percolation limits. From this follows that the components are dispersive over different broad frequency ranges. At \({\overline{\uplambda }} = 0.35\), the dispersion at low frequencies results from the transition at \({\overline{\uplambda }} = 1/3\), the dispersion at intermediate frequencies from the transition at \({\overline{\uplambda }} = 1/2\), and the dispersion at high frequencies from the transition at \({\overline{\uplambda }} = 1\).

Fig. 5
figure 5

Numerical solution of Eq. (109) for \({\widehat{\upsigma }}^\prime _x({\overline{\phi }})\) (solid curve), \({\widehat{\upsigma }}^\prime _y({\overline{\phi }})\) (dashed curve) and \({\widehat{\upsigma }}^\prime _z({\overline{\phi }})\) (dash-dotted curve) as functions of average porosity \({\overline{\phi }}\) for anisotropy parameters \({{\overline{d}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{e}}}_\mathbbm {m}=0.58\) and local percolation probability \(\uplambda ({\overline{\phi }})\) as in Fig. 2. The percolation threshold is at \({\overline{\phi }}_c=1/27\approx 3.7\%\) as in Figure 2. Experimental data for shaly sands, normalized as in Eq. (108), are displayed as symbols. Square symbols are the data from Table 1, while circles are from Table 3 in Waxman and Smits [95]

11 Application to experiment

A simple application of the new approximations to experiment can be made for anisotropic shaly quartz sands. It has recently been suggested by Nguyen, Vu, and Vu [94] that the electrical conductivities of shaly sands measured by Waxman and Smits [95] and Hill and Milburn [96] are anisotropic. Anisotropy ratios as large as 5 between normal and transversal conductivities are suggested in Nguyen, Vu, and Vu [94], Table 1. For details on the samples and their description the reader is referred to Waxman and Smits [95], p. 113. The electrical conductivity of the anisotropic shaly sands was measured using an impedance bridge for only one particular orientation. The data with an accuracy of 0.1 percent will be viewed here as resulting from randomly orienting the sample. Figure 5 shows the logarithm of the measured d.c. electrical conductivities normalized according to Eq. (108). Square symbols are the data reported in Table 1, while circles are from Table 3 in Waxman and Smits [95].

Detailed sample specific information for the anisotropy of pore \(\mathbbm {p}\) or matrix \(\mathbbm {m}\) in the shaly sands are not available from the publication [95]. To apply the theory nevertheless, it will be assumed for simplicity that \({{\overline{d}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{e}}}_\mathbbm {m}\) holds true. For the local percolation probability \(\uplambda ({\overline{\phi }})\) the function underlying Fig. 2 is assumed again, so that the percolation threshold occurs again at porosity \({\overline{\phi }}_c=1/27\approx 3.7\%\). The numerical solution of Eq. (109) yields the three eigenvalues \({\widehat{\upsigma }}^\prime _x({\overline{\phi }}), {\widehat{\upsigma }}^\prime _y({\overline{\phi }})\) and \({\widehat{\upsigma }}^\prime _z({\overline{\phi }})\) of the conductivity tensor as functions of average porosity. The three curves for \({{\overline{d}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{e}}}_\mathbbm {m}=0.58\) are displayed in Fig. 5.

The figure suggests that the scatter in the experimental data may indeed arise from an anisotropy of the samples. In conclusion, the theoretical assumptions underlying Eq. (109) and the experimental indications, that the measured conductivities of shaly sands are anisotropic, can be consistently quantified with the help of the approximations developed in this paper.

12 Conclusion

Generalized effective medium approximations for anisotropic media have been derived from local porosity theory. A purely geometric and experimentally measurable characterization of the anisotropy in (local) connectivities was introduced as a prerequisite for the generalized theory. The resulting anisotropic local porosity theory is based on the familiar local geometry distributions from the isotropic case. These distributions characterize the complex pore-space geometry. The generalized theory contains isotropic local porosity theory as a special case. The local geometry distributions characterizing local anisotropic connectivity are obtained from local percolating directions. Local percolating directions were defined based on the generalized percolation criterion introduced in Eq. (21). The geometric quantities enter directly into the generalized anisotropic self-consistent effective medium equation. No adjustable parameters appear in the generalized anisotropic local porosity theory. Several hitherto unknown generalized effective medium approximation formulae are reported. Among them the analytical results (125) and (132) for the cylindrical cases are of greatest interest. The generalized theory recovers previously known anisotropic effective medium approximations in the literature as special cases. The practical applicability of the theoretical predictions has been tested by a simple comparison with experiment.