# Nonlocal Effective Electromagnetic Wave Characteristics of Composite Media:
Beyond the Quasistatic Regime

Salvatore Torquato<EMAIL_ADDRESS>http://torquato.princeton.edu
Department of Physics, Princeton University, Princeton, New Jersey 08544, USA
Department of Chemistry, Princeton University, Princeton, New Jersey 08544,
USA  Princeton Institute for the Science and Technology of Materials,
Princeton University, Princeton, New Jersey 08544, USA  Program in Applied and
Computational Mathematics, Princeton University, Princeton, New Jersey 08544,
USA  Jaeuk Kim  Department of Physics, Princeton University, Princeton, New
Jersey 08544, USA

###### Abstract

We derive exact nonlocal homogenized constitutive relations for the effective
electromagnetic wave properties of disordered two-phase composites and
metamaterials from first principles. This exact formalism enables us to extend
the long-wavelength limitations of conventional homogenization estimates of
the effective dynamic dielectric constant tensor
$\boldsymbol{\varepsilon}_{e}({\bf k}_{I},\omega)$ for arbitrary
microstructures so that it can capture spatial dispersion well beyond the
quasistatic regime (where $\omega$ and ${\bf k}_{I}$ are the frequency and
wave vector of the incident radiation). We accomplish this task by deriving
nonlocal strong-contrast expansions that exactly account for complete
microstructural information (infinite set of $n$-point correlation functions)
and hence multiple scattering to all orders for the range of wave numbers for
which our extended homogenization theory applies, i.e., $0\leq|{\bf
k}_{I}|\ell\lesssim 1$ (where $\ell$ is a characteristic heterogeneity length
scale). Because of the fast-convergence properties of such expansions, their
lower-order truncations yield accurate closed-form approximate formulas for
$\boldsymbol{\varepsilon}_{e}({\bf k}_{I},\omega)$ that apply for a wide class
of microstructures. These nonlocal formulas are resummed representations of
the strong-contrast expansions that still accurately capture multiple
scattering to all orders via the microstructural information embodied in the
spectral density, which is easy to compute for any composite. The accuracy of
these microstructure-dependent approximations is validated by comparison to
full-waveform simulation computations for both 2D and 3D ordered and
disordered models of composite media. Thus, our closed-form formulas enable
one to predict accurately and efficiently the effective wave characteristics
well beyond the quasistatic regime for a wide class of composite
microstructures without having to perform full-blown simulations. We find that
disordered hyperuniform media are generally less lossy than their
nonhyperuniform counterparts. We also show that certain disordered
hyperuniform particulate composites exhibit novel wave characteristics,
including the capacity to act as low-pass filters that transmit waves
“isotropically” up to a selected wave number and refractive indices that
abruptly change over a narrow range of wave numbers. Our results demonstrate
that one can design the effective wave characteristics of a disordered
composite by engineering the microstructure to possess tailored spatial
correlations at prescribed length scales. Thus, our findings can accelerate
the discovery of novel electromagnetic composites.

##  I Introduction

The theoretical problem of estimating the effective properties of multiphase
composite media is an outstanding one and dates back to work by some of the
luminaries of science, including Maxwell Maxwell (1873), Lord Rayleigh Lord
Rayleigh (1892), and Einstein Einstein (1906). The preponderance of previous
theoretical studies have focused on the determination of static effective
properties (e.g., dielectric constant, elastic moduli and fluid permeability)
using a variety of methods, including approximation schemes Maxwell (1873);
Bruggeman (1935); Brinkman (1949); Budiansky (1965), bounding techniques
Prager (1961); Hashin and Shtrikman (1962); Beran (1965); Kohn and Lipton
(1988); Torquato (2002); Milton (2002) and exact series-expansion procedures
Brown (1955); Felderhof  _et al._ (1982); Sen and Torquato (1989); Torquato
(1997). The latter set of investigations teaches us that an exact
determination of an effective property, given the phase properties of the
composite, generally requires an infinite set of correlation functions that
characterizes the composite microstructure.

Our focus in this paper is the determination of the effective dynamic
dielectric constant tensor $\boldsymbol{\varepsilon}_{e}({\bf k}_{I},\omega)$
of a two-phase dielectric composite, which depends on the frequency $\omega$
or wave vector ${\bf k}_{I}$ of the incident radiation Sheng (1995); Sihvola
(1999). From this effective property, one can determine the corresponding
effective wave speed $c_{e}$ and attenuation coefficient $\gamma_{e}$. The
preponderance of previous homogenization theories of
$\boldsymbol{\varepsilon}_{e}({\bf k}_{I},\omega)$ apply only in the
quasistatic or long-wavelength regime Keller (1964); Frisch (1968); Tsang and
Kong (1981); Jing  _et al._ (1992)111 A notable exception is the generalized
coherent potential approximation for particle suspensions presented in Ref.
Jing  _et al._ (1992). However, this is not a closed-form formula for the
effective dielectric constant and requires, as input, certain numerical
simulations of the electric fields., i.e., applicable when $|{\bf
k}_{I}|\ell\ll 1$, where $\ell$ is a characteristic heterogeneity length
scale. This spectral range is the realm of nonresonant dielectric behavior.
Virtually all earlier formulas $\boldsymbol{\varepsilon}_{e}({\bf
k}_{I},\omega)$ apply to very special microstructures, namely, dielectric
scatterers that are well-defined inclusions, e.g., nonoverlapping spheres in a
matrix (see Fig. 1). Examples of such popular closed-form approximation
formulas devised for spherical scatterers in a matrix include the Maxwell-
Garnett Garnett (1904); Ruppin (2000) and quasicrystalline Lax (1952); Sihvola
(1999); Ao and Kong (2002) estimates, among others Sheng (1995).

Figure 1: Left panel: The preponderance of previous theoretical treatments of
the effective dynamic dielectric constant are local in nature and restricted
to dispersions of well-defined dielectric scatterers (inclusions) in a matrix,
as illustrated in the leftmost panel. In contrast to the nonlocal strong-
contrast formalism presented here, earlier studies cannot treat the more
general two-phase microstructures shown in the middle panel (spinodal
decomposition pattern Ma and Torquato (2017)) and rightmost panel (Debye
random medium Yeong and Torquato (1998); Ma and Torquato (2020)), both of
which have “phase-inversion” symmetry Torquato (2002). Our formalism can in
principle treat two-phase media of arbitrary microstructures.

In the present investigation, we derive nonlocal homogenized constitutive
relations from Maxwell’s equations to obtain exact expressions for the
effective dynamic dielectric constant tensor
$\boldsymbol{\varepsilon}_{e}({\bf k}_{I},\omega)$ of a macroscopically
anisotropic two-phase medium of arbitrary microstructure that is valid well
beyond the quasistatic regime, i.e., from the infinite-wavelength limit down
to intermediate wavelengths ($0\leq|{\bf k}_{I}|\ell\lesssim 1$). This task is
accomplished by extending the strong-contrast expansion formalism, which has
been used in the past exclusively for the static limit Torquato (2002) and
quasistatic regime Rechtsman and Torquato (2008), and establishing that the
resulting homogenized constitutive relations are nonlocal in space (Sec. III),
i.e., the average polarization field at position $\bf x$ depends on the
average electric field at other positions around $\bf x$. (Such nonlocal
relations are well known in the context of crystal optics in order to account
for “spatial dispersion,” i.e., the dependence of dielectric properties on a
wave vector Agranovich and Ginzburg (1984).) The terms of the strong-contrast
expansion are explicitly given in terms of integrals over products of Green’s
functions and the $n$-point correlation functions of the random two-phase
medium (defined in Sec. II.1) to infinite order. This representation exactly
treats multiple scattering to all orders for the range of wave numbers for
which our extended homogenization theory applies, i.e., $0\leq|{\bf
k}_{I}|\ell\lesssim 1$.

It is noteworthy that the strong-contrast formalism is a significant departure
from standard multiple-scattering theory Frisch (1968); Sheng (1995); Tsang
(2001); Cazé and Schotland (2015), as highlighted in Sec. III. Moreover, as we
show there, our strong-contrast formalism has a variety of “tuning knobs” that
enable one to obtain distinctly different expansions and approximations
designed for different classes of microstructures.

Because of the fast-convergence properties of strong-contrast expansions,
elaborated in Sec. III.2, their lower-order truncations yield accurate closed-
form approximate formulas for the effective dielectric constant that apply for
a wide class of microstructures over the aforementioned broad range of
incident wavelengths, volume fractions, and contrast ratios (Sec. IV). Thus,
we are able to accurately account for multiple scattering in the resonant
realm (e.g., Bragg diffraction for periodic media), in contrast to the
Maxwell-Garnett and quasicrystalline approximations, which are known to break
down in this spectral range. These nonlocal strong-contrast formulas can be
regarded as approximate resummations of the expansions that still accurately
capture multiple-scattering effects to all orders via the nonlocal attenuation
function $F({\bf Q})$ (Sec. VI). The key quantity $F({\bf Q})$ is a functional
of the spectral density $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left({\bf
Q}\right)}}$ (Sec. II.2), which is straightforward to determine for general
microstructures either theoretically, computationally or via scattering
experiments. We employ precise full-waveform simulation methods (Sec. VII) to
show that these microstructure-dependent approximations are accurate for both
two-dimensional (2D) and three-dimensional (3D) ordered and disordered models
of particulate composite media (Sec. VIII). This validation means that they
can be used to predict accurately the effective wave characteristics well
beyond the quasistatic regime for a wide class of composite microstructures
(Sec. IX) without having to perform full-blown simulations. This broad
microstructure class includes particulate media consisting of identical or
polydisperse particles of general shape (ellipsoids, cubes, cylinders,
polyhedra) that may or not overlap, and cellular networks as well as media
without well-defined inclusions (see Sec. III.2 for details). Thus, our
nonlocal formulas can be employed to accelerate the discovery of novel
electromagnetic composites by appropriate tailoring of the spectral densities
Torquato (2016); Chen and Torquato (2018) and then generating the
microstructures that achieve them Chen and Torquato (2018).

Although our strong-contrast formulas for the effective dynamic dielectric
constant apply to periodic two-phase media, the primary applications that we
have in mind are correlated disordered microstructures because they can
provide advantages over periodic ones with high crystallographic symmetries
López (2018); Yu  _et al._ (2015) which include perfect isotropy and
robustness against defects Florescu  _et al._ (2013); Man  _et al._ (2013). We
are interested in both “garden-variety” models of disordered two-phase media
Torquato (2002) as well as exotic hyperuniform forms Torquato and Stillinger
(2003); Zachary and Torquato (2009); Torquato (2018a). Hyperuniform two-phase
systems are characterized by an anomalous suppression of volume-fraction
fluctuations in the infinite-wavelength limit Torquato and Stillinger (2003);
Zachary and Torquato (2009); Torquato (2018a), i.e., the spectral density
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left({\bf Q}\right)}}$ obeys the
condition

$\lim_{\absolutevalue{{\bf Q}}\to
0}\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left({\bf Q}\right)}}=0.$ (1)

Such two-phase media encompass all periodic systems, many quasiperiodic media
and exotic disordered ones; see Ref. Torquato (2018a) and references therein.
Disordered hyperuniform systems lie between liquids and crystals; they are
like liquids in that they are statistically isotropic without any Bragg peaks,
and yet behave like crystals in the manner in which they suppress the large-
scale density fluctuations Torquato and Stillinger (2003); Zachary and
Torquato (2009); Torquato (2018a). Hyperuniform systems have attracted great
attention over the last decade because of their deep connections to a wide
range of topics that arise in physics Torquato  _et al._ (2015); Zhang  _et
al._ (2016a); Hexner  _et al._ (2017); Ricouvier  _et al._ (2017); Oğuz  _et
al._ (2017); Ma and Torquato (2017); López (2018); Yu  _et al._ (2018); Wang
_et al._ (2018); Torquato (2018a); Lei and Ni (2019); Gorsky  _et al._ (2019);
Klatt  _et al._ (2020), materials science Ma  _et al._ (2016); Xu  _et al._
(2017); Torquato and Chen (2018); Chen and Torquato (2018); Kim and Torquato
(2019a), mathematics Ghosh and Lebowitz (2018); Brauchart  _et al._ (2019);
Torquato  _et al._ (2019), and biology Jiao  _et al._ (2014); Torquato (2018a)
as well as for their emerging technological importance in the case of the
disordered varieties Florescu  _et al._ (2009); Man  _et al._ (2013); Ma  _et
al._ (2016); Leseur  _et al._ (2016); Froufe-Pérez  _et al._ (2017); Klatt and
Torquato (2018); Zhang  _et al._ (2018); Torquato (2018a); Gorsky  _et al._
(2019).

We apply our strong-contrast formulas to predict the real and imaginary parts
of the effective dielectric constant for model microstructures with typical
disorder (nonhyperuniform) as well as those with exotic hyperuniform disorder
(Sec. V). We are particularly interested in exploring the dielectric
properties of a special class of hyperuniform composites called disordered
stealthy hyperuniform media, which are defined to be those that possess zero-
scattering intensity for a set of wave vectors around the origin Uche  _et
al._ (2004); Batten  _et al._ (2008); Torquato  _et al._ (2015); Zhang  _et
al._ (2015); Chen and Torquato (2018). Such materials have recently been shown
to be endowed with novel optical, acoustic, mechanical, and transport
properties Rechtsman and Torquato (2008); Leseur  _et al._ (2016); Zhang  _et
al._ (2016b); Degl’Innocenti  _et al._ (2016); Gkantzounis  _et al._ (2017);
Xu  _et al._ (2017); Torquato and Chen (2018); Kim and Torquato (2020). Among
other findings, we show that disordered hyperuniform media are generally less
lossy than their nonhyperuniform counterparts. We also demonstrate that
disordered stealthy hyperuniform particulate composites exhibit singular wave
characteristics, including the capacity to act as low-pass filters that
transmit waves “isotropically” up to a selected wave number. They also can be
engineered to exhibit refractive indices that abruptly change over a narrow
range of wave numbers by tuning the spectral density. Our results demonstrate
that one can design the effective wave characteristics of a disordered
composite, hyperuniform or not, by engineering the microstructure to possess
tailored spatial correlations at prescribed length scales.

##  II Background

###  II.1 $n$-point correlation functions

A two-phase random medium is a domain of space
$\mathcal{V}\subseteq\mathbb{R}^{d}$ that is partitioned into two disjoint
regions that make up $\mathcal{V}$: a phase 1 region $\mathcal{V}_{1}$ of
volume fraction $\phi_{1}$ and a phase 2 region $\mathcal{V}_{2}$ of volume
fraction $\phi_{2}$ Torquato (2002). The phase indicator function
$\mathcal{I}^{(i)}(\mathbf{x})$ of phase $i$ for a given realization is
defined as

$\mathcal{I}^{(i)}(\mathbf{x})=\left\\{\begin{aligned}
&1,&\mathbf{x}\in\mathcal{V}_{i},\\\
&0,&\mathbf{x}\notin\mathcal{V}_{i}.\end{aligned}\right.$ (2)

The $n$-point correlation function $S_{n}^{(i)}$ for phase $i$ is defined by
Torquato and Stell (1982); Torquato (2002):

$\displaystyle
S_{n}^{(i)}(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})=\left\langle\prod_{j=1}^{n}{\cal
I}^{(i)}(\mathbf{x}_{j})\right\rangle,$ (3)

where angular brackets denote an ensemble average over realizations. The
function $S_{n}^{(i)}(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$ has a
probabilistic interpretation: It gives the probability of finding the ends of
the vectors $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ all in phase $i$. For
statistically homogeneous media,
$S_{n}^{(i)}(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$ is translationally
invariant and, in particular, the one-point function is position independent,
i.e., $S_{1}^{(i)}(\mathbf{x_{1}})=\phi_{i}$.

###  II.2 Two-point statistics

For statistically homogeneous media, the two-point correlation function for
phase 2 is simply related to that for phase 1 via the expression
$S^{(2)}_{2}({\bf r})=S^{(1)}_{2}({\bf r})-2\phi_{1}+1$, and hence the
autocovariance function is given by

$\chi_{{}_{V}}({\bf r})\equiv S^{(1)}_{2}({\bf
r})-{\phi_{1}}^{2}=S^{(2)}_{2}({\bf r})-{\phi_{2}}^{2},$ (4)

which we see is the same for phase 1 and phase 2. Thus, $\chi_{{}_{V}}({\bf
r}=0)=\phi_{1}\phi_{2}$ and, assuming the medium possesses no long-range
order, $\lim_{|{\bf r}|\rightarrow\infty}\chi_{{}_{V}}({\bf r})=0$. For
statistically homogeneous and isotropic media, ${\chi}_{{}_{V}}({\bf r})$
depends only on the magnitude of its argument $r=|\bf r|$, and hence is a
radial function. In such cases, its slope at the origin is directly related to
the specific surface $s$ (interface area per unit volume); i.e.,
asymptotically, we have Torquato (2002)

$\chi_{{}_{V}}({\bf r})=\phi_{1}\phi_{2}-\beta(d)s\;|{\bf r}|+{\cal O}(|{\bf
r}|^{2}),$ (5)

where

$\beta(d)=\frac{\Gamma(d/2)}{2\sqrt{\pi}\Gamma((d+1)/2)},$ (6)

and $\Gamma(x)$ is the gamma function.

The nonnegative spectral density ${\tilde{\chi}}_{{}_{V}}({\bf Q})$, which is
proportional to scattering intensity Debye  _et al._ (1957), is the Fourier
transform of $\chi_{{}_{V}}({\bf r})$, i.e.,

${\tilde{\chi}}_{{}_{V}}({\bf Q})=\int_{\mathbb{R}^{d}}\chi_{{}_{V}}({\bf
r})\,e^{-i{\bf Q\cdot r}}{\rm d}{\bf r}\geq 0,\qquad\mbox{for all}\;{\bf Q},$
(7)

where ${\bf Q}$ represents the momentum-transfer wave vector. For
statistically homogeneous media, the spectral density must obey the following
sum rule Torquato (2020):

$\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}{\tilde{\chi}}_{{}_{V}}({\bf
Q})\,\differential{\bf Q}=\chi_{{}_{V}}({\bf r}=0)=\phi_{1}\phi_{2}.$ (8)

For isotropic media, the spectral density depends only on $Q=|{\bf Q}|$ and,
as a consequence of Eq. (5), its large-$k$ behavior is controlled by the
following power-law form:

${\tilde{\chi}}_{{}_{V}}({\bf Q})\sim\frac{\gamma(d)\,s}{Q^{d+1}},\qquad
Q\rightarrow\infty,$ (9)

where

$\gamma(d)=2^{d}\,\pi^{(d-1)/2}\,\Gamma((d+1)/2)\,\beta(d)$ (10)

is a $d$-dimensional constant and $\beta(d)$ is given by Eq. (6).

###  II.3 Packings

We call a packing in $\mathbb{R}^{d}$ a collection of nonoverlapping particles
Torquato (2018b). In the case of a packing of identical spheres of radius $a$
at number density $\rho$, the spectral density
$\tilde{\chi}_{{}_{V}}(\textbf{Q})$ is directly related to the structure
factor $S({\bf Q})$ of the sphere centers (Torquato, 2002, 2016):

$\displaystyle\tilde{\chi}_{{}_{V}}(\textbf{Q})=\phi_{2}\,\tilde{\alpha}_{2}(Q;a)\,S(\textbf{Q}),$
(11)

where

$\displaystyle\tilde{\alpha}_{2}(Q;a)=\frac{1}{v_{1}(a)}\left(\frac{2\pi
a}{Q}\right)^{d}\,J^{2}_{d/2}(Qa),$ (12)

$\mathinner{J_{\nu}\mathopen{\left(x\right)}}$ is the Bessel function of the
first kind of order $\nu$, $\phi_{2}=\rho v_{1}(a)$ is the packing fraction
(fraction of space covered by the spheres), and

$v_{1}(a)=\frac{\pi^{d/2}a^{d}}{\Gamma(1+d/2)}$ (13)

is the $d$-dimensional volume of a sphere of radius $a$.

###  II.4 Hyperuniformity and volume-fraction fluctuations

Originally introduced in the context of point configurations Torquato and
Stillinger (2003), the hyperuniformity concept was generalized to treat two-
phase media Zachary and Torquato (2009). Here the phase volume fraction
fluctuates within a spherical window of radius $R$, which can be characterized
by the local volume-fraction variance $\sigma_{{}_{V}}^{2}(R)$. This variance
is directly related to integrals involving either the autocovariance function
$\chi_{{}_{V}}({\bf r})$ Lu and Torquato (1990) or the spectral density
$\tilde{\chi}_{{}_{V}}({\bf Q})$ Zachary and Torquato (2009).

For typical disordered two-phase media, the variance $\sigma_{{}_{V}}^{2}(R)$
for large $R$ goes to zero like $R^{-d}$ Lu and Torquato (1990); Quintanilla
and Torquato (1997); Torquato (2002). However, for hyperuniform disordered
media, $\sigma_{{}_{V}}^{2}(R)$ goes to zero asymptotically more rapidly than
the inverse of the window volume, i.e., faster than $R^{-d}$, which is
equivalent to the condition (1) on the spectral density. Stealthy hyperuniform
two-phase media are a subclass of hyperuniform systems in which
$\tilde{\chi}_{{}_{V}}(\mathbf{Q})$ is zero for a range of wave vectors around
the origin, i.e.,

$\tilde{\chi}_{{}_{V}}(\mathbf{Q})=0\qquad\mbox{for}\;0\leq|{\bf Q}|\leq
Q_{\text{U}},$ (14)

where $Q_{\text{U}}$ is some positive number.

As in the case of hyperuniform point configurations Torquato and Stillinger
(2003); Zachary and Torquato (2009, 2011); Torquato (2018a), there are three
different scaling regimes (classes) that describe the associated large-$R$
behaviors of the volume-fraction variance when the spectral density goes to
zero as a power law ${\tilde{\chi}}_{{}_{V}}({\bf Q})\sim|{\bf Q}|^{\alpha}$
as $|\bf Q|\to 0$:

$\displaystyle\sigma^{2}_{{}_{V}}(R)\sim\begin{cases}R^{-(d+1)},\quad\quad\quad\alpha>1\qquad&\text{(Class
I)},\\\ R^{-(d+1)}\ln R,\quad\alpha=1\qquad&\text{(Class II)},\\\
R^{-(d+\alpha)},\quad 0<\alpha<1\qquad&\text{(Class III)},\end{cases}$ (15)

where the exponent $\alpha$ is a positive constant. Thus, the characteristics
length of the representative elementary volume for a hyperuniform medium will
depend on the hyperuniformity class (scaling). Class I is the strongest form
of hyperuniformity, which includes all perfect periodic packings as well as
some disordered packings, such as disordered stealthy packings described in
Sec. V.4.

###  II.5 Popular effective-medium approximations

Here we explicitly state the specific functional forms of an extended Maxwell-
Garnett approximation and quasicrystalline approximation for the effective
dynamic dielectric constant $\varepsilon_{e}(k_{1})$ of isotropic media
composed of identical spheres of dielectric constant $\varepsilon_{2}$
embedded in a matrix phase of dielectric constant $\varepsilon_{1}$. In Sec.
VIII, we compare the predictions of these formulas to those of our nonlocal
approximations. The small-wave-number expansions of these popular
approximations are provided in the Supplementary Material SM . There we also
provide the corresponding asymptotic behaviors of our strong-contrast
approximations.

####  II.5.1 Maxwell-Garnett approximation

Maxwell-Garnett approximations (MGAs) Garnett (1904); Ruppin (2000) are
derived by substituting the dielectric polarizability of a single dielectric
sphere into the Clausius-Mossotti equation Ruppin (2000), which consequently
ignores the spatial correlations between the particles. In three dimensions,
we utilize the following extended MGA that makes use of the exact electric
dipole polarizability $\mathinner{\alpha_{e}\mathopen{\left(k_{1}\right)}}$ of
a single dielectric sphere of radius $a$ Doyle (1989):

$\displaystyle\frac{\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}-\varepsilon_{1}}{\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}+2\varepsilon_{1}}=$
$\displaystyle\phi_{2}\mathinner{\alpha_{e}\mathopen{\left(k_{1}\right)}}/a^{3},$
(16)

where

$\mathinner{\alpha_{e}\mathopen{\left(k_{1}\right)}}=\frac{3i}{2{k_{1}}^{3}}\frac{m\mathinner{\psi_{1}\mathopen{\left(mk_{1}a\right)}}\mathinner{\psi_{1}^{\prime}\mathopen{\left(k_{1}a\right)}}-\mathinner{\psi_{1}\mathopen{\left(k_{1}a\right)}}\mathinner{\psi_{1}^{\prime}\mathopen{\left(mk_{1}a\right)}}}{m\mathinner{\psi_{1}\mathopen{\left(mk_{1}a\right)}}\mathinner{\xi_{1}^{\prime}\mathopen{\left(k_{1}a\right)}}-\mathinner{\xi_{1}\mathopen{\left(k_{1}a\right)}}\mathinner{\psi_{1}^{\prime}\mathopen{\left(mk_{1}a\right)}}},$

$m\equiv\sqrt{\varepsilon_{2}/\varepsilon_{1}}$,
$\mathinner{\psi_{1}\mathopen{\left(x\right)}}\equiv
x\mathinner{j_{1}\mathopen{\left(x\right)}}$,
$\mathinner{\xi_{1}\mathopen{\left(x\right)}}\equiv
x\mathinner{h_{1}^{(1)}\mathopen{\left(x\right)}}$, the prime symbol
$(^{\prime})$ denotes the derivative of a function,
$\mathinner{h_{1}^{(1)}\mathopen{\left(x\right)}}$ is the spherical Hankel
function of the first kind of order 1, and
$\mathinner{j_{1}\mathopen{\left(x\right)}}$ is the spherical Bessel function
of the first kind of order 1.

The 2D analog of Eq. (16) can be obtained by using the dynamic dielectric
polarizability $\alpha_{e}$ of a dielectric cylinder of radius $a$ given in
Ref. Silveirinha and Engheta (2007):

$\displaystyle\frac{\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}-\varepsilon_{1}}{\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}+\varepsilon_{1}}=$
$\displaystyle\frac{\phi_{2}}{2\pi}\mathinner{\alpha_{e}\mathopen{\left(k_{1}\right)}}/a^{2},$
(17)

where

$\displaystyle\mathinner{\alpha_{e}\mathopen{\left(k_{1}\right)}}=$
$\displaystyle\frac{4(\varepsilon_{2}-\varepsilon_{1})}{i{k_{1}}^{2}m\varepsilon_{1}}\mathinner{J_{1}\mathopen{\left(mk_{1}a\right)}}\big{[}\mathinner{J_{1}^{\prime}\mathopen{\left(mk_{1}a\right)}}\mathinner{\mathcal{H}_{1}^{(1)}\mathopen{\left(k_{1}a\right)}}$
$\displaystyle-m\mathinner{J_{1}\mathopen{\left(mk_{1}a\right)}}\mathinner{{\mathcal{H}_{1}^{(1)}}^{\prime}\mathopen{\left(k_{1}a\right)}}\big{]}^{-1},$

where $\mathinner{\mathcal{H}_{\nu}^{(1)}\mathopen{\left(x\right)}}$ is the
Hankel function of the first kind of order $\nu$. Here, only transverse-
electric (TE) polarization is considered; i.e., the electric field is
perpendicular to the axis of the cylinder.

As with all MGA theories, formulas (16) and (17) neglect spatial correlations
between the particles and hence are only valid for low inclusion packing
fractions. In the static limit, Eqs. (16) and (17) reduce to the Hashin-
Shtrikman estimates $\varepsilon_{\mathrm{HS}}$; see relation (72).

####  II.5.2 Quasicrystalline approximations

The quasicrystalline approximation (QCA) for the quasistatic effective dynamic
dielectric constant $\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}$
employs the “effective” Green’s function of spherical scatterers up to the
level of the pair correlation function $\mathinner{g_{2}\mathopen{\left({\bf
r}\right)}}$ Lax (1952); Sihvola (1999); Ao and Kong (2002). However, the QCA
accounts only for the structure factor in the infinite-wavelength limit Ao and
Kong (2002) [i.e.,
$\mathinner{S\mathopen{\left(0\right)}}=1+\rho\int(\mathinner{g_{2}\mathopen{\left({\bf
r}\right)}}-1)\differential{{\bf r}}$] and consequently, spatial correlations
at finite wavelengths are ignored. The QCA for $d=3$ can be explicitly written
as follows Sihvola (1999):

$\displaystyle{\phi_{2}}^{2}\beta_{21}\quantity[\frac{\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}-\varepsilon_{1}}{\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}+2\varepsilon_{1}}]^{-1}=\phi_{2}-i\Bigg{\\{}\frac{2}{3}\phi_{2}\mathinner{S\mathopen{\left(0\right)}}(k_{1}a)^{3}$
$\displaystyle\times\quantity[1+i\frac{2}{3(1-\beta_{21}\phi_{2})}(k_{1}a)^{3}\mathinner{S\mathopen{\left(0\right)}}]^{-1}\Bigg{\\}}\beta_{21},$
(18)

where $\beta_{21}$ is defined in Eq. (33). Interestingly, the QCA predicts
that the hyperuniform composites [$S(0)=0$] will be transparent for all wave
numbers, which cannot be true for stealthy hyperuniform media [cf. Eq. (14)],
since the transparency interval must be finite; see Sec. VIII. Note that Eq.
(18) is the complex conjugate of the one given in Ref. Sihvola (1999) so that
it is consistent with the sign convention for the imaginary part
$\imaginary[\varepsilon_{e}]$ used here.

##  III Exact Strong-Contrast Expansions and Nonlocality

The original strong-contrast expansions for the effective dynamic dielectric
constant obtained by Rechtsman and Torquato Rechtsman and Torquato (2008) were
derived from homogenized constitutive relations that are local in space and
hence are strictly valid only in the quasistatic regime. In the present work,
we follow the general strong-contrast formalism of Torquato Torquato (2002)
that was devised for the purely static problem Torquato (2002) and show that
it naturally leads to exact homogenized constitutive relations for the
averaged fields that are nonlocal in space. The crucial consequences of this
development are exact expressions for the effective dynamic dielectric
constant tensor $\boldsymbol{\varepsilon}_{e}({\bf k}_{I})$ for a
macroscopically anisotropic medium of arbitrary microstructure into which a
plane wave of wave vector ${\bf k}_{I}$ is incident. These expressions for
$\boldsymbol{\varepsilon}_{e}({\bf k}_{I})$ are valid from the infinite-
wavelength limit down to wavelengths ($\lambda=2\pi/\absolutevalue{{\bf
k}_{I}}$) on the order of the heterogeneity length scale $\ell$. We then
briefly explain how our theory departs substantially from standard multiple-
scattering theory Frisch (1968); Sheng (1995); Tsang (2001); Cazé and
Schotland (2015). We explicitly show they necessarily require complete
microstructural information, as embodied in the infinite set of $n$-point
correlation functions (Sec. II.1) of the composite. We also describe the fast-
convergence properties of strong-contrast expansions and their consequences
for extracting accurate approximations for $\boldsymbol{\varepsilon}_{e}({\bf
k}_{I})$.

###  III.1 Strong-contrast expansions

Figure 2: (a) Schematic of a large $d$-dimensional ellipsoidal,
macroscopically anisotropic two-phase composite medium embedded in an infinite
reference phase of dielectric constant tensor $\boldsymbol{\varepsilon}_{I}$
(gray regions) under an applied electric field $\mathinner{{\bf
E}_{0}\mathopen{\left({\bf x}\right)}}=\tilde{{\bf
E}}_{0}\mathinner{\exp\mathopen{\left(i({\bf k}_{I}\cdot{\bf x}-\omega
t)\right)}}$ of a frequency $\omega$ and a wave vector ${\bf k}_{I}$ at
infinity. The wavelength $\lambda$ associated with the applied field can span
from the quasistatic regime ($2\pi\ell/\lambda\ll 1$) down to the
intermediate-wavelength regime ($2\pi\ell/\lambda\lesssim 1$), where $\ell$ is
a characteristic heterogeneity length scale. (b) After homogenization, the
same ellipsoid can be regarded as a homogeneous specimen with an effective
dielectric constant $\boldsymbol{\varepsilon_{e}}({{\bf k}_{I},\omega})$,
which depends on $\omega$ and ${\bf k}_{I}$. As noted in the main text, we
omit the $\omega$ dependence of $\boldsymbol{\varepsilon_{e}}$ because
(without loss of generally) we assume a linear dispersion relation between
$|{\bf k}_{I}|$ and $\omega$.

Consider a macroscopically large ellipsoidal two-phase statistically
homogeneous but anisotropic composite specimen in $\mathbb{R}^{d}$ embedded
inside an infinitely large reference phase $I$ with a dielectric constant
tensor $\boldsymbol{\varepsilon}_{I}$. The microstructure is perfectly
general, and it is assumed that a characteristic heterogeneity length scale
$\ell$ is much smaller than the specimen size, i.e., $\ell\ll L$. The shape of
this specimen is purposely chosen to be nonspherical since any rigorously
correct expression for the effective property must ultimately be independent
of the shape of the composite specimen in the infinite-volume limit. It is
assumed that the applied or incident electric field $\mathinner{{\bf
E}_{0}\mathopen{\left({\bf x}\right)}}$ is a plane wave of an angular
frequency $\omega$ and wave vector ${\bf k}_{I}$ in the reference phase, i.e.,

$\mathinner{{\bf E}_{0}\mathopen{\left({\bf x}\right)}}=\tilde{{\bf
E}}_{0}\mathinner{\exp\mathopen{\left(i({\bf k}_{I}\cdot{\bf x}-\omega
t)\right)}}.$ (19)

Our interest is the exact expression for the effective dynamic dielectric
constant tensor $\boldsymbol{\varepsilon}_{e}({\bf k}_{I},\omega)$. Without
loss of generality, we assume a linear dispersion relation in the reference
phase, i.e., $k_{I}\equiv\absolutevalue{{\bf
k}_{I}}=\sqrt{\varepsilon_{I}}\omega/c$, where $c$ is the speed of light in
vacuum, and thus we henceforth do not explicitly indicate the dependence of
functions on $\omega$. The composite is assumed to be nonmagnetic, implying
that the phase magnetic permeabilities are identical, i.e.,
$\mu_{1}=\mu_{2}=\mu_{0}$, where $\mu_{0}$ is the magnetic permeability of the
vacuum. For simplicity, we assume real-valued, frequency-independent isotropic
phase dielectric constants $\varepsilon_{1}$ and $\varepsilon_{2}$.
Nonetheless, the composite can be generally lossy (i.e.,
$\boldsymbol{\varepsilon}_{e}$ is complex-valued) due to scattering from the
inhomogeneities in the local dielectric constant. It is noteworthy that our
results can be straightforwardly extended to phase dielectric constants that
are complex-valued (dissipative media), but this is not done in the present
work.

Here we present a compact derivation of strong-contrast expansions. It follows
the general formalism of Torquato Torquato (2002) closely but departs from it
at certain key steps when establishing the nonlocality of the homogenized
constitutive relation. (A detailed derivation is given in the Supplementary
Material SM .) For simplicity, we take the reference phase $I$ to be phase $q$
(equal to 1 or 2). Under the aforementioned assumptions, the local electric
field $\mathinner{{\bf E}\mathopen{\left({\bf x}\right)}}$ solves the time-
harmonic Maxwell equation Rechtsman and Torquato (2008):

$\curl{\curl{\mathinner{{\bf E}\mathopen{\left({\bf
x}\right)}}}}-{k_{q}}^{2}\mathinner{{\bf E}\mathopen{\left({\bf
x}\right)}}=\quantity(\frac{\omega}{c})^{2}\mathinner{{\bf
P}\mathopen{\left({\bf x}\right)}},$ (20)

where $\mathinner{{\bf P}\mathopen{\left({\bf x}\right)}}$ is the polarization
field given by

$\mathinner{{\bf P}\mathopen{\left({\bf
x}\right)}}\equiv\quantity[\mathinner{{\varepsilon}\mathopen{\left({\bf
x}\right)}}-{\varepsilon}_{q}]\mathinner{{\bf E}\mathopen{\left({\bf
x}\right)}}$ (21)

and

$\varepsilon({\bf x})=(\varepsilon_{p}-\varepsilon_{q}){\cal I}^{(p)}({\bf
x})+\varepsilon_{q}$ (22)

is the local dielectric constant, and ${\cal I}^{(p)}({\bf x})$ is the
indicator function for phase $p$ [cf. Eq. (2)]. The vector $\mathinner{{\bf
P}\mathopen{\left({\bf x}\right)}}$ is the induced flux field relative to
reference phase $q$ due to the presence of phase $p$, and hence is zero in the
reference phase $q$ and nonzero in the “polarized” phase $p$ ($p\neq q$).

Using the Green’s function formalism, the local electric field can be
expressed in terms of the following integral equation Rechtsman and Torquato
(2008); Torquato (2002):

$\displaystyle\mathinner{{\bf E}\mathopen{\left({\bf x}\right)}}=$
$\displaystyle\mathinner{{\bf E}_{0}\mathopen{\left({\bf
x}\right)}}+\int\mathinner{\boldsymbol{G}^{(q)}\mathopen{\left({\bf{\bf
x}-{\bf x}^{\prime}}\right)}}\cdot\mathinner{{\bf P}\mathopen{\left({\bf
x}^{\prime}\right)}}\differential{{\bf x}^{\prime}},$ (23)

where the second-rank tensor Green’s function
$\mathinner{\boldsymbol{G}^{(q)}\mathopen{\left({\bf r}\right)}}$ associated
with the reference phase $q$ is given by 222 The Green’s function in this work
(24) differs from the one given in Ref. Rechtsman and Torquato (2008) by a
multiplicative factor $(\omega/c)^{2}$. Thus, Eq. (24) converges to its static
counterpart in the static limit (i.e., $\omega\to 0$).

$\mathinner{\boldsymbol{G}^{(q)}\mathopen{\left({\bf
r}\right)}}=-\boldsymbol{D}^{(q)}\mathinner{\delta\mathopen{\left({\bf
r}\right)}}+\mathinner{\boldsymbol{H}^{(q)}\mathopen{\left({\bf r}\right)}},$
(24)

$\boldsymbol{D}^{(q)}$ is a constant second-rank tensor that arises when one
excludes an infinitesimal region around the position of the singularity ${\bf
x}^{\prime}={\bf x}$ in the Green’s function, and
$\mathinner{\boldsymbol{H}^{(q)}\mathopen{\left({\bf r}\right)}}$ represents
the contribution outside of the “exclusion” region:

$\displaystyle\mathinner{H^{(q)}_{ij}\mathopen{\left({\bf r}\right)}}=$
$\displaystyle i\frac{\pi}{2\varepsilon_{q}}\quantity(\frac{k_{q}}{2\pi
r})^{d/2}\big{\\{}\quantity[k_{q}r\mathinner{\mathcal{H}_{d/2-1}^{(1)}\mathopen{\left(k_{q}r\right)}}-\mathinner{\mathcal{H}_{d/2}^{(1)}\mathopen{\left(k_{q}r\right)}}]\delta_{ij}$
$\displaystyle+k_{q}r\mathinner{\mathcal{H}_{d/2+1}^{(1)}\mathopen{\left(k_{q}r\right)}}\hat{{\bf
r}}_{i}\hat{{\bf r}}_{j}\big{\\}},$ (25)

where $\hat{{\bf r}}\equiv{\bf r}/\absolutevalue{{\bf r}}$ is a unit vector
directed to ${\bf r}$, and
$\mathinner{\mathcal{H}_{\nu}^{(1)}\mathopen{\left(x\right)}}$ is the Hankel
function of the first kind of order $\nu$. The Fourier transform of Eq. (24)
is particularly simple and concise:

$\mathinner{\tilde{G}^{(q)}_{ij}\mathopen{\left({\bf
k}\right)}}=\frac{1}{\varepsilon_{q}}\frac{{k_{q}}^{2}\delta_{ij}-k_{i}k_{j}}{k^{2}-{k_{q}}^{2}}.$
(26)

Note that Eq. (26) is independent of the shape of the exclusion region, which
stands in contrast to the shape-dependent Fourier transform of
$\mathinner{H^{(q)}_{ij}\mathopen{\left({\bf r}\right)}}$; see Supplementary
Material SM  for details.

The use of Eq. (21) and (23) leads to an integral equation for the generalized
cavity intensity field $\mathinner{{\bf F}\mathopen{\left({\bf x}\right)}}$:

$\displaystyle\mathinner{{\bf F}\mathopen{\left({\bf x}\right)}}=$
$\displaystyle\mathinner{{\bf E}_{0}\mathopen{\left({\bf
x}\right)}}+\int_{\epsilon}\mathinner{\boldsymbol{H}^{(q)}\mathopen{\left({\bf{\bf
x}-{\bf x}^{\prime}}\right)}}\cdot\mathinner{{\bf P}\mathopen{\left({\bf
x}^{\prime}\right)}}\differential{{\bf x}^{\prime}},$ (27)

where the integral subscript $\epsilon$ indicates that the integral is to be
carried out by omitting the exclusion region and then allowing it to uniformly
shrink to zero. Here $\mathinner{{\bf F}\mathopen{\left({\bf x}\right)}}$ is
related directly to $\mathinner{{\bf E}\mathopen{\left({\bf x}\right)}}$ via

$\mathinner{{\bf F}\mathopen{\left({\bf
x}\right)}}=\quantity{\boldsymbol{I}+\boldsymbol{D}^{(q)}[\mathinner{\varepsilon\mathopen{\left({\bf
x}\right)}}-\varepsilon_{q}]}\cdot{\bf E}({\bf x}).$ (28)

Using the definitions (21) and (28), we obtain a linear constitutive relation
between $\mathinner{{\bf P}\mathopen{\left({\bf x}\right)}}$ and
$\mathinner{{\bf F}\mathopen{\left({\bf x}\right)}}$:

$\mathinner{{\bf P}\mathopen{\left({\bf
x}\right)}}=\mathinner{\boldsymbol{L}^{(q)}\mathopen{\left({\bf
x}\right)}}\cdot\mathinner{{\bf F}\mathopen{\left({\bf x}\right)}},$ (29)

where

$\mathinner{\boldsymbol{L}^{(q)}\mathopen{\left({\bf
x}\right)}}=[\mathinner{\varepsilon\mathopen{\left({\bf
x}\right)}}-\varepsilon_{q}]\cdot\quantity{\boldsymbol{I}+\boldsymbol{D}^{(q)}[\mathinner{\varepsilon\mathopen{\left({\bf
x}\right)}}-\varepsilon_{q}]}^{-1}.$ (30)

It is noteworthy that one is free to choose any convenient exclusion-region
shape, provided that its boundary is sufficiently smooth. The choice of the
exclusion-region shape is crucially important because it determines the type
of modified electric field that results as well as the corresponding expansion
parameter in the series expansion for the effective dynamic dielectric
constant tensor. For example, when a spheroidal-shaped exclusion region at a
position ${\bf x}$ in $\mathbb{R}^{d}$ is chosen to be aligned with the
polarization vector $\mathinner{{\bf P}\mathopen{\left({\bf x}\right)}}$ at
that position, we have

$\boldsymbol{D}^{(q)}=\frac{A^{*}}{\varepsilon_{q}}\boldsymbol{I},$ (31)

where $\boldsymbol{I}$ is the second-rank identity tensor and $A^{*}\in[0,1]$
is the depolarization factor for a spheroid Torquato (2002) with the
aforementioned alignment. In the special cases of a sphere, disk-like limit,
and needle-like limit, $A^{*}=1/d,1,0$, respectively. Thus, for these three
cases, Eq. (30) yields the following shape-dependent tensor:

$\displaystyle\mathinner{\boldsymbol{L}^{(q)}\mathopen{\left({\bf
x};A^{*}\right)}}=\mathinner{\mathcal{I}^{(p)}\mathopen{\left({\bf
x}\right)}}\boldsymbol{I}$
$\displaystyle\times\begin{cases}d\varepsilon_{q}\beta_{pq},&A^{*}=1/d~{}(\text{spherical})\\\
\varepsilon_{q}\quantity(1-\varepsilon_{q}/\varepsilon_{p}),&A^{*}=1~{}(\text{disk-
like})\\\
\varepsilon_{q}\quantity(\varepsilon_{p}/\varepsilon_{q}-1),&A^{*}=0~{}(\text{needle-
like})\end{cases},$ (32)

where $\beta_{pq}$ is the dielectric polarizability defined by

$\beta_{pq}\equiv\frac{\varepsilon_{p}-\varepsilon_{q}}{\varepsilon_{p}+(d-1)\varepsilon_{q}}.$
(33)

Moreover, in these three cases, the generalized cavity intensity field
$\mathinner{{\bf F}\mathopen{\left({\bf x};A^{*}\right)}}$ reduces to

$\mathinner{{\bf F}\mathopen{\left({\bf
x};A^{*}\right)}}\to\begin{cases}\mathinner{{\bf E}\mathopen{\left({\bf
x}\right)}}+\frac{\displaystyle\mathinner{{\bf P}\mathopen{\left({\bf
x}\right)}}}{\displaystyle
d\varepsilon_{q}},&A^{*}=1/d~{}(\text{spherical})\\\
\frac{\displaystyle\mathinner{{\bf D}\mathopen{\left({\bf
x}\right)}}}{\displaystyle\varepsilon_{q}},&A^{*}=1~{}(\text{disk-like})\\\
\mathinner{{\bf E}\mathopen{\left({\bf x}\right)}},&A^{*}=0~{}(\text{needle-
like}),\end{cases}$ (34)

respectively, where $\mathinner{{\bf D}\mathopen{\left({\bf x}\right)}}$ is
the displacement field. Importantly, among these three cases, the original and
modified strong-contrast expansions arise only when the exclusion region is
taken to be a sphere, as we elaborate below.

We now show that our formalism yields an exact relation between the
polarization field ${\bf P}({\bf x})$ and the applied field ${\bf E}_{0}({\bf
x})$ that is nonlocal in space. It is more convenient at this stage to utilize
a compact linear operator notation, which enables us to express the integral
equation (27) as

$\displaystyle{\bf F}={\bf E}_{0}+\boldsymbol{H}{\bf P},$ (35)

where we temporarily drop the superscript $q$. A combination of this equation
with Eq. (29) yields the following integral equation for the polarization
field

${\bf P}=\boldsymbol{L}{\bf E}_{0}+\boldsymbol{L}\boldsymbol{H}{\bf P}.$ (36)

The desired nonlocal relation is obtained from Eq. (36) by successive
substitutions:

$\displaystyle{\bf P}=\boldsymbol{S}{\bf E}_{0}.$ (37)

where

$\boldsymbol{S}=[\boldsymbol{I}-\boldsymbol{L}\boldsymbol{H}]^{-1}\boldsymbol{L},$
(38)

is a generalized scattering operator that has superior mathematical properties
compared to the scattering operator $\boldsymbol{\mathcal{T}}$ that arises in
standard multiple-scattering theory Frisch (1968); Sheng (1995); Tsang (2001),
as we elaborate below in Remark viii.. More explicitly, the nonlocal relation
(37) can be expressed as

$\displaystyle\mathinner{{\bf P}\mathopen{\left({\bf 1}\right)}}=$
$\displaystyle\int_{\epsilon}\differential{{\bf
2}}\mathinner{\boldsymbol{S}\mathopen{\left({\bf 1},{\bf
2}\right)}}\cdot\mathinner{{\bf E}_{0}\mathopen{\left({\bf 2}\right)}},$ (39)

where boldface numbers ${\bf 1},{\bf 2}$ are shorthand notations for position
vectors ${\bf r}_{1},{\bf r}_{2}$. Ensemble averaging Eq. (39) and invoking
statistical homogeneity yields the convolution relation

$\mathinner{\left\langle{\bf P}\right\rangle\mathopen{\left({\bf
1}\right)}}=\int_{\epsilon}\differential{{\bf
2}}\mathinner{\left\langle\boldsymbol{S}\right\rangle\mathopen{\left({\bf
1}-{\bf 2}\right)}}\cdot\mathinner{{\bf E}_{0}\mathopen{\left({\bf
2}\right)}},$ (40)

where the operator $\left\langle\boldsymbol{S}\right\rangle$ depends on
relative positions, i.e.,
$\mathinner{\left\langle\boldsymbol{S}\right\rangle\mathopen{\left({\bf
1},{\bf
2}\right)}}=\mathinner{\left\langle\boldsymbol{S}\right\rangle\mathopen{\left({\bf
1}-{\bf 2}\right)}}$, and angular brackets denote an ensemble average.
Formally, the nonlocal relation (40) is the same as the one given in Torquato
Torquato (2002) for the static problem, but nonlocality was not explicitly
invoked there. Taking the Fourier transform of (40) yields a compact Fourier
representation of this nonlocal relation, namely,

$\mathinner{\widetilde{\left\langle{\bf P}\right\rangle}\mathopen{\left({\bf
k}\right)}}=\mathinner{\widetilde{\left\langle\boldsymbol{S}\right\rangle}\mathopen{\left({\bf
k}\right)}}\cdot\mathinner{\tilde{{\bf E}}_{0}\mathopen{\left({\bf
k}\right)}},$ (41)

where $\mathinner{\widetilde{\left\langle f\right\rangle}\mathopen{\left({\bf
k}\right)}}\equiv\int\mathinner{\left\langle
f\right\rangle\mathopen{\left({\bf
x}\right)}}\mathinner{\exp\mathopen{\left(-i{\bf k}\cdot{\bf
x}\right)}}\differential{{\bf x}}$. From Eq. (19), $\mathinner{\tilde{{\bf
E}}_{0}\mathopen{\left({\bf k}\right)}}=\tilde{{\bf
E}}_{0}\mathinner{\delta\mathopen{\left({\bf k}-{\bf k}_{q}\right)}}$,
implying that the wave vector ${\bf k}$ in Eq. (41) must be identical to ${\bf
k}_{q}$.

As in the static case Sen and Torquato (1989); Torquato (2002) and quasistatic
regime Rechtsman and Torquato (2008), the ensemble-averaged operator
$\mathinner{\left\langle\boldsymbol{S}\right\rangle\mathopen{\left({\bf
r}\right)}}$, which is given explicitly in the Supplementary Material SM  in
terms of the $n$-point correlation functions and products of the tensor
$\boldsymbol{H}({\bf r})$, depends on the shape of the macroscopic ellipsoidal
composite specimen (see Fig. 2). This shape-dependence arises from because
$\boldsymbol{H}$ decays like $r^{-d}$ for large $r$, and hence,
$\mathinner{\left\langle\boldsymbol{S}\right\rangle\mathopen{\left({\bf
r}\right)}}$ involves conditionally convergent integrals Torquato (2002). To
avoid such conditional convergence issues, we follow previous strong-contrast
formulations by seeking to eliminate the applied field ${\bf E}_{0}$ in Eq.
(40) in favor of the average cavity field $\mathinner{\left\langle{\bf
F}\right\rangle\mathopen{\left({\bf r}\right)}}$ in order to get a
corresponding nonlocal homogenized constitutive relation between
$\mathinner{\left\langle{\bf P}\right\rangle\mathopen{\left({\bf r}\right)}}$
and $\mathinner{\left\langle{\bf F}\right\rangle\mathopen{\left({\bf
r}\right)}}$ or vice versa. Thus, solving for ${\bf E}_{0}$ in Eq. (40) and
substituting into the ensemble average of (35) yields

$\langle{\bf
F}\rangle=[\left\langle\boldsymbol{S}\right\rangle^{-1}+\boldsymbol{H}]\left\langle{\bf
P}\right\rangle.$ (42)

Inverting this expression leads to the following nonlocal constitutive
relation:

$\mathinner{\left\langle{\bf P}\right\rangle\mathopen{\left({\bf
1}\right)}}=\int\differential{{\bf 2}}\boldsymbol{L}_{e}^{(q)}({\bf 1}-{\bf
2})\cdot\mathinner{\left\langle{\bf F}\right\rangle\mathopen{\left({\bf
2}\right)}},$ (43)

where $\boldsymbol{L}_{e}^{(q)}({\bf r})$ is a kernel that is derived
immediately below and explicitly given by

$\displaystyle\mathinner{\boldsymbol{L}_{e}^{(q)}\mathopen{\left({\bf
r}\right)}}\equiv$ $\displaystyle\int\differential{{\bf
r}^{\prime}}\quantity[\mathinner{\boldsymbol{\varepsilon}_{e}\mathopen{\left({\bf
r}^{\prime}\right)}}-\varepsilon_{q}\boldsymbol{I}\mathinner{\delta\mathopen{\left({\bf
r}^{\prime}\right)}}]\cdot\Big{[}\boldsymbol{I}\mathinner{\delta\mathopen{\left({\bf
r}-{\bf r}^{\prime}\right)}}$
$\displaystyle+\boldsymbol{D}^{(q)}\cdot\quantity(\mathinner{\boldsymbol{\varepsilon}_{e}\mathopen{\left({\bf
r}-{\bf
r}^{\prime}\right)}}-\varepsilon_{q}\boldsymbol{I}\mathinner{\delta\mathopen{\left({\bf
r}-{\bf r}^{\prime}\right)}})\Big{]}^{-1}.$ (44)

We are not aware of any previous work that derives such an exact nonlocal
homogenized constitutive relation (43) from first principles.

Note that the support $\ell_{s}$ of the kernel $\boldsymbol{L}_{e}^{(q)}({\bf
r})$ relative to the incident wavelength $\lambda$ determines the degree of
spatial dispersion. When $\lambda$ is finite, the relation between
$\mathinner{\left\langle{\bf P}\right\rangle\mathopen{\left({\bf x}\right)}}$
and $\mathinner{\left\langle{\bf F}\right\rangle\mathopen{\left({\bf
x}\right)}}$ in Eq. (43) is nonlocal in space. In the regime
$\ell_{s}\ll\lambda$, the nonlocal relation (43) can be well approximated by
the local relation $\mathinner{\left\langle{\bf
P}\right\rangle\mathopen{\left({\bf
x}\right)}}\approx\quantity[\int\mathinner{\boldsymbol{L}_{e}^{(q)}\mathopen{\left({\bf
x}^{\prime}\right)}}\differential{{\bf
x}^{\prime}}]\cdot\mathinner{\left\langle{\bf
F}\right\rangle\mathopen{\left({\bf x}\right)}}$. Indeed, in the static limit,
$\boldsymbol{L}_{e}^{(q)}({\bf r})$ tends to a Dirac delta function
$\mathinner{\delta\mathopen{\left({\bf r}\right)}}$, expression (43) becomes
the position-independent local relation

$\left\langle{\bf
P}\right\rangle=\boldsymbol{L}_{e}^{(q)}\cdot\left\langle{\bf F}\right\rangle$
(45)

derived earlier Torquato (2002).

The nonlocal constitutive relation in direct space, Eq. (43), can be reduced
to a linear product form in Fourier space by taking the Fourier transform of
Eq. (43):

$\displaystyle\mathinner{\widetilde{\left\langle{\bf
P}\right\rangle}\mathopen{\left({\bf
k}_{q}\right)}}=\mathinner{\boldsymbol{L}_{e}^{(q)}\mathopen{\left({\bf
k}_{q}\right)}}\cdot\mathinner{\widetilde{\left\langle{\bf
F}\right\rangle}\mathopen{\left({\bf k}_{q}\right)}}.$ (46)

The wave-vector-dependent effective tensor
$\mathinner{\boldsymbol{L}_{e}^{(q)}\mathopen{\left({\bf k}_{q}\right)}}$ is
postulated (see discussion in the Supplementary Material SM ) to be given by

$\mathinner{\boldsymbol{L}_{e}^{(q)}\mathopen{\left({\bf
k}_{q}\right)}}\equiv[\mathinner{\boldsymbol{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}-\varepsilon_{q}\boldsymbol{I}]\cdot\quantity{\boldsymbol{I}+\boldsymbol{D}^{(q)}\cdot[\mathinner{\boldsymbol{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}-\varepsilon_{q}\boldsymbol{I}]}^{-1}.$ (47)

This linear fractional form for
$\mathinner{\boldsymbol{L}_{e}^{(q)}\mathopen{\left({\bf k}\right)}}$ is
consistent with the one derived for the static limit Torquato (2002) and for
the quasistatic regime Rechtsman and Torquato (2008). Taking the inverse
Fourier transform of Eq. (47) yields its corresponding direct-space
representation Eq. (44). Taking the Fourier transform of Eq. (42) yields

$\displaystyle\mathinner{\widetilde{\left\langle{\bf
F}\right\rangle}\mathopen{\left({\bf k}_{q}\right)}}=$
$\displaystyle\quantity[{\mathinner{\widetilde{\left\langle\boldsymbol{S}\right\rangle}\mathopen{\left({\bf
k}_{q}\right)}}}^{-1}+\mathinner{\boldsymbol{\tilde{H}}\mathopen{\left({\bf
k}_{q}\right)}}]\cdot\mathinner{\widetilde{\left\langle{\bf
P}\right\rangle}\mathopen{\left({\bf k}_{q}\right)}}.$ (48)

Comparing Eq. (46) to Eq. (48), and specifically choosing a spherical
exclusion region, as discussed in Eq. (32), yields the desired exact strong-
contrast expansions for general macroscopically anisotropic two-phase media:

$\displaystyle{\phi_{p}}^{2}\beta^{2}_{pq}[\boldsymbol{\varepsilon}_{e}({\bf
k}_{q})+(d-1)\varepsilon_{q}\boldsymbol{I}]\cdot[\boldsymbol{\varepsilon}_{e}({\bf
k}_{q})-\varepsilon_{q}\boldsymbol{I}]^{-1}=\phi_{p}\beta_{pq}\boldsymbol{I}-\sum_{n=2}^{\infty}\mathinner{\boldsymbol{A}_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{n},$ (49)

where $\mathinner{\boldsymbol{A}_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}$ is a wave-vector-dependent second-rank tensor that is a
functional involving the set of correlation functions
$S_{1}^{(p)},S_{2}^{(p)},\cdots,S_{n}^{(p)}$ (defined in Sec. II.1) and
products of the second-rank tensor field
$\mathinner{\boldsymbol{H}^{(q)}\mathopen{\left({\bf r}\right)}}$, which is
explicitly given as [see also Eq. (25)]

$\displaystyle\mathinner{H^{(q)}_{ij}\mathopen{\left({\bf r}\right)}}=$
$\displaystyle\begin{cases}\frac{i}{4}\quantity[{k_{q}}^{2}\mathinner{\mathcal{H}_{0}^{(1)}\mathopen{\left(k_{q}r\right)}}-\frac{k_{q}}{r}\mathinner{\mathcal{H}_{1}^{(1)}\mathopen{\left(k_{q}r\right)}}]\delta_{ij}+\frac{i{k_{q}}^{2}}{4}\mathinner{\mathcal{H}_{2}^{(1)}\mathopen{\left(k_{q}r\right)}}\hat{{\bf
r}}_{i}\hat{{\bf r}}_{j},&d=2\\\
\frac{\mathinner{\exp\mathopen{\left(ik_{q}r\right)}}}{\varepsilon_{q}4\pi
r^{3}}\big{\\{}\quantity[-1+ik_{q}r+(k_{q}r)^{2}]\delta_{ij}+\quantity[3-3ik_{q}r-(k_{q}r)^{2}]\hat{{\bf
r}}_{i}\hat{{\bf r}}_{j}\big{\\}},&d=3\end{cases}$ (50)

where $k_{q}\equiv\absolutevalue{{\bf k}_{q}}$. Specifically, for $n=2$ and
$n\geq 3$, these $n$-point tensors associated with the polarized phase $p$
are, respectively, given by

$\displaystyle\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}=$ $\displaystyle
d\varepsilon_{q}\int_{\epsilon}\differential{{\bf
r}}\mathinner{\boldsymbol{H}^{(q)}\mathopen{\left({\bf r}\right)}}e^{-i{\bf
k}_{q}\cdot{\bf r}}\mathinner{\chi_{{}_{V}}\mathopen{\left({\bf r}\right)}},$
(51) $\displaystyle\mathinner{\boldsymbol{A}_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}=$ $\displaystyle
d\varepsilon_{q}\quantity(\frac{-d\varepsilon_{q}}{\phi_{p}})^{n-2}\int_{\epsilon}\differential{{\bf
r}_{1}}\cdots\differential{{\bf
r}_{n-1}}\mathinner{\boldsymbol{H}^{(q)}\mathopen{\left({\bf r}_{1}-{\bf
r}_{2}\right)}}e^{-i{\bf k}_{q}\cdot({\bf r}_{1}-{\bf
r}_{2})}\cdot\mathinner{\boldsymbol{H}^{(q)}\mathopen{\left({\bf r}_{2}-{\bf
r}_{3}\right)}}e^{-i{\bf k}_{q}\cdot({\bf r}_{2}-{\bf r}_{3})}\cdot$
$\displaystyle~{}\cdots~{}\mathinner{\boldsymbol{H}^{(q)}\mathopen{\left({\bf
r}_{n-1}-{\bf r}_{n}\right)}}e^{-i{\bf k}_{q}\cdot({\bf r}_{n-1}-{\bf
r}_{n})}\mathinner{\Delta_{n}^{(p)}\mathopen{\left({\bf r}_{1},\cdots,{\bf
r}_{n}\right)}},$ (52)

where $\int_{\epsilon}\differential{{\bf r}}\equiv\lim_{\epsilon\to
0^{+}}\int_{\absolutevalue{{\bf r}}>\epsilon}\differential{{\bf r}}$ and
$\Delta_{n}^{(p)}$ is a position-dependent determinant involving correlation
functions of the polarized phase $p$ up to the $n$-point level:

$\mathinner{\Delta_{n}^{(p)}\mathopen{\left({\bf r}_{1},\cdots,{\bf
r}_{n}\right)}}=\begin{vmatrix}\mathinner{S_{2}^{(p)}\mathopen{\left({\bf
r}_{1},{\bf r}_{2}\right)}}&\mathinner{S_{1}^{(p)}\mathopen{\left({\bf
r}_{1}\right)}}&\cdots&0\\\ \mathinner{S_{3}^{(p)}\mathopen{\left({\bf
r}_{1},{\bf r}_{2},{\bf
r}_{3}\right)}}&\mathinner{S_{2}^{(p)}\mathopen{\left({\bf r}_{2},{\bf
r}_{3}\right)}}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
\mathinner{S_{n}^{(p)}\mathopen{\left({\bf r}_{1},\cdots,{\bf
r}_{n}\right)}}&\mathinner{S_{n-1}^{(p)}\mathopen{\left({\bf
r}_{2},\cdots,{\bf
r}_{n}\right)}}&\cdots&\mathinner{S_{2}^{(p)}\mathopen{\left({\bf
r}_{n-1},{\bf r}_{n}\right)}}\end{vmatrix}.$ (53)

For macroscopically isotropic media, the effective dielectric tensor is
isotropic, i.e., $\mathinner{\boldsymbol{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}=\mathinner{\varepsilon_{e}\mathopen{\left({\bf
k}_{q}\right)}}\boldsymbol{I}$. The corresponding strong-contrast expansion
for $\mathinner{\varepsilon_{e}\mathopen{\left({\bf k}_{q}\right)}}$ is
obtained by taking the trace of both sides of Eq. (49):

$\displaystyle{\phi_{p}}^{2}\beta^{2}_{pq}[\varepsilon_{e}({\bf
k}_{q})+(d-1)\varepsilon_{q}][\varepsilon_{e}({\bf
k}_{q})-\varepsilon_{q}]^{-1}$ $\displaystyle=$
$\displaystyle\phi_{p}\beta_{pq}-\sum_{n=2}^{\infty}\mathinner{A_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{n},$ (54)

where $\mathinner{\varepsilon_{e}\mathopen{\left({\bf
k}_{q}\right)}}=\Tr[\mathinner{\boldsymbol{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}]/d$, $\mathinner{A_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}=\Tr[\mathinner{\boldsymbol{A}_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}]/d$ for $n\geq 2$ and $\Tr[\;]$ denotes the trace operation.
Furthermore, for statistically isotropic media, the effective dielectric
constant becomes independent of the direction of the wave vector, i.e.,
$\mathinner{\varepsilon_{e}\mathopen{\left({\bf
k}_{q}\right)}}=\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}$.

Remarks:

  1. i.

Importantly, the strong-contrast expansion (49) is a series representation of
a linear fractional transformation of the variable
$\boldsymbol{\varepsilon}_{e}({\bf k}_{q})$ (left-hand side of the equation),
rather than the effective dielectric constant tensor itself. The series
expansion in powers of the polarizability $\beta_{pq}$ of this particular
rational function of $\boldsymbol{\varepsilon}_{e}({\bf k}_{q})$ has important
consequences for the predictive power of approximations derived from the
expansion, as detailed in Sec. III.2.

  2. ii.

The fact that the exact expansion (49), extracted from our nonlocal relation
(46), is explicitly given in terms of integrals over products of the relevant
Green’s functions and the $n$-point correlation functions to infinite order
implies that multiple-scattering to all orders is exactly treated for the
range of wave numbers for which our extended homogenization theory applies,
i.e., $0\leq|{\bf k}_{q}|\ell\lesssim 1$.

  3. iii.

Note that Eq. (49) represents two different series expansions: one for $q=1$
and $p=2$ and the other for $q=2$ and $p=1$.

  4. iv.

The exact expansions represented by (49) are independent of the reference
phase $q$ and hence independent of the wave vector ${\bf k}_{q}$.

  5. v.

For $d=2$, the strong-contrast expansion applies for TE polarization only.
This implies that the electric field and wave vector are parallel to the plane
or transverse to an axis of symmetry in a 3D system whose cross sections are
identical.

  6. vi.

Formally, the original strong-contrast expansions that apply in the
quasistatic regime Rechtsman and Torquato (2008) can be obtained from the
nonlocal strong-contrast expansions (49) by simply replacing the exponential
functions that appear in the expressions for the second-rank tensors
$\mathinner{\boldsymbol{A}_{n}^{(p)}\mathopen{\left({\bf k}_{q}\right)}}$,
defined by Eq. (52), by unity.

  7. vii.

For statistically isotropic media, the effective phase speed
$\mathinner{c_{e}\mathopen{\left(k_{q}\right)}}$ and attenuation coefficient
$\mathinner{\gamma_{e}\mathopen{\left(k_{q}\right)}}$ are determined by the
scalar effective dielectric constant
$\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}$:

$\displaystyle\mathinner{c_{e}\mathopen{\left(k_{q}\right)}}/c=$
$\displaystyle{n_{e}(k_{q})}^{-1}=\real\quantity[1/\sqrt{\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}}],$
(55) $\displaystyle\mathinner{\gamma_{e}\mathopen{\left(k_{q}\right)}}/c=$
$\displaystyle{\kappa_{e}(k_{q})}^{-1}=\imaginary\quantity[1/\sqrt{\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}}],$
(56)

where ${n_{e}}(k_{q})$ and
$\mathinner{\kappa_{e}\mathopen{\left(k_{q}\right)}}$ are the effective
refractive and extinction indices, respectively. The quantity
$\mathinner{\exp\mathopen{\left(-2\pi\gamma_{e}/c_{e}\right)}}$ is the factor
by which the incident-wave amplitude is attenuated for a period $2\pi/\omega$.

  8. viii.

The strong-contrast formalism is a significant departure from perturbative
expansions obtained from standard multiple-scattering theory Frisch (1968);
Sheng (1995); Tsang (2001); Cazé and Schotland (2015). The operator
$\boldsymbol{S}$ defined in Eq. (38) is a generalization of the standard
scattering operator
$\boldsymbol{\mathcal{T}}=[\boldsymbol{I}-\boldsymbol{V}\boldsymbol{G}]^{-1}\boldsymbol{V}$
Frisch (1968); Sheng (1995); Tsang (2001), where
$\boldsymbol{V}\equiv\quantity[\mathinner{\varepsilon\mathopen{\left({\bf
r}\right)}}-\varepsilon_{q}]\boldsymbol{I}$ is the scattering potential. Thus,
the perturbation series resulting from $\boldsymbol{\mathcal{T}}$ is a weak-
contrast expansion with the inherent limitations that it converges only for
small contrast ratios (see also Sec. III.2). By contrast, due to the different
possible choices for the exclusion regions and reference phases, there is an
infinite variety of series expansions that result from the strong-contrast
formalism with generally fast-convergence properties. The particular strong-
contrast expansion can be designed for different classes of microstructures
(see Appendixes A and D). The corresponding strong-contrast “self-energy”
$\boldsymbol{\Sigma}$ is a linear fractional transform of $\boldsymbol{L}_{e}$
[see Eq. (48)], namely,
$\boldsymbol{\Sigma}={\boldsymbol{L}_{e}}^{-1}\quantity[\boldsymbol{I}-\boldsymbol{D}{\boldsymbol{L}_{e}}^{-1}]^{-1}$,
which is a generalization of the self-energy in standard multiple-scattering
theory Frisch (1968); Sheng (1995); Tsang (2001); Cazé and Schotland (2015).
Thus, it is highly nontrivial to relate diagrammatic expansions of the strong-
contrast formalism to those of multiple-scattering theory. An elaboration of
how strong-contrast expansion generalize those from multiple-scattering theory
is presented in the Supplementary Material SM .

###  III.2 Convergence properties and accuracy of truncated series

The form of the strong-contrast expansion parameter $\beta_{pq}$ in Eq. (49)
is a direct consequence of the choice of a spherical region excluded from the
volume integrals in Eq. (49) due to singularities in the Green’s functions
Rechtsman and Torquato (2008). It is bounded by

$-\frac{1}{d-1}\leq\beta_{pq}\equiv\frac{\varepsilon_{p}-\varepsilon_{q}}{\varepsilon_{p}+(d-1)\varepsilon_{q}}\leq
1,$ (57)

which implies that the strong-contrast expansion (49) can converge rapidly,
even for infinite contrast ratio $\varepsilon_{p}/\varepsilon_{q}\to\infty$.
Other choices for the shape of the exclusion region will lead to different
expansion parameters that will generally be bounded but can lead to expansions
with significantly different convergence properties Torquato (2002). In
Appendix A, we present the corresponding expansions for disk-like and needle-
like exclusion regions, which are exceptional cases that lead to slowly
converging weak-contrast expansions with expansion parameter
$(\varepsilon_{p}-\varepsilon_{q})/\varepsilon_{q}$, and thus are unbounded
when $\varepsilon_{p}/\varepsilon_{q}\to\infty$.

Figure 3: Schematic of the optimal multiscale “coated-spheres” model that
realizes the isotropic Hashin-Shtrikman bounds on $\varepsilon_{e}$ Hashin and
Shtrikman (1963). Each composite sphere is composed of a spherical inclusion
of one phase (dispersed phase) that is surrounded by a concentric spherical
shell of the other phase such that the fraction of space occupied by the
dispersed phase is equal to its overall phase volume fraction. The composite
spheres fill all space, implying that their sizes range down to the
infinitesimally small. When phase 2 is the disconnected inclusion (dispersed)
phase, this two-phase medium minimizes and maximizes the effective static
dielectric constant $\varepsilon_{e}$ for prescribed volume fraction and
contrast ratio, when $\varepsilon_{2}/\varepsilon_{1}>1$ and
$\varepsilon_{2}/\varepsilon_{1}<1$, respectively. It has recently been proved
that these highly degenerate optimal Hashin-Shtrikman multiscale distributions
of spheres are hyperuniform Kim and Torquato (2019a, b).

Importantly, in the purely static case, the expansion (49) becomes identical
to one derived by Sen and Torquato Sen and Torquato (1989) and its truncation
after second-order terms [i.e., setting $\boldsymbol{A}_{n}^{(p)}=0$ for all
$n\geq 3$] yields the generalized Hashin-Shtrikman bounds Hashin and Shtrikman
(1963) derived by Willis Willis (1981) that are optimal since they are
realized by certain statistically anisotropic composites in which there is a
disconnected, dispersed phase in a connected matrix phase Milton (1981). In
the case of an isotropic effective dielectric constant $\varepsilon_{e}$, the
optimal Hashin-Shtrikman upper and lower bounds for any phase-contrast ratio
$\varepsilon_{2}/\varepsilon_{1}$ are exactly realized by the multiscale
“coated-spheres” model, which is depicted in Fig. 3 in two dimensions. Affine
transformations of the coated spheres in the $d$ orthogonal directions lead to
oriented coated ellipsoids that are optimal for the macroscopically
anisotropic case. The lower bound corresponds to the case when the high-
dielectric-constant phase is the dispersed, disconnected phase and the upper
bound corresponds to the instance in which the high-dielectric-constant phase
is the connected matrix. Thus, Torquato Torquato (1985, 2002) observed that
the strong-contrast expansions (49) in the static limit can be regarded as
ones that perturb around such optimal composites, implying that the first few
terms of the expansion can yield accurate approximations of the effective
property for a class of particulate composites as well as more general
microstructures, depending on whether the high-dielectric phase percolates or
not. For example, even when $\varepsilon_{2}/\varepsilon_{1}\gg 1$, the
dispersed phase 2 can consist of identical or polydisperse particles of
general shape (ellipsoids, cubes, cylinders, polyhedra) with prescribed
orientations that may or not overlap, provided that the particles are
prevented from forming large clusters compared to the specimen size. Moreover,
when $\varepsilon_{2}/\varepsilon_{1}\ll 1$, the matrix phase can be a
cellular network Torquato and Chen (2018). Finally, for moderate values of the
contrast ratio $\varepsilon_{2}/\varepsilon_{1}$, even more general
microstructures (e.g., those without well-defined inclusions) can be
accurately treated. Importantly, we show that for the dynamic problem under
consideration, the first few terms of the expansion (49) yield accurate
approximations of $\boldsymbol{\varepsilon}_{e}({\bf k}_{q})$ for a similar
wide class of two-phase media (see Sec. VIII). Analogous approximations were
derived and applied for the quasistatic regime Rechtsman and Torquato (2008);
Chen and Torquato (2018).

We now show how lower-order truncations of the series (49) can well
approximate higher-order functionals (i.e., higher-order diagrams) of the
exact series to all orders in terms of lower-order diagrams. Such truncations
of strong-contrast expansions are tantamount to approximate but resummations
of the strong expansions, which enables multiple-scattering and spatial
dispersion effects to be accurately captured to all orders. Solving the left-
hand side of Eq. (49) for $\boldsymbol{\varepsilon}_{e}$ yields the rational
function in $\beta_{pq}$:

$\displaystyle\frac{\mathinner{\boldsymbol{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}}{\varepsilon_{q}}=$
$\displaystyle\boldsymbol{I}+d{\phi_{p}}^{2}\beta_{pq}\Big{[}\phi_{p}(1-\phi_{p}\beta_{pq})\boldsymbol{I}$
$\displaystyle-\sum_{n=2}^{\infty}\mathinner{\boldsymbol{A}_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{n-1}\Big{]}^{-1}.$ (58)

Expanding Eq. (58) in powers of the scalar polarizability $\beta_{pq}$ yields
the series

$\frac{\mathinner{\boldsymbol{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}}{\varepsilon_{q}}=\sum_{n=0}^{\infty}\mathinner{\boldsymbol{B}_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{n},$ (59)

where the first several functionals
$\mathinner{\boldsymbol{B}_{n}^{(p)}\mathopen{\left({\bf k}_{q}\right)}}$ are
explicitly given in terms of
$\boldsymbol{A}_{0}^{(p)},~{}\boldsymbol{A}_{1}^{(p)},\ldots,~{}\boldsymbol{A}_{n}^{(p)}$
as

$\displaystyle\mathinner{\boldsymbol{B}_{0}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}=$ $\displaystyle\boldsymbol{I}$
$\displaystyle\mathinner{\boldsymbol{B}_{1}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}=$ $\displaystyle d\phi_{p}\boldsymbol{I}$
$\displaystyle\mathinner{\boldsymbol{B}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}=$ $\displaystyle
d\quantity[\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}+\phi_{p}^{2}\boldsymbol{I}]$
$\displaystyle\mathinner{\boldsymbol{B}_{3}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}=$
$\displaystyle\frac{d}{\phi_{p}}\Big{[}{\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}}^{2}+\phi_{p}\mathinner{\boldsymbol{A}_{3}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}$
$\displaystyle\quad+2\phi_{p}^{2}\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}+\phi_{p}^{4}\boldsymbol{I}\Big{]}$
$\displaystyle\mathinner{\boldsymbol{B}_{4}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}=$
$\displaystyle\frac{d}{\phi_{p}^{2}}\Bigg{[}\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}^{3}+2\phi_{p}\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}\cdot\mathinner{\boldsymbol{A}_{3}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}$
$\displaystyle~{}+3\phi_{p}^{2}{\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}}^{2}+\phi_{p}^{2}\mathinner{\boldsymbol{A}_{4}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}$
$\displaystyle~{}+2\phi_{p}^{3}\mathinner{\boldsymbol{A}_{3}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}+3\phi_{p}^{4}\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}+\phi_{p}^{6}\boldsymbol{I}\Bigg{]},$

where $\boldsymbol{T}^{n}$ stands for $n$ successive inner products of a
second-rank tensor $\boldsymbol{T}$.

Let us now compare the exact expansion (59) to the one that results when
expanding the truncation of the exact expression (49) for
$[\boldsymbol{\varepsilon}_{e}+(d-1)\boldsymbol{\varepsilon}_{q}]\cdot(\boldsymbol{\varepsilon}_{e}-\boldsymbol{\varepsilon}_{q})^{-1}$
at the two-point level:

$\displaystyle\frac{\mathinner{\boldsymbol{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}}{\varepsilon_{q}}\approx$
$\displaystyle\boldsymbol{I}+d{\phi_{p}}\beta_{pq}\quantity[{(1-\phi_{p}\beta_{pq})\boldsymbol{I}-\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}\beta_{pq}/\phi_{p}}]^{-1}$ (60) $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\mathinner{\boldsymbol{C}^{(p)}_{n}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{n},$ (61)

where the $n$th-order functional
$\mathinner{\boldsymbol{C}^{(p)}_{n}\mathopen{\left({\bf k}_{q}\right)}}$ for
any $n$ is given in terms of the volume fraction $\phi_{p}$ and
$\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf k}_{q}\right)}}$,
which has the following diagrammatic representation:

$\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}=\begin{minipage}[h]{72.26999pt}
\includegraphics[width=72.26999pt]{Diagram.eps} \end{minipage}.$ (62)

Here the solid and wavy lines joining two nodes represent the spatial
correlation via $\mathinner{\chi_{{}_{V}}\mathopen{\left({\bf r}\right)}}$ and
a wave vector-dependent Green’s function
$\mathinner{\boldsymbol{H}^{(q)}\mathopen{\left({\bf r}\right)}}e^{-i{\bf
k}_{q}\cdot{\bf r}}$ between the nodes, respectively. The black node indicates
a volume integral and carries a factor of $d\varepsilon_{q}$. The first
several functionals $\mathinner{\boldsymbol{C}^{(p)}_{n}\mathopen{\left({\bf
k}_{q}\right)}}$ are explicitly given as

$\displaystyle\mathinner{\boldsymbol{C}^{(p)}_{0}\mathopen{\left({\bf
k}_{q}\right)}}$ $\displaystyle=\boldsymbol{I},$
$\displaystyle\mathinner{\boldsymbol{C}^{(p)}_{1}\mathopen{\left({\bf
k}_{q}\right)}}$ $\displaystyle=d\phi_{p}\boldsymbol{I},$
$\displaystyle\mathinner{\boldsymbol{C}^{(p)}_{2}\mathopen{\left({\bf
k}_{q}\right)}}$
$\displaystyle=d\quantity[\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}+{\phi_{p}}^{2}\boldsymbol{I}],$
$\displaystyle\mathinner{\boldsymbol{C}^{(p)}_{3}\mathopen{\left({\bf
k}_{q}\right)}}$
$\displaystyle=\frac{d}{\phi_{p}}\quantity[\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}+{\phi_{p}}^{2}\boldsymbol{I}]^{2},$
$\displaystyle\mathinner{\boldsymbol{C}^{(p)}_{4}\mathopen{\left({\bf
k}_{q}\right)}}$
$\displaystyle=\frac{d}{{\phi_{p}}^{2}}\quantity[\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}+{\phi_{p}}^{2}\boldsymbol{I}]^{3}.$

Thus, comparing (59) to (61), we see that truncation of the expansion (49) for
$[\boldsymbol{\varepsilon}_{e}+(d-1)\boldsymbol{\varepsilon}_{q}]\cdot(\boldsymbol{\varepsilon}_{e}-\boldsymbol{\varepsilon}_{q})^{-1}$
at the two-point level actually translates into approximations of the higher-
order functionals to all orders in terms of the first-order diagram $\phi_{p}$
and the second-order diagram (62). This two-point truncation can be thought of
as an approximate but accurate resummed representation of the exact expansion
(59). Note that the approximate expansion (61) is exact through second order
in $\beta_{pq}$. Clearly, truncation of Eq. (49) at the three-point level (see
Appendix C) will yield even better approximations of the higher-order
functionals.

##  IV Strong-Contrast Approximation Formulas

Here we describe lower-order truncations of the strong-contrast expansions
that are expected to yield accurate closed-form formulas for
$\boldsymbol{\varepsilon}_{e}({\bf k}_{q})$ that apply over a broad range of
wavelengths ($k_{q}\ell\lesssim 1$), volume fractions and contrast ratios for
a wide class of microstructures.

###  IV.1 Macroscopically anisotropic media

For the ensuing treatment, it is convenient to rewrite the expansions (49),
valid for macroscopically anisotropic media in $\mathbb{R}^{d}$, in the
following manner:

$\displaystyle{\phi_{p}}^{2}\beta^{2}_{pq}\frac{\boldsymbol{\varepsilon}_{e}({\bf
k}_{q})+(d-1)\varepsilon_{q}\boldsymbol{I}}{\boldsymbol{\varepsilon}_{e}({\bf
k}_{q})-\varepsilon_{q}\boldsymbol{I}}$ $\displaystyle=$
$\displaystyle\phi_{p}\beta_{pq}\boldsymbol{I}-\sum_{n=2}^{M}\mathinner{\boldsymbol{A}_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{n}+\mathinner{\mathcal{R}_{M}\mathopen{\left({\bf
k}_{q}\right)}},$ (63)

where the $M$th-order remainder term is defined as

$\mathinner{\mathcal{R}_{M}\mathopen{\left({\bf
k}_{q}\right)}}\equiv\sum_{n=M+1}^{\infty}\mathinner{\boldsymbol{A}_{n}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{n}.$ (64)

Truncating the exact nonlocal expansion (63) at the two- and three-point
levels, i.e., setting $\mathinner{\mathcal{R}_{2}\mathopen{\left({\bf
k}_{q}\right)}}=0$ and $\mathinner{\mathcal{R}_{3}\mathopen{\left({\bf
k}_{q}\right)}}=0$, respectively, yields

$\displaystyle{\phi_{p}}^{2}\beta^{2}_{pq}\frac{\boldsymbol{\varepsilon}_{e}({\bf
k}_{q})+(d-1)\varepsilon_{q}\boldsymbol{I}}{\boldsymbol{\varepsilon}_{e}({\bf
k}_{q})-\varepsilon_{q}\boldsymbol{I}}=$
$\displaystyle\phi_{p}\beta_{pq}\boldsymbol{I}-\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{2},$ (65)
$\displaystyle{\phi_{p}}^{2}\beta^{2}_{pq}\frac{\boldsymbol{\varepsilon}_{e}({\bf
k}_{q})+(d-1)\varepsilon_{q}\boldsymbol{I}}{\boldsymbol{\varepsilon}_{e}({\bf
k}_{q})-\varepsilon_{q}\boldsymbol{I}}=$
$\displaystyle\phi_{p}\beta_{pq}\boldsymbol{I}-\Big{[}\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{2}$
$\displaystyle\quad+\mathinner{\boldsymbol{A}_{3}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}{\beta_{pq}}^{3}\Big{]},$ (66)

Compared to the quasistatic approximation Rechtsman and Torquato (2008), these
nonlocal approximations substantially extend the range of applicable wave
number, namely, $0\leq\absolutevalue{{\bf k}_{q}}\ell\lesssim 1$.

###  IV.2 Macroscopically isotropic media

All of the applications considered in this paper, will focus on the case of
macroscopically isotropic media, i.e., they are described by the scalar
effective dielectric constant $\mathinner{\varepsilon_{e}\mathopen{\left({\bf
k}_{q}\right)}}=\mbox{Tr}\,[\boldsymbol{\varepsilon}_{e}({\bf k}_{q})]/d$ but
depend on the direction of the wave vector ${\bf k}_{q}$.

####  IV.2.1 Strong-contrast approximation at the two-point level

Solving Eq. (54) for the effective dielectric constant
$\mathinner{\varepsilon_{e}\mathopen{\left({\bf k}_{q}\right)}}$ yields the
strong-contrast approximation for macroscopically isotropic media:

$\displaystyle\frac{\mathinner{\varepsilon_{e}\mathopen{\left({\bf
k}_{q}\right)}}}{\varepsilon_{q}}$
$\displaystyle=1+\frac{d{\beta_{pq}}{\phi_{p}}^{2}}{\phi_{p}(1-\beta_{pq}\phi_{p})-\beta_{pq}\mathinner{A_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}}$
$\displaystyle=1+\frac{d{\beta_{pq}}{\phi_{p}}^{2}}{\phi_{p}(1-\beta_{pq}\phi_{p})+\frac{(d-1)\pi\beta_{pq}}{2^{d/2}\mathinner{\Gamma\mathopen{\left(d/2\right)}}}\mathinner{F\mathopen{\left({\bf
k}_{q}\right)}}},$ (67)

where $\beta_{pq}$ is defined in Eq. (33),
$\mathinner{A_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}\equiv\Tr[\mathinner{\boldsymbol{A}_{2}^{(p)}\mathopen{\left({\bf
k}_{q}\right)}}]/d$, and $\mathinner{F\mathopen{\left({\bf Q}\right)}}$ is
what we call the nonlocal attenuation function of a composite for reasons we
describe below. The direct- and Fourier-space representations of
$\mathinner{F\mathopen{\left({\bf Q}\right)}}$ are given as

$\displaystyle\mathinner{F\mathopen{\left({\bf Q}\right)}}$
$\displaystyle\equiv-\frac{2^{d/2}\mathinner{\Gamma\mathopen{\left(d/2\right)}}}{\pi}{Q}^{2}\int_{\epsilon}\frac{i}{4}\quantity(\frac{Q}{2\pi
r})^{d/2-1}$
$\displaystyle\quad\quad\times\mathinner{\mathcal{H}_{d/2-1}^{(1)}\mathopen{\left(Qr\right)}}e^{-i{\bf
Q}\cdot{\bf r}}\mathinner{\chi_{{}_{V}}\mathopen{\left({\bf
r}\right)}}\differential{{\bf r}}$ (68)
$\displaystyle=-\frac{\mathinner{\Gamma\mathopen{\left(d/2\right)}}}{2^{d/2}\pi^{d+1}}{Q}^{2}\int\frac{\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left({\bf
q}\right)}}}{\absolutevalue{{\bf q}+{\bf Q}}^{2}-Q^{2}}\differential{{\bf
q}}.$ (69)

The exponential $\mathinner{\exp\mathopen{\left(-i{\bf Q}\cdot{\bf
r}\right)}}$ in Eq. (68) arises from the phase difference associated with the
incident waves at positions separated by ${\bf r}$. In the quasistatic regime,
this phase factor is negligible, and Eq. (68) reduces to the local attenuation
function $\mathinner{\mathcal{F}\mathopen{\left(Q\right)}}$ (derived in Ref.
Rechtsman and Torquato (2008) and summarized in the Supplementary Material SM
) because it is barely different from unity over the correlation length
associated with the autocovariance function
$\mathinner{\chi_{{}_{V}}\mathopen{\left(r\right)}}$. The strong-contrast
approximation (67) was postulated in Ref. Kim and Torquato (2020) on physical
grounds. By contrast, the present work derives it as a consequence of our
exact nonlocal formalism (Sec. III).

For statistically isotropic media, the effective dielectric constant as well
as the attenuation function are independent of the direction of the incident
wave vector ${\bf k}_{q}$, and thus, they can be considered as functions of
the wave number, i.e.,
$\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}=\mathinner{\varepsilon_{e}\mathopen{\left({\bf
k}_{q}\right)}}$ and
$\mathinner{F\mathopen{\left(k_{q}\right)}}=\mathinner{F\mathopen{\left({\bf
k}_{q}\right)}}$. Then, the real and imaginary parts of Eq. (68) can be
simplified as

$\displaystyle\imaginary[\mathinner{F\mathopen{\left(Q\right)}}]=$
$\displaystyle\begin{cases}-\frac{Q^{2}}{\pi^{2}}\int_{0}^{\pi/2}\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(2Q\cos\phi\right)}}\differential{\phi},&d=2\\\
-\frac{Q}{2(2\pi)^{3/2}}\int_{0}^{2Q}q\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(q\right)}}\differential{q},&d=3\end{cases}$
(70) $\displaystyle\real[\mathinner{F\mathopen{\left(Q\right)}}]=$
$\displaystyle-\frac{2Q^{2}}{\pi}\mathrm{p.v.}\int_{0}^{\infty}\differential{q}\frac{1}{q(Q^{2}-q^{2})}\imaginary[\mathinner{F\mathopen{\left(q\right)}}],$
(71)

where Eq. (71) is valid for $d=2,3$, and $\mathrm{p.v.}$ stands for the Cauchy
principal value. Following conventional usage, we say that a composite
attenuates waves at a given wave number if the imaginary part of the effective
dielectric constant is positive. Recall that attenuation in the present study
occurs only because of multiple-scattering effects (not absorption). While it
is the imaginary part of $F(Q)$ that determines directly the degree of
attenuation or, equivalently, $\imaginary[\varepsilon_{e}]$, we see from Eq.
(71) that the real part of $F(Q)$ is directly related to its imaginary part.
It is for this reason that we refer to the complex function $F(Q)$ as the
(nonlocal) attenuation function.

####  IV.2.2 Modified strong-contrast approximation at the two-point level

Here we extend the validity of the strong-contrast approximation (67) so that
it is accurate at larger wave numbers and hence better captures spatial
dispersion. This is done by an appropriate rescaling of the wave number in the
reference phase, $k_{q}$, which we show is tantamount to approximately
accounting for higher-order contributions in the remainder term
$\mathinner{\mathcal{R}_{2}\mathopen{\left({k}_{q}\right)}}$. Given that the
strong-contrast expansion for isotropic media perturbs around the Hashin-
Shtrikman structures (see Fig. 3) in the static limit, it is natural to use
the scaling $\sqrt{\varepsilon_{\mathrm{HS}}/\varepsilon_{q}}k_{q}$, where
$\varepsilon_{\mathrm{HS}}$ is the Hashin-Shtrikman estimate, i.e.,

$\varepsilon_{\mathrm{HS}}\equiv\varepsilon_{q}\quantity[1+\frac{d\phi_{p}\beta_{pq}}{1-\phi_{p}\beta_{pq}}],$
(72)

which gives the Hashin-Shtrikman lower bound and upper bound if
$\varepsilon_{p}>\varepsilon_{q}$ and $\varepsilon_{p}<\varepsilon_{q}$,
respectively. This scaling yields the following scaled strong-contrast
approximation for statistically isotropic media:

$\frac{\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}}{\varepsilon_{q}}=1+\frac{d{\beta_{pq}}{\phi_{p}}^{2}}{\phi_{p}(1-\beta_{pq}\phi_{p})+\frac{(d-1)\pi\beta_{pq}}{2^{d/2}\mathinner{\Gamma\mathopen{\left(d/2\right)}}}\mathinner{F\mathopen{\left(\sqrt{\frac{\varepsilon_{\mathrm{HS}}}{\varepsilon_{q}}}k_{q}\right)}}}.$
(73)

We now show that the scaled approximation (73) indeed provides good estimates
of leading-order corrections of
$\mathinner{\mathcal{R}_{2}\mathopen{\left({k}_{q}\right)}}$ in powers of
$k_{q}$. To do so, we employ the concept of the averaged (effective) Green’s
function $\left\langle\mathinner{\boldsymbol{G}^{(q)}\mathopen{\left({\bf
q}\right)}}\right\rangle$ of an inhomogeneous medium which in principle
accounts for the all multiple-scattering events

$\displaystyle\left\langle\mathinner{\boldsymbol{G}^{(q)}\mathopen{\left({\bf
q}\right)}}\right\rangle=\quantity(\frac{\omega}{c})^{2}\quantity{\quantity[q^{2}-{\mathinner{k_{e}\mathopen{\left(\omega\right)}}}^{2}]\boldsymbol{I}-{\bf
q}{\bf q}}^{-1},$ (74)

where
${\mathinner{k_{e}\mathopen{\left(\omega\right)}}}\equiv\sqrt{\mathinner{\varepsilon_{e}\mathopen{\left(\omega\right)}}}\omega/c=\sqrt{\mathinner{\varepsilon_{e}\mathopen{\left(\omega\right)}}/\varepsilon_{q}}k_{q}$
is the effective wave number at a frequency $\omega$, and
$\mathinner{\varepsilon_{e}\mathopen{\left(\omega\right)}}$ is the exact
effective dynamic dielectric constant, assuming a well-defined homogenization
description Sheng (1995); Carminati and Sáenz (2009). Since exact complete
microstructural information is, in principle, accounted for with the effective
Green’s function (74), the exact strong-contrast expansion can be
approximately equated to the one truncated at the two-point level with an
attenuation function given in terms of the effective Green’s function, i.e.,

$\displaystyle{\phi_{p}}^{2}\beta^{2}_{pq}\frac{\varepsilon_{e}(k_{q})+(d-1)\varepsilon_{q}}{\varepsilon_{e}(k_{q})-\varepsilon_{q}}$
$\displaystyle=$
$\displaystyle\phi_{p}\beta_{pq}-\mathinner{A_{2}^{(p)}\mathopen{\left(k_{q}\right)}}{\beta_{pq}}^{2}+\mathinner{\mathcal{R}_{2}\mathopen{\left(k_{q}\right)}},$
(75) $\displaystyle\approx$
$\displaystyle\phi_{p}\beta_{pq}-\mathinner{A_{2}^{(p)}\mathopen{\left(\mathinner{k_{e}\mathopen{\left(\omega\right)}}\right)}}{\beta_{pq}}^{2}.$
(76)

When the functional form of $\mathinner{A_{2}^{(p)}\mathopen{\left(Q\right)}}$
or, equivalently, $\mathinner{F\mathopen{\left(Q\right)}}$ is available, it is
possible to solve Eq. (76) for
$\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}$ in a self-
consistent manner. Instead, we show that by assuming
$k_{e}(\omega)\approx\sqrt{\varepsilon_{\mathrm{HS}}/\varepsilon_{q}}k_{q}$,
which results in Eq. (73), we obtain good estimates of the leading-order
corrections of $\mathinner{\mathcal{R}_{2}\mathopen{\left({k}_{q}\right)}}$ in
powers of $k_{q}$. Thus, the scaled approximation (73) provides better
estimates of the higher-order three-point approximation (given in Appendix C)
than the unmodified strong-contrast approximation (60). This can be easily
confirmed in the quasistatic regime from the small-$k_{q}$ expansions of
$\mathinner{A_{2}^{(p)}\mathopen{\left(k_{q}\right)}}$ and
$\mathinner{A_{3}^{(p)}\mathopen{\left(k_{q}\right)}}$ given in Ref. Rechtsman
and Torquato (2008). We also confirm the improved predictive capacity of the
scaled strong-contrast approximation to better capture dispersive
characteristics for ordered and disordered models via comparison to finite-
difference time-domain (FDTD) simulations; see Figs. 8 and 9, and Sec. VIII in
the Supplementary Material SM .

##  V Model Microstructures

Figure 4: Representative images of configurations of the four models of 2D
disordered particulate media described in this section. These include (a)
overlapping spheres, (b) equilibrium packings, (c) class I hyperuniform
polydisperse packings, and (d) stealthy hyperuniform packings. For all models,
the volume fraction of the dispersed phase (shown in black) is
$\phi_{2}=0.25$. Note that (a) and (b) are not hyperuniform. While these
models consist of distributions of particles, both overlapping and
nonoverlapping, the formulas derived in Sec. IV can be applied to any two-
phase microstructure. Indeed, Appendix D describes applications to media with
phase-inversion symmetry.

We consider four models of 2D and 3D disordered media to understand the effect
of microstructure on the effective dynamic dielectric constant, two of which
are nonhyperuniform (overlapping spheres and equilibrium hard-sphere packings)
and two of which are hyperuniform (hyperuniform polydisperse packings and
stealthy hyperuniform packings of identical spheres). The particles of
dielectric constant $\varepsilon_{2}$ are distributed throughout a matrix of
dielectric constant $\varepsilon_{1}$. We also compute the spectral density
for each model, which is the required microstructural information to evaluate
the nonlocal strong-contrast approximations discussed in Sec. VB.

Representative images of configurations of the four aforementioned models of
2D disordered particulate media are depicted in Fig. 4. Note that the degree
of volume-fraction fluctuations decreases from the leftmost image to the
rightmost one.

###  V.1 Overlapping spheres

Overlapping spheres (also called the fully-penetrable-sphere model) refer to
an uncorrelated (Poisson) distribution of spheres of radius $a$ throughout a
matrix Torquato (2002). For such nonhyperuniform models at number density
$\rho$ in $d$-dimensional Euclidean space $\mathbb{R}^{d}$, the autocovariance
function is known analytically Torquato (2002):

$\mathinner{\chi_{{}_{V}}\mathopen{\left(r\right)}}=\exp\quantity(-\rho
v_{2}(r;a))-{\phi_{1}}^{2},$ (77)

where $\phi_{1}=\exp(-\rho\mathinner{v_{1}\mathopen{\left(a\right)}})$ is the
volume fraction of the matrix phase (phase 1), $v_{1}(a)$ is given by (13),
and $\mathinner{v_{2}\mathopen{\left(r;a\right)}}$ represents the union volume
of two spheres whose centers are separated by a distance $r$. In two and three
dimensions, the latter quantity is explicitly given, respectively, by

$\displaystyle\frac{\mathinner{v_{2}\mathopen{\left(r;a\right)}}}{\mathinner{v_{1}\mathopen{\left(a\right)}}}=\begin{cases}\begin{aligned}
2\mathinner{\Theta\mathopen{\left(x-1\right)}}&+\frac{2}{\pi}\Big{[}\pi+x\quantity(1-x^{2})^{1/2}\\\
&-\mathinner{\cos^{-1}\mathopen{\left(x\right)}}\Big{]}\mathinner{\Theta\mathopen{\left(1-x\right)}}\end{aligned},&d=2\\\
2\mathinner{\Theta\mathopen{\left(x-1\right)}}+\quantity(1+\frac{3x}{2}-\frac{x^{3}}{2})\mathinner{\Theta\mathopen{\left(1-x\right)}},&d=3\end{cases}$

where $x\equiv r/2a$, and $\Theta(x)$ (equal to 1 for $x>0$ and zero
otherwise) is the Heaviside step function. For $d=2$ and $d=3$, the particle
phase (phase 2) percolates when $\phi_{2}\approx 0.68$ (Ref. Quintanilla  _et
al._ (2000)) and $\phi_{2}\approx 0.29$ (Ref. Rintoul and Torquato (1997)),
respectively. The corresponding spectral densities are easily found
numerically by performing the Fourier transforms indicated in Eq. (7). In this
work, we apply this model for $\phi_{2}$’s well below the percolation
thresholds.

###  V.2 Equilibrium packings

Another disordered nonhyperuniform model we treat is distributions of
equilibrium (Gibbs) of identical hard spheres of radius $a$ along the stable
fluid branch Hansen and McDonald (2013); Torquato (2002). The structure
factors of such packings are well approximated by the Percus-Yevick solution
Hansen and McDonald (2013); Torquato (2002), which is analytically solvable
for odd values of $d$. For $d=3$, the Percus-Yevick solution gives the
following expression for the structure factor
$\mathinner{S\mathopen{\left(Q\right)}}$ Torquato (2002):

$\displaystyle\mathinner{S\mathopen{\left(Q\right)}}=$
$\displaystyle\Big{(}1-\rho\frac{16\pi
a^{3}}{q^{6}}\Big{\\{}\big{[}24a_{1}\phi_{2}-12(a_{1}+2a_{2})\phi_{2}q^{2}$
$\displaystyle\quad+(12a_{2}\phi_{2}+2a_{1}+a_{2}\phi_{2})q^{4}]\cos(q)$
$\displaystyle+[24a_{1}\phi_{2}q-2(a_{1}+2a_{1}\phi_{2}+12a_{2}\phi_{2})q^{3}\big{]}\sin(q)$
$\displaystyle-24\phi_{2}(a_{1}-a_{2}q^{2})\Big{\\}}\Big{)}^{-1},$ (78)

where $q=2Qa$, $a_{1}=(1+2\phi_{2})^{2}/(1-\phi_{2})^{4}$, and
$a_{2}=-(1+0.5\phi_{2})^{2}/(1-\phi_{2})^{4}$. Using this solution in
conjunction with Eq. (11) yields the corresponding spectral density
${\tilde{\chi}}_{{}_{V}}(Q)$. For $d=2$, we obtain the spectral density from
Monte Carlo generated disk packings Torquato (2002).

###  V.3 Hyperuniform polydisperse packings

Class I hyperuniform packings of spheres with a polydispersity in size can be
constructed from nonhyperuniform progenitor point patterns via a tessellation-
based procedure Kim and Torquato (2019a, b). Specifically, we employ the
centers of 2D and 3D configurations of identical hard spheres in equilibrium
at a packing fraction $0.45$ and particle number $N=1000$ as the progenitor
point patterns. One begins with the Voronoi tessellation Torquato (2002) of
these progenitor point patterns. We then rescale the particle in the $j$th
Voronoi cell $C_{j}$ without changing its center such that the packing
fraction inside this cell is identical to a prescribed value $\phi_{2}<1$. The
same process is repeated over all cells. The final packing fraction is
$\phi_{2}=\sum_{j=1}^{N}\mathinner{v_{1}\mathopen{\left(a_{j}\right)}}/V_{\mathfrak{F}}=\rho\mathinner{v_{1}\mathopen{\left(a\right)}}$,
where $\rho$ is the number density of particle centers and $a$ represents the
mean sphere radius. In the small-$\absolutevalue{{\bf Q}}$ regime, the
spectral densities of the resulting particulate composites exhibit a power-law
scaling $\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left({\bf
Q}\right)}}\sim\absolutevalue{{\bf Q}}^{4}$, which are of class I.

###  V.4 Stealthy hyperuniform packings

Stealthy hyperuniform particle systems, which are also class I, are defined by
the spectral density vanishing around the origin, i.e.,
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left({\bf Q}\right)}}=0$ for
$0<\absolutevalue{{\bf Q}}\leq Q_{\mathrm{U}}$; see Eq.(14). We obtain the
spectral density from realizations of disordered stealthy hyperuniform
packings for $d=2,3$ that are numerically generated via the following two-step
procedure. Specifically, we first generate such point configurations
consisting of $N$ particles in a fundamental cell $\mathfrak{F}$ under
periodic boundary conditions via the collective-coordinate optimization
technique Uche  _et al._ (2004); Batten  _et al._ (2008); Zhang  _et al._
(2015), which amounts to finding numerically the ground-state configurations
for the following potential energy;

$\mathinner{\Phi\mathopen{\left({\bf
r}^{N}\right)}}=\frac{1}{V_{\mathfrak{F}}}\sum_{{\bf
Q}}\mathinner{\tilde{v}\mathopen{\left({\bf
Q}\right)}}\mathinner{S\mathopen{\left({\bf
Q}\right)}}+\sum_{i<j}\mathinner{u\mathopen{\left(r_{ij}\right)}},$ (79)

where

$\mathinner{\tilde{v}\mathopen{\left({\bf
Q}\right)}}=\begin{cases}1,&Q_{\text{L}}<\absolutevalue{{\bf Q}}\leq
Q_{\text{U}},\\\ 0,&\mathrm{otherwise},\end{cases}$ (80)

and a soft-core repulsive term Zhang  _et al._ (2017)

$\mathinner{u\mathopen{\left(r\right)}}=\begin{cases}(1-r/\sigma)^{2},&r<\sigma,\\\
0,&\mathrm{otherwise}.\end{cases}$ (81)

In contrast to the usual collective-coordinate procedure Uche  _et al._
(2004); Batten  _et al._ (2008); Zhang  _et al._ (2015), the interaction (79)
used here also includes a soft-core repulsive energy (81), as done in Ref.
Zhang  _et al._ (2017). Thus, the associated ground-state configurations are
still disordered, stealthy and hyperuniform, and their nearest-neighbor
distances are larger than the length scale $\sigma$ due to the soft-core
repulsion $\mathinner{u\mathopen{\left(r\right)}}$. Finally, to create
packings, we follow Ref. Zhang  _et al._ (2016b) by circumscribing the points
by identical spheres of radius $a<\sigma/2$ under the constraint that they
cannot overlap (See the Supplementary Material SM  for certain results
concerning stealthy “nonhyperuniform” packings in which $Q_{\text{L}}>0$.) The
parameters used to generate these disordered stealthy packings are summarized
in the Supplementary Material SM .

###  V.5 Spectral Densities for the Four Models

Here, we plot the spectral density
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(Q\right)}}$ for the four
models at the selected particle-phase volume fraction of $\phi_{2}=0.25$; see
Fig. 5. From the long- to intermediate-wavelength regimes ($Qa\lesssim 4$),
their spectral densities are considerably different from one another.
Overlapping spheres depart the most from hyperuniformity, followed by
equilibrium packings. Stealthy packings suppress volume-fraction fluctuations
to a greater degree than hyperuniform polydisperse packings over a wider range
of wavelengths. In the small-wavelength regime ($Qa\gg 4$), all four curves
tend to collapse onto a single curve, reflecting the fact that all four models
are composed of spheres of similar sizes.

Figure 5: The spectral density
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(Q\right)}}$ as a function of
dimensionless wave number $Qa$ for the four models of 3D disordered media:
overlapping spheres, equilibrium packings, class I hyperuniform polydisperse
packings, and stealthy hyperuniform packings. In all cases, the volume
fraction of the dispersed phase is $\phi_{2}=0.25$. For hyperuniform
polydisperse packings, $a$ is the mean sphere radius. The other three models
consist of identical spheres of radius $a$. Corresponding graphs of the
spectral densities for the 2D models are provided in the Supplementary
Material SM .

##  VI Results for the Nonlocal Attenuation Function

We report some general behaviors of the nonlocal attenuation function
$\mathinner{F\mathopen{\left({\bf Q}\right)}}$ [cf. Eqs. (68) or (69)] for
nonhyperuniform and hyperuniform media for long and intermediate wavelengths.
We also provide plots of both the real and imaginary parts of
$\mathinner{F\mathopen{\left(Q\right)}}$ for the four models of disordered
two-phase media considered in this work, which depends on wave number $Q$.

The function $\mathinner{F\mathopen{\left(Q\right)}}$ depends on the
microstructure via the spectral density
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(Q\right)}}$. Thus, assuming
that the latter quantity has the power-law scaling
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(Q\right)}}\sim Q^{\alpha}$ as
$Q\to 0$, the asymptotic behavior of $\mathinner{F\mathopen{\left(Q\right)}}$
in the long-wavelength limit ($Q\to 0$) follows as

$\displaystyle\imaginary[\mathinner{F\mathopen{\left(Q\right)}}]$
$\displaystyle\sim\begin{cases}Q^{d},&\text{nonhyperuniform}\\\
Q^{d+\alpha},&\text{hyperuniform}\end{cases},$ $\displaystyle\text{as~{}}Q\to
0,$ (82) $\displaystyle\real[\mathinner{F\mathopen{\left(Q\right)}}]$
$\displaystyle\sim Q^{2},$ $\displaystyle~{}~{}\text{as~{}}Q\to 0,$ (83)

where we use Eqs. (68) and (70). Recall that the exponent $\alpha$ lies in the
open interval $(0,\infty)$ for hyperuniform systems (see Sec. II.4). For
nonhyperuniform systems studied here, we take $\alpha=0$. The reader is
referred to the Supplementary Material SM  for derivations of Eqs. (82) and
(83). Importantly, in the quasistatic regime, the imaginary parts of the
effective dielectric constant for both strong-contrast approximations, (67)
and (73), are determined by the asymptotic behaviors
$\mathinner{F\mathopen{\left(Q\right)}}$ indicated in Eq. (82), i.e.,

$\displaystyle\imaginary[\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}]\sim\imaginary[\mathinner{F\mathopen{\left(k_{q}\right)}}]$
$\displaystyle\sim$
$\displaystyle\begin{cases}{k_{q}}^{d},&\text{nonhyperuniform}\\\
{k_{q}}^{d+\alpha},&\text{hyperuniform}\end{cases},\quad\text{as~{}}k_{q}\to
0.$ (84)

Thus, hyperuniform media are less lossy than their nonhyperuniform
counterparts in the quasistatic regime.

In the case of stealthy hyperuniform media [i.e.,
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(Q\right)}}=0$ for $0\leq
Q<Q_{\mathrm{U}}$], the imaginary part of
$\mathinner{F\mathopen{\left(Q\right)}}$ defined by Eq. (70) is identically
zero (transparent or lossless) for any space dimension $d$ for a range of wave
numbers; specifically,

$\displaystyle\imaginary[\mathinner{F\mathopen{\left(Q\right)}}]=0,~{}~{}\text{for
}0\leq Q<Q_{\mathrm{U}}/2.$ (85)

(The Supplementary Material SM  describes how the local attenuation function
$\mathinner{\mathcal{F}\mathopen{\left(Q\right)}}$ derived in Ref. Rechtsman
and Torquato (2008) generally differs from its nonlocal counterpart.) The
transparency interval in which
$\imaginary[\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}]=0$
predicted by the two strong-contrast approximations [Eqs. (67) and (73)] is
thus given by

$\displaystyle\imaginary[\mathinner{\varepsilon_{e}\mathopen{\left(k_{q}\right)}}]=0,$
$\displaystyle\quad\text{for~{}}\begin{cases}0\leq
k_{q}<\frac{Q_{\mathrm{U}}}{2},&\text{[from Eq. \eqref{eq:modified-strong-
contrast-approximation_2pt}],}\\\ 0\leq
k_{q}<\frac{Q_{\mathrm{U}}}{2(\varepsilon_{\mathrm{HS}}/\varepsilon_{q})^{1/2}},&\text{[from
Eq. \eqref{eq:scaled-2pt-approx}}],\end{cases}$ (86)

where $\varepsilon_{\mathrm{HS}}$ is given in Eq. (72). When
$\varepsilon_{p}>\varepsilon_{q}$, since
$\varepsilon_{\mathrm{HS}}>\varepsilon_{q}$, the scaled approximation
accurately predicts a narrower transparency interval than the unscaled
variant, as verified in Sec. VIII. Interestingly, the transparency interval
obtained from the less accurate formula (67) agrees with the one obtained from
previous simulation results for stealthy hyperuniform “point” scatterers
Leseur  _et al._ (2016), not the finite-sized particles considered here.

Figure 6 shows $\mathinner{F\mathopen{\left(Q\right)}}$ for the four distinct
models of disordered particulate media in $\mathbb{R}^{3}$: overlapping
spheres, equilibrium packings, stealthy hyperuniform packings, and
hyperuniform polydisperse packings. (Its 2D counterpart is provided in the
Supplementary Material SM .) We clearly see that these attenuation functions
exhibit common large-$Q$ behaviors, regardless of the microstructures. From
the quasistatic to the intermediate-wavelength regimes ($Qa<1$), however, the
attenuation characteristics (imaginary parts
$\imaginary[\mathinner{F\mathopen{\left(Q\right)}}]$) are considerably
different from one model to another. For example, stealthy hyperuniform media
are transparent up to a finite wavelength, and hyperuniform polydisperse
packings exhibit much less attenuation than nonhyperuniform systems.

Figure 6: The negative values of (a) the real and (b) imaginary parts of the
nonlocal attenuation function $\mathinner{F\mathopen{\left(Q\right)}}$ as a
function of the dimensionless wave number $Qa$ [defined in Eq. (68)] for the
four models of 3D disordered composite media considered in this paper. The
inset in (b) is the log-log plot of the larger panel. The volume fraction of
the dispersed phase for each model is $\phi_{2}=0.25$. The first three models
consist of identical spheres of radius $a$. For class I hyperuniform
polydisperse particulate media, $a$ is the mean sphere radius.

##  VII Simulation Procedure to Compute Effective Dynamic Dielectric Constant

In Ref. Kim and Torquato (2020), we established preliminary comparisons of
strong-contrast approximation (67) and numerical simulations via the extended
version of the fast-Fourier-Transform-based technique Moulinec and Suquet
(1998); Eyre and Milton (1999). Because of convergence issues, however, in
this paper, we employ a more reliable numerical technique, i.e., the FDTD
method Taflove  _et al._ (2013), using an open-source software package Oskooi
_et al._ (2010). We focus here on particulate media and take the matrix to be
the reference medium (phase 1) and the particles to be the polarized phase
(phase 2).

Figure 7:  Schematic of the general simulation setup for either (a) periodic
or (b) nonperiodic composites consisting of $N$ spheres of radius $a$ in a
matrix. In both cases, Gaussian pulses of electric fields propagate from the
planar sources (shown in red lines) to the packings (shown in black circles).
The wave number (spectrum) of the pulses spans between $\min[k_{1}]$ and
$\max[k_{1}]$. Periodic boundary conditions are applied along all directions,
except for the propagation direction $\hat{{\bf x}}$. The perfectly matched
layers (PML, shown in blue) of thickness $L_{\mathrm{PML}}$ are placed at both
ends of the simulation box to absorb any reflected and transmitted waves. To
estimate the effective dielectric constant as described in step 2 below, we
consider the subregion $V$ (shown in gray) that excludes from the composite
within the simulation box two relatively thin slabs of thickness
$L_{\mathrm{boundary}}$ along the propagation direction.

The general simulation setup is schematically illustrated in Fig. 7. The
simulation procedure for macroscopically isotropic media is carried out in
three steps:

  1. 1.

For a given medium, we obtain the steady-state spatial distributions of
electric field $\mathinner{E_{y}\mathopen{\left({\bf r},\omega\right)}}$ and
electric displacement field $\mathinner{D_{y}\mathopen{\left({\bf
r},\omega\right)}}$ at a given frequency $\omega$ (or the corresponding wave
number $k_{1}$). Specifically, the planar source generates Gaussian pulses of
electric fields that propagate along the $\hat{{\bf x}}$ direction with wave
number $k_{1}$ that spans between $\min[k_{1}]$ and $\max[k_{1}]$. Using the
aforementioned MEEP package Oskooi  _et al._ (2010), we compute time evolution
of electric field $\mathinner{E_{y}\mathopen{\left({\bf r},t\right)}}$ and
electric displacement field $\mathinner{D_{y}\mathopen{\left({\bf
r},t\right)}}$ for a period of time
$6\pi\sqrt{\varepsilon_{1}}/\quantity{c\min[k_{1}]}$, where $c$ is the speed
of light in vacuum. We then compute the temporal Fourier transforms of these
fields inside packings. The values of the simulation parameters (indicated in
Fig. 7) for the 2D and 3D ordered and disordered models studied in this
article are summarized in the Supplementary Material SM .

  2. 2.

At each value of $k_{1}$, we postprocess $\mathinner{E_{y}\mathopen{\left({\bf
r},\omega\right)}}$ and $\mathinner{D_{y}\mathopen{\left({\bf
r},\omega\right)}}$ to estimate the effective dielectric constant
$\mathinner{\varepsilon_{*}\mathopen{\left(k_{1},\omega\right)}}$ of a single
configuration by solving the following self-consistent equation:

$\mathinner{\varepsilon_{*}\mathopen{\left(k_{1}\right)}}=\mathinner{\overline{\varepsilon_{*}}\mathopen{\left(k_{1}\right)}},$
(87)

where

$\displaystyle\mathinner{\varepsilon_{*}\mathopen{\left(k_{1}\right)}}\equiv\frac{\mathinner{\tilde{D}_{y}\mathopen{\left(k_{e},\omega\right)}}}{\mathinner{\tilde{E}_{y}\mathopen{\left(k_{e},\omega\right)}}},\quad\mathinner{\overline{\varepsilon_{*}}\mathopen{\left(k_{1}\right)}}\equiv\quantity(\frac{k_{e}}{\omega/c})^{2},$

$k_{e}$ is a complex-valued effective wave number, and

$\displaystyle\mathinner{\tilde{D}_{y}\mathopen{\left(q,\omega\right)}}\equiv$
$\displaystyle\frac{1}{\absolutevalue{V}}\int_{V}\mathinner{D_{y}\mathopen{\left({\bf
r},\omega\right)}}e^{-iq\hat{{\bf x}}\cdot{\bf r}}\differential{{\bf r}},$
$\displaystyle\mathinner{\tilde{E}_{y}\mathopen{\left(q,\omega\right)}}\equiv$
$\displaystyle\frac{1}{\absolutevalue{V}}\int_{V}\mathinner{E_{y}\mathopen{\left({\bf
r},\omega\right)}}e^{-iq\hat{{\bf x}}\cdot{\bf r}}\differential{{\bf r}},$

where $V$ is a rectangular parallelepiped subregion within the composite
(shown in gray in Fig. 6) that is slightly smaller than the simulation box and
is used to reduce undesired boundary effects. The homogenization task is
carried out by numerically finding the minimizer of
$\absolutevalue{\varepsilon_{*}-\overline{\varepsilon_{*}}}^{2}$ with an
initial guess $\overline{\varepsilon_{*}}=\varepsilon_{\mathrm{HS}}$ via the
Broyden-Fletcher-Goldfarb-Shanno (BFGS) nonlinear optimization algorithm Liu
and Nocedal (1989), where $\varepsilon_{\mathrm{HS}}$ is the Hashin-Shtrikman
estimate given by Eq. (72). Details of step 2 are provided in the
Supplementary Material SM .

  3. 3.

Steps 1-2 are repeated for a sufficient number of configurations for
disordered media. Then, we compute the effective dielectric constant at a
given $k_{1}$ by ensemble averaging
$\mathinner{\varepsilon_{*}\mathopen{\left(k_{1}\right)}}$, i.e.,
$\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}=\left\langle\mathinner{\varepsilon_{*}\mathopen{\left(k_{1}\right)}}\right\rangle$.

It is important to note that we ensure that the run times employed in step 1
are sufficiently long such that the computed effective dielectric constants
achieve stable and accurate steady-state values. We emphasize that the result
$\mathinner{\varepsilon_{*}\mathopen{\left(k_{1}\right)}}$ obtained from Eq.
(87) in step 2 is nonlocal in space because it is calculated from the nonlocal
constitutive relation
$\mathinner{\tilde{D}_{y}\mathopen{\left(k_{e},\omega\right)}}=\mathinner{\varepsilon_{*}\mathopen{\left(k_{1}\right)}}\mathinner{\tilde{E}_{y}\mathopen{\left(k_{e},\omega\right)}}$.
For periodic media, step 3 is unnecessary because all configurations are
identical.

##  VIII Comparison of Simulations of
$\mathinner{\varepsilon_{e}\mathopen{\left({\bf k}_{1}\right)}}$ to Various
Approximations Formulas

In this section, we compare our simulations of the effective dynamic
dielectric constant $\mathinner{\varepsilon_{e}\mathopen{\left({\bf
k}_{1}\right)}}$ for various 2D and 3D ordered and disordered model
microstructures to the predictions of the strong-contrast formulas as well as
to conventional approximations, such as MGA (17) and QCA (18). Most of these
models provide stringent tests of the predictive power of the approximations
at finite wave numbers because they are characterized by nontrivial spatial
correlations at intermediate length scales.

###  VIII.1 2D and 3D periodic media

We first carry out our FDTD simulations for the effective dynamic dielectric
constant $\mathinner{\varepsilon_{e}\mathopen{\left({\bf k}_{1}\right)}}$ of
2D and 3D periodic packings (square and simple-cubic lattice packings), which
necessarily depends on the direction of the incident wave ${\bf k}_{1}$. While
these periodic packings are macroscopically isotropic, due to their cubic
symmetry, they are statistically anisotropic (see Supplementary Material SM
for details). For simplicity, we consider only the case where ${\bf k}_{1}$ is
aligned with one of the minimal lattice vectors, i.e., the $\Gamma$-$X$
direction in the first Brillouin zone. Such periodic models enable us to
validate our simulations because $\varepsilon_{e}({\bf k}_{1})$ also can be
accurately extracted from the lowest two photonic bands that are calculated
via MPB, an open-source software package Johnson and Joannopoulos (2001). The
results from the band-structure calculations and our FDTD simulations show
excellent agreement. In particular, our simulations accurately predict two
salient dielectric characteristics that must be exhibited by periodic
packings: transparency up to a finite wave number associated with the edge of
the first Brillouin zone (i.e., $\imaginary[\varepsilon_{e}]=0$ for
$0\leq\absolutevalue{{\bf k}_{1}}\lesssim\pi$), and resonancelike attenuation
due to Bragg diffraction within the photonic band gap 333Here, the large value
of the imaginary part $\imaginary[\varepsilon_{e}]$ is due to a small
penetration depth of evanescent waves. (i.e., a peak in
$\imaginary[\varepsilon_{e}]$ or, equivalently, a sharp transition in
$\real[\varepsilon_{e}]$ 444The Kramers-Kronig relations (90) and (91) dictate
that a resonance phenomenon in the dielectric response, that is, a sharp peak
in $\imaginary[\varepsilon_{e}]$, must correspond to a sharp transition in
$\real[\varepsilon_{e}]$, and vice versa.). Thus, our numerical homogenization
scheme is valid down to intermediate wavelengths (see the Supplementary
Material SM  for comparison of the band-structure and FDTD computations).

Importantly, while our strong-contrast approximations Eqs. (67) and (73)
account for directionality of the incident waves, the MGA and QCA are
independent of the direction of ${\bf k}_{1}$. In Fig. 8, the FDTD simulation
results are compared with the MGAs for $d=2,3$ [Eqs. (17) and (16)], QCA (18)
for $d=3$, as well as the unscaled and scaled strong-contrast approximations
(67) and (73) for $d=2,3$. While all approximations agree with the FDTD
simulations in the quasistatic regime, the MGA and QCA fail to capture
properly two key features: no loss of energy up to a finite wave number and
resonancelike attenuation in the band gaps. Each strong-contrast approximation
captures both of these salient characteristics. However, it is noteworthy that
the scaled strong-contrast approximation [Eq. (73)] agrees very well with the
FDTD simulations. For contrast ratios $\varepsilon_{2}/\varepsilon_{1}<1$,
FDTD simulations are also in very good agreement with the predictions of
strong-contrast approximations for a wide range of wave numbers, as detailed
in the Supplementary Material SM .

(a)  (b)

Figure 8: Comparison of the predictions of the strong-contrast formulas, Eqs.
(67) and (73), to the Maxwell-Garnett [Eqs. (17) and (16)] and QCA (18)
approximations for the effective dynamic dielectric constant
$\mathinner{\varepsilon_{e}\mathopen{\left({\bf k}_{1}\right)}}$ as a function
of the dimensionless wave number $k_{1}L$ of periodic packings to our
corresponding computer simulation results. We consider (a) 3D simple cubic
lattice and (b) 2D square lattice of packing fraction $\phi_{2}=0.25$ and
contrast ratio $\varepsilon_{2}/\varepsilon_{1}=4$. Here, $k_{1}$ is the wave
number in the reference (matrix) phase along the $\Gamma$-$X$ direction, and
$L$ is the side length of a unit cell.

###  VIII.2 Disordered nonhyperuniform and hyperuniform media

To test the predictive capacity of approximation formulas for
$\varepsilon_{e}(k_{1})$ for disordered media as measured against simulations,
we choose to study two distinctly different models: disordered nonhyperuniform
packings (equilibrium packings) and disordered stealthy hyperuniform
disordered packings [i.e.,
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(Q\right)}}=0$ for $0\leq
Qa<1.5$] for both 2D and 3D. Again, we compare our simulation results to the
MGAs [Eqs. (17) and (16)], QCA (18), strong-contrast approximation (67), and
the scaled counterpart (73). The conventional approximations fail to capture
spatial dispersion effects. Specifically, the MGA neglects any microstructural
information, except for the particle shape, and thus cannot account for long-
range correlations, such as the lossless property of stealthy hyperuniform
media. By contrast, while the QCA formula yields better estimates of
$\imaginary[\varepsilon_{e}]$ for nonhyperuniform systems, it cannot generally
capture the correct transparency characteristics of hyperuniform systems,
e.g., it incorrectly predicts
$\imaginary[\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}]=0$ for
all wave numbers, regardless of whether the medium is stealthy hyperuniform or
nonstealthy hyperuniform; see Fig. 9 (b).

(a)  (b)

Figure 9:  Comparison of the predictions of the strong-contrast formulas Eqs.
(67) and (73) to the MGA(16) and QCA (18) approximations for the effective
dynamic dielectric constant
$\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}$ as a function of
the dimensionless wave number $k_{1}a$ of 3D disordered sphere packings to our
corresponding computer simulation results. We consider (a) equilibrium
packings and (b) stealthy hyperuniform packings
[$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(Q\right)}}=0$ for $0\leq
Qa<1.5$] of sphere radius $a$, packing fraction $\phi_{2}=0.25$, and phase-
contrast ratio $\varepsilon_{2}/\varepsilon_{1}=4$. Here, $k_{1}$ is the wave
number in the reference (matrix) phase, and the error bars in the FDTD
simulations represent the standard errors over independent configurations.

On the other hand, the scaled strong-contrast approximation provides excellent
estimates of $\varepsilon_{e}(k_{1})$ for both disordered models, even for
large wave numbers ($0\leq k_{1}a\leq 1$); see Fig. 9. Moreover, the
predictions of both strong-contrast approximations accurately capture the
salient microstructural differences between the nonhyperuniform and
hyperuniform models because they incorporate spatial correlations at finite
wavelengths via the spectral density
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left(Q\right)}}$. For example,
they properly predict that stealthy hyperuniform media are lossless up to a
finite wave number, even if at different cutoff values; see Eq. (86).
Corresponding 2D results are presented in the Supplementary Material SM
because they are qualitatively the same as the 3D results.

##  IX Predictions of Strong-Contrast Approximations for Disordered
Particulate Media

Having established the accuracy of the scaled strong-contrast approximation
(73) for ordered and disordered media in the previous section, we now apply it
to the four different disordered models discussed in Sec. V in order to study
how $\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}$ varies with the
microstructure. We first study how
$\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}$ varies with $k_{1}$
at a fixed contrast ratio $\varepsilon_{2}/\varepsilon_{1}=10$ for the four
models; see Fig. 10. According to Eq. (84), nonhyperuniform and hyperuniform
media in the quasistatic regime have the different scalings, i.e.,
$\imaginary[\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}]\sim{k_{1}}^{d}$
and
$\imaginary[\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}]\sim{k_{1}}^{d+\alpha}$,
respectively, where $\alpha>0$ for hyperuniform systems. This implies that
hyperuniform media are less lossy than their nonhyperuniform counterparts as
$k_{1}$ tends to zero, as seen in the insets of Fig. 10. Moreover, beyond the
quasistatic regime, each model exhibits “effective” transparency for a range
of wave numbers that depends on the microstructure. For 2D and 3D models,
hyperuniform polydisperse packings tend to be effectively transparent for a
wide range of wave numbers compared to the nonhyperuniform ones, while the
stealthy hyperuniform systems are perfectly transparent for the widest range
of wave numbers, as established in Eq. (86) and Sec. VIII. For each model, the
“effective” transparency spectral range must be accompanied by normal
dispersion [i.e., an increase in
$\real[\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}]$ with
$k_{1}$] Agranovich and Ginzburg (1984) because our strong-contrast
approximation is consistent with the Kramers-Kronig relations (see Appendix
B). Moreover, we see that anomalous dispersion [i.e., a decrease in
$\real[\mathinner{\varepsilon_{e}\mathopen{\left(k_{1}\right)}}]$ with
$k_{1}$] occurs at wave numbers larger but near the respective transition
between the effective transparency and appreciable attenuation, which again is
dictated by the Kramers-Kronig relations. The specific anomalous dispersion
behavior is microstructure dependent.

(a)  (b)

Figure 10: Predictions of the scaled strong-contrast approximation (73) for
the effective dynamic dielectric constant $\varepsilon_{e}(k_{1})$ as a
function of the dimensionless wave number $k_{1}a$ of the four models of
disordered media at volume fraction $\phi_{2}=0.25$ and contrast ratio
$\varepsilon_{2}/\varepsilon_{1}=10$: (a) three dimensions and (b) two
dimensions. The inset in the lower panel is the log-log plot of the larger
panel.

(a)  (b)

Figure 11:  Predictions of the strong-contrast approximation (73) for the
effective dynamic dielectric constant $\varepsilon_{e}(k_{1})$ as function of
dielectric-contrast ratio $\varepsilon_{2}/\varepsilon_{1}$ for the four
disordered models, as per Fig. 10, at packing fraction $\phi_{2}=0.25$ and
wave number (a)
$\sqrt{\varepsilon_{\mathrm{HS}}/\varepsilon_{1}}~{}k_{1}a=0.7$ in three
dimensions and (b)
$\sqrt{\varepsilon_{\mathrm{HS}}/\varepsilon_{1}}~{}k_{1}a=0.6$ in two
dimensions.

We now examine how the imaginary part $\imaginary[\varepsilon_{e}]$ varies
with the contrast ratio $\varepsilon_{2}/\varepsilon_{1}$ for the disordered
models for a given large wave number $k_{1}$ inside the transparency interval
for 2D and 3D stealthy hyperuniform systems. These results are summarized in
Fig. 11. The disparity in the attenuation characteristics across
microstructures widens significantly as the contrast ratio increases. Clearly,
overlapping spheres are the lossiest systems. Hyperuniform polydisperse
packings can be nearly as lossless as stealthy hyperuniform ones.

We also study the effect of packing fraction $\phi_{2}$ on the effective phase
speed $\mathinner{c_{e}\mathopen{\left(k_{1}\right)}}$ and effective
attenuation coefficient $\mathinner{\gamma_{e}\mathopen{\left(k_{1}\right)}}$,
as defined by Eqs. (55) and (56), respectively. For concreteness, we focus on
3D stealthy hyperuniform packings. We first generate such packings at a
packing fraction $\phi_{2}=0.4$ and $Q_{\mathrm{U}}a=1.5$, as described in
Sec. V.4. Without changing particle positions, we then shrink particle radii
to attain a packing fraction $\phi_{2}=0.25$, whose stealthy regime is now
$Q_{\mathrm{U}}a\approx 1.33$. The coefficients
$\mathinner{c_{e}\mathopen{\left(k_{1}\right)}}$ and
$\mathinner{\gamma_{e}\mathopen{\left(k_{1}\right)}}$ for these packings with
$\varepsilon_{2}/\varepsilon_{1}=4$ are estimated from the scaled
approximation (73); see Fig. 12. It is seen that the waves propagate
significantly more slowly through the denser medium due to an increase in
multiple-scattering events. Moreover, the transparency intervals (wave-number
ranges where the effective attenuation coefficients vanish) are larger for the
packing with the higher stealthy cutoff value $Q_{U}a=1.5a$ ($\phi_{2}=0.4$),
as predicted by Eq. (86).

Figure 12:  Predictions of the scaled strong-contrast approximation (73) for
the effective wave speed $c_{e}$ and the negative of the attenuation
coefficient $\gamma_{e}$ as a function of the dimensionless wave number
$k_{1}a$ for 3D stealthy hyperuniform sphere packings of contrast ratio
$\varepsilon_{2}/\varepsilon_{1}=4$ at two different packing fractions:
$\phi_{2}=0.4$ with $Q_{\mathrm{U}}a=1.5$ and $\phi_{2}=0.25$ with
$Q_{\mathrm{U}}a\approx 1.33$. The inset is a magnification of the lower
panel.

##  X Discussion

All previous closed-form homogenization estimates of the effective dynamic
dielectric constant apply only at long wavelengths (quasistatic regime) and
for very special macroscopically isotropic disordered composite
microstructures, namely, nonoverlapping spheres or spheroids in a matrix. In
this work, we lay the theoretical foundation that enables us to substantially
extend previous work in both its generality and applicability. First, we
derive exact homogenized constitutive relations for the effective dynamic
dielectric constant tensor $\boldsymbol{\varepsilon}_{e}({\bf k}_{q})$ that
are nonlocal in space from first principles. Second, our strong-contrast
representation of $\boldsymbol{\varepsilon}_{e}({\bf k}_{q})$ exactly accounts
for complete microstructural information (infinite set of $n$-point
correlation functions) for arbitrary microstructures and hence multiple
scattering to all orders for the range of wave numbers for which our extended
homogenization theory applies, i.e., $0\leq|{\bf k}_{q}|\ell\lesssim 1$ (where
$\ell$ is a characteristic heterogeneity length scale). Third, we extract from
the exact expansions accurate nonlocal closed-form approximate formulas for
$\boldsymbol{\varepsilon}_{e}({\bf k}_{q})$, relations (67) and (73), which
are resummed representations of the exact expansions that incorporate
microstructural information through the spectral density
$\mathinner{\tilde{\chi}_{{}_{V}}\mathopen{\left({\bf Q}\right)}}$, which is
easily ascertained for general microstructures either theoretically,
computationally, or via scattering experiments. Depending on whether the high-
dielectric phase percolates or not, the wide class of microstructures that we
can treat includes particulate media consisting of identical or polydisperse
particles of general shape (ellipsoids, cubes, cylinders, polyhedra) with
prescribed orientations that may or not overlap, and cellular networks as well
as media without well-defined inclusions (Sec. III.2). Our approximations
account for multiple scattering across a range of wave numbers. Fourth, we
carry out precise full-waveform simulations for various 2D and 3D models of
ordered and disordered media to validate the accuracy of our nonlocal
microstructure-dependent approximations for wave numbers well beyond the
quasistatic regime.

Having established the accuracy of the scaled strong-contrast approximation
(73), we then apply it to four models of 2D and 3D disordered media (both
nonhyperuniform and hyperuniform) to investigate the effect of microstructure
on the effective wave characteristics. Among other findings, we show that
disordered hyperuniform media are generally less lossy than their
nonhyperuniform counterparts. We also find that our scaled formula (73)
accurately predicts that disordered stealthy hyperuniform media possess a
transparency wave number interval
$[0,0.5\,Q_{\mathrm{U}}(\varepsilon_{\mathrm{HS}}/\varepsilon_{q})^{-1/2})$
[cf. (86)], where most nonhyperuniform disordered media are opaque. Note that,
using multiple-scattering simulations, Leseur et al. Leseur  _et al._ (2016)
were the first to show that stealthy hyperuniform systems should exhibit a
transparency interval, but for “point” scatterers, not the finite-sized
scatterers considered here. Interestingly, their transparency-interval
prediction coincides with the one predicted by our less accurate strong-
contrast formula [cf. Eq. (86)].

The accuracy of our nonlocal closed-form formulas has important practical
implications, since one can now use them to accurately and efficiently predict
the effective wave characteristics well beyond the quasistatic regime of a
wide class of composite microstructures without having to perform
computationally expensive full-blown simulations. Thus, our nonlocal formulas
can be used to accelerate the discovery of novel electromagnetic composites by
appropriate tailoring of the spectral densities and then constructing the
corresponding microstructures by using the Fourier-space inverse methods Chen
and Torquato (2018). For example, from our findings in the present study, it
is clear that stealthy disordered particulate media can be employed as low-
pass filters that transmit waves “isotropically” up to a selected wave number.
Moreover, using the spectral densities of the type found by Chen and Torquato
Chen and Torquato (2018) for stealthy hyperuniform packings (characterized by
a peak value at $Q=Q_{\text{U}}$ with intensities that rapidly decay to zero
for larger wave numbers) and formula (73), one can design materials with
refractive indices that abruptly change over a narrow range of wave numbers.
Of course, one could also explore the design space of effective wave
properties of nonhyperuniform disordered composite media for potential
applications.

Previously, disordered media were often described using cluster expansions
Frisch (1968); Sheng (1995); Tsang (2001), while ordered media were often
studied through dispersion relations and band-structures calculations.
Accordingly, it has been of primary importance to bridge the gap between the
treatments of ordered and disordered to better understand the optical
properties of correlated media. Thus, our work represents an initial step
toward a unified theory to describe the effective optical properties of both
ordered and disordered microstructures over a wide range of incident
wavelengths.

There are a variety of directions for future research. First, our formalism
can be straightforwardly extended to heterogeneous materials composed of more
than two phases or continuous media. Second, it is also of interest to extend
our formalism to applications and phenomena (e.g., magnetic effects) relevant
to the smaller wavelengths noted for metamaterials Linden  _et al._ (2004);
Smith  _et al._ (2004); Wang  _et al._ (2012).

##  Appendix A Different Expansions as a Result of Different Exclusion-Region
Shapes

To get a sense of how the resulting expansions change due to the choice of the
exclusion-region shape, we consider the aforementioned oriented spheroidal
exclusion region in the two limiting disk-like and needle-like cases.
Comparing the expansion parameters in the limit cases given in Eq. (32) to the
strong-contrast expansion with a spherical exclusion region given in Eq. (54),
one can obtain the counterparts of Eq. (54) with disk-like and needle-like
exclusion regions. Specifically, we replace parameters in Eq. (54) according
to the following mappings:

$\displaystyle\beta_{pq}\to$
$\displaystyle(\varepsilon_{p}-\varepsilon_{q})/(d\varepsilon_{p}),$
$\displaystyle\frac{\mathinner{{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}+(d-1)\varepsilon_{q}{}}{\mathinner{{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}-\varepsilon_{q}{}}\to$
$\displaystyle\frac{d\mathinner{{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}}{\mathinner{{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}-\varepsilon_{q}{}},$ $\displaystyle~{}~{}~{}~{}\text{disk-
like},$ $\displaystyle\beta_{pq}\to$
$\displaystyle(\varepsilon_{p}-\varepsilon_{q})/(d\varepsilon_{q}),$
$\displaystyle\frac{\mathinner{{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}+(d-1)\varepsilon_{q}{}}{\mathinner{{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}-\varepsilon_{q}{}}\to$
$\displaystyle\frac{d\varepsilon_{q}{}}{\mathinner{{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}-\varepsilon_{q}{}},$ $\displaystyle~{}~{}~{}~{}\text{needle-
like},$

resulting in the following expansions, respectively,

$\displaystyle{\phi_{p}}^{2}\quantity(\frac{\varepsilon_{p}-\varepsilon_{q}}{d\varepsilon_{p}})^{2}\frac{d\mathinner{{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}}{\mathinner{{\varepsilon}_{e}\mathopen{\left({\bf
k}_{q}\right)}}-\varepsilon_{q}{}}=$