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

We derive exact nonlocal expressions for the effective dielectric constant tensor ${\boldsymbol \varepsilon}_e({\bf k}_I, \omega)$ of disordered two-phase composites and metamaterials from first principles. This formalism extends the long-wavelength limitations of conventional homogenization estimates of ${\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 frequency and wavevector of the incident radiation). This is done by deriving nonlocal strong-contrast expansions that exactly account for multiple scattering for the range of wavenumbers for which our extended homogenization theory applies, i.e., $0 \le |{\bf k}_I| \ell \lesssim 1$ (where $\ell$ is a characteristic heterogeneity length scale). Due to the fast-convergence properties of such expansions, their lower-order truncations yield accurate closed-form approximate formulas for ${\varepsilon}_e({\bf k}_I,\omega)$ that incorporate microstructural information via 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 methods 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 without having to perform full-blown simulations. Among other results, we show that certain disordered hyperuniform particulate composites exhibit novel wave characteristics. 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.


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 [1], Lord Rayleigh [2] and Einstein [3]. 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 [1,[4][5][6], bounding techniques [7][8][9][10][11][12] and exact series-expansion procedures [13][14][15][16]. 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 ε e (k I , ω) of a two-phase dielectric composite, which depends on the frequency ω or wavenumber k I of the incident radiation [17,18]. From this effective property, one can determine the corresponding effective wave speed c e and attenuation coefficient γ e . Virtually all previous homogenization estimates of ε e (k I , ω) apply only in the quasistatic or long-wavelength regime [19] [20], i.e., applicable when |k I |ℓ ≪ 1, where ℓ is a characteristic heterogeneity length scale. This spectral range is the realm of non-resonant dielectric behavior. Examples of such popular closed-form approximation formulas devised for spherical scatterers in a matrix include the Maxwell-Garnett [21,22] and quasicrystalline [18,23,24] estimates, among others [17].
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 ε e (k I , ω) of a macroscopically anisotropic two-phase medium of arbitrary microstructure that are valid well beyond the quasistatic regime, i.e., from the infinite-wavelength limit down to intermediate wavelengths (0 ≤ |k I |ℓ 1). This is accomplished by extending the strong-contrast expansion formalism, which has been used in the past exclusively for the static limit [11] and quasistatic regime [25], and establishing that the resulting homogenized constitutive relations are nonlocal in space (Sec. IV), i.e., the average polarization field at position x depends on the average electric field at other positions around x. (Such nonlocal relations are well-known in the context of crystal optics in order to account for "spatial dispersion," i.e., dependence of dielectric properties on wavevector [26].) 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 A) to infinite order. This implies that multiple scattering to all orders is exactly treated for the range of wavenumbers for which our extended homogenization theory applies, i.e., 0 ≤ |k I |ℓ 1.
Due to the fast-convergence properties of such expansions, elaborated in Sec. IV B, 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. V). 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 to be approximate resummations of the expansions that still accurately capture multiple scattering effects to all orders via the nonlocal attenuation function F (Q) (Sec. VI). The key quantity F (Q) is a functional of the spectral densityχ V (Q) (Sec. II B), 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 microstructuredependent approximations are accurate for both twodimensional (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, cellular networks as well as media without well-defined inclusions (see Sec. IV B for details). Thus, our nonlocal formulas can be employed to accelerate the discovery of novel electromagnetic composites by appropriate tailoring of the spectral densities [27,28] and then generating the microstructures that achieve them [28].
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. III). 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 wavevectors around the origin [28,36,[60][61][62]. Such materials have recently been shown to be endowed with novel optical, acoustic, mechanical, and transport properties [25,48,49,56,[63][64][65][66]. 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 wavenumber. They also can be engineered to exhibit refractive indices that abruptly change over a narrow range of wavenumbers 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.

A. n-Point Correlation Functions
A two-phase random medium is a domain of space V ⊆ R d that is partitioned into two disjoint regions that make up V: a phase 1 region V 1 of volume fraction φ 1 and a phase 2 region V 2 of volume fraction φ 2 [11]. The phase indicator function I (i) (x) of phase i for a given realization is defined as The n-point correlation function S (i) n for phase i is defined by [11,67]: where angular brackets denote an ensemble average over realizations. The function S (i) n (x 1 , . . . , x n ) has a probabilistic interpretation: it gives the probability of finding the ends of the vectors x 1 , . . . , x n all in phase i. For statistically homogeneous media, S (i) n (x 1 , . . . , x n ) is translationally invariant and, in particular, the one-point function is position-independent, i.e., S

B. 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 (r) − 2φ 1 + 1, and hence the autocovariance function is given by which we see is the same for phase 1 and phase 2. Thus, χ V (r = 0) = φ 1 φ 2 and, assuming the medium possesses no long-range order, lim |r|→∞ χ V (r) = 0. For statistically homogeneous and isotropic media, χ V (r) depends only on the magnitude of its argument r = |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 [11] where and Γ(x) is the gamma function. The nonnegative spectral densityχ V (Q), which is proportional to scattering intensity [68], is the Fourier transform of χ V (r), i.e., where Q represents the momentum-transfer wavevector. For statistically homogeneous media, the spectral density must obey the following sum rule [69] For isotropic media, the spectral density only depends on Q = |Q| and, as a consequence of (5), its large-k behavior is controlled by the following power-law form: where is a d-dimensional constant and β(d) is given by (6).

C. Packings
We call a packing in R d a collection of nonoverlapping particles [70]. In the case of a packing of identical spheres of radius a at number density ρ, the spectral densityχ V (Q) is directly related to the structure factor S(Q) of the sphere centers [11,27]: J ν (x) is the Bessel function of the first kind of order ν, φ 2 = ρv 1 (a) is the packing fraction (fraction of space covered by the spheres), and is the d-dimensional volume of a sphere of radius a.

D. Hyperuniformity and Volume-Fraction Fluctuations
Originally introduced in the context of point configurations [33], the hyperuniformity concept was generalized to treat two-phase media [34]. Here the phase volume fraction fluctuates within a spherical window of radius R, which can be characterized by the local volume-fraction variance σ 2 V (R). This variance is directly related to integrals involving either the autocovariance function χ V (r) [71] or the spectral densityχ V (Q) [34].
For typical disordered two-phase media, the variance σ 2 V (R) for large R goes to zero like R −d [11,71,72]. However, for hyperuniform disordered media, σ 2 V (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χ V (Q) is zero for a range of wavevectors around the origin, i.e., where Q U is some positive number. As in the case of hyperuniform point configurations [33][34][35]73], 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χ V (Q) ∼ |Q| α as |Q| → 0: where the exponent α 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. III D.

E. 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 ε e (k 1 ) of isotropic media composed of identical spheres of dielectric constant ε 2 embedded in a matrix phase of dielectric constant ε 1 . In Sec. VIII, we compare the predictions of these formulas to those of our nonlocal approximations. The small-wavenumber expansions of these popular approximations are provided in the SM. There we also provide the corresponding asymptotic behaviors of our strong-contrast approximations.

Maxwell-Garnett approximation
Maxwell-Garnett approximations (MGAs) [21,22] are derived by substituting the dielectric polarizability of a single dielectric sphere into the Clausius-Mossotti equation [22], 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 α e (k 1 ) of a single dielectric sphere of radius a [74]: where 1 (x), the prime symbol ( ′ ) denotes the derivative of a function, h 1 (x) is the spherical Hankel function of the first kind of order 1, and j 1 (x) 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 α e of a dielectric cylinder of radius a given in Ref. [75]: where where H ν (x) is the Hankel function of the first kind of order ν. Here, only 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 spheres and hence are only valid for low inclusion packing fractions. In the static limit, Eqs. (16) and (17) reduce to the Hashin-Shtrikman estimates ε HS ; see relation (77).

Quasicrystalline approximations
The quasicrystalline approximation (QCA) for the quasistatic effective dynamic dielectric constant ε e (k 1 ) employs the "effective" Green's function of spherical scatterers up to the level of the pair correlation function g 2 (r) [18,23,24]. However, this only accounts for the structure factor in the infinite-wavelength limit [24] [i.e., S(0) = 1 + ρ (g 2 (r) − 1) dr] and consequently, spatial correlations at finite wavelengths are ignored. The QCA for d = 3 can be explicitly written as follows [18]: where β 21 is defined in Eq. (38). Interestingly, the QCA predicts that the hyperuniform composites [S(0) = 0] will be transparent for all wavenumbers, which cannot be true for stealthy hyperuniform media (cf. (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. [18] so that it is consistent with the sign convention for the imaginary part Im[ε e ] used here.

III. MODEL MICROSTRUCTURES
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 ε 2 are distributed throughout a matrix of dielectric constant ε 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. 1. Note that the degree of volume-fraction fluctuations decreases from the leftmost image to the rightmost one.

A. Overlapping Spheres
Overlapping spheres (also called fully-penetrablesphere model) refer to an uncorrelated (Poisson) distribution of spheres of radius a throughout a matrix [11]. For such nonhyperuniform models at number density ρ in d-dimensional Euclidean space R d , the autocovariance function is known analytically [11]: where φ 1 = exp(−ρ v 1 (a)) is the volume fraction of the matrix phase (phase 1), v 1 (a) is given by (13), and v 2 (r; a) represents the union volume of two spheres whose centers are separated by a distance r. In two and three dimensions, the latter is explicitly given respectively by where x ≡ r/2a, and Θ(x) (equal to 1 for x > 0 and zero otherwise) is Heaviside step function. For d = 2 and d = 3, the particle phase (phase 2) percolates when φ 2 ≈ 0.68 (Ref. [76]) and φ 2 ≈ 0.29 (Ref. [77]), respectively. The corresponding spectral densities are easily found numerically by performing the Fourier transforms indicated in (7). In this work, we apply this model for φ 2 's well below the percolation thresholds.

B. Equilibrium Packings
Another disordered nonhyperuniform model we treat are distributions of equilibrium (Gibbs) of identical hard spheres of radius a along the stable fluid branch [11,78]. The structure factors of such packings are well approximated by the Percus-Yevick solution [11,78], which is analytically solvable for odd values of d. For d = 3, the Percus-Yevick solution gives the following expression for the structure factor S(Q) [11]: where q = 2Qa, Using this solution in conjunction with (11) yields the corresponding spectral densitỹ χ V (Q). For d = 2, we obtain the spectral density from Monte Carlo generated disk packings [11].

C. 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 [50,79]. 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 [11] of these progenitor point patterns. We then rescale the particle in the jth Voronoi cell C j without changing its center such that the packing fraction inside this cell is identical to a prescribed value φ 2 < 1. The same process is repeated over all cells. The final packing fraction is where ρ is the number density of particle centers and a represents the mean sphere radius. In the small-|Q| regime, the spectral densities of the resulting particulate composites exhibit a power-law scalingχ V (Q) ∼ |Q| 4 , which are of class I.

D. Stealthy Hyperuniform Packings
Stealthy hyperuniform particle systems, which are also class I, are defined by the spectral density vanishing around the origin, i.e.,χ V (Q) = 0 for 0 < |Q| ≤ Q 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 F under periodic boundary conditions via the collective-coordinate optimization technique [60][61][62], which amounts to finding numerically the ground-state configurations for the following potential energy; whereṽ and a soft-core repulsive term [80] u(r) = In contrast to the usual collective-coordinate procedure [60][61][62], the interaction (21) used here also includes a soft-core repulsive energy (23), as done in Ref. [80]. Thus, the associated ground-state configurations are still disordered, stealthy and hyperuniform, and their nearestneighbor distances are larger than the length scale σ due to the soft-core repulsion u(r). Finally, to create packings, we follow Ref. [63] by circumscribing the points by identical spheres of radius a < σ/2 under the constraint that they cannot overlap (See the SM for certain results concerning stealthy "nonhyperuniform" packings in which Q L > 0.) The parameters used to generate these disordered stealthy packings are summarized in the SM.

E. Spectral Densities for the Four Models
Here, we plot the spectral densityχ V (Q) for the four models at the selected particle-phase volume fraction of φ 2 = 0.25 From the long-to intermediate-wavelength regimes (Qa 4), their spectral densities are considerably different from one another. Overlapping spheres departs the most from hyperuniformity, followed by equilibrium packings. Stealthy packings suppress volumefraction fluctuations to a greater degree than hyperuniform polydisperse packings over a wider range of wavelengths. In the small-wavelength regime (Qa ≫ 4), all four curves tend to collapse onto a single curve, reflecting that fact that all four models are composed of spheres of similar sizes.

FIG. 2. (Color online)
Plots of the spectral densityχ V (Q) 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 φ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 SM.

IV. EXACT STRONG-CONTRAST EXPANSIONS AND NONLOCALITY
The original strong-contrast expansions for the effective dynamic dielectric constant obtained by Rechtsman and Torquato [25] 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 [11] that was devised for the purely static problem [11] 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 ε e (k R ) for a macroscopically anisotropic medium of arbitrary microstructure into which a plane wave of wavevector k I is incident. These expressions for ε e (k R ) are valid from the infinite-wavelength limit down to wavelengths (λ = 2π/|k I |) on the order of the heterogeneity length scale ℓ. We explicitly show they necessarily require complete microstructural information, as embodied in the infinite set of n-point correlation functions (Sec. II A) of the composite. We also describe the fast-convergence properties of strong-contrast expansions and their consequences for extracting accurate approximations for ε e (k R ).

A. Strong-Contrast Expansions
Consider a macroscopically large ellipsoidal two-phase statistically homogeneous but anisotropic composite specimen in R d embedded inside an infinitely large reference phase I with a dielectric constant tensor ε I . The microstructure is perfectly general, and it is assumed that the inhomogeneity length scale ℓ is much smaller than the specimen size, i.e., ℓ ≪ 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 E 0 (x) is a plane wave of a angular frequency ω and wavevector k I in the reference phase, i.e., Our interest is the exact expression for the effective dynamic dielectric constant tensor ε e (k I , ω). Without loss of generality, we assume a linear dispersion relation in the reference phase, i.e., k I ≡ |k I | = √ ε I ω/c, where c is the speed of light in vacuum, and thus we henceforth do not explicitly indicate the dependence of functions on ω. The composite is assumed to be nonmagnetic, implying that the phase magnetic permeabilities are identical, i.e., µ 1 = µ 2 = µ 0 . For simplicity, we assume realvalued, frequency-independent isotropic phase dielectric constants ε 1 and ε 2 . Nonetheless, the composite can be generally lossy (i.e., ε 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 will not be done in the present work.
Here we present a compact derivation of strongcontrast expansions. It follows the general formalism of Torquato [11] 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 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 E(x) solves the time-harmonic Maxwell equation [25]: ) of a frequency ω and a wavevector kI at infinity. The wavelength λ associated with the applied field can span from the quasistatic regime (2πℓ/λ ≪ 1) down to the intermediate-wavelength regime (2πℓ/λ 1), where ℓ is the inhomogeneity length scale. (b) After homogenization, the same ellipsoid can be regarded to be a homogeneous specimen with an effective dielectric constant εe(kI , ω), which depends on ω and kI . As noted in the main text, we omit the ω dependence of εe because (without loss of generally) we assume a linear dispersion relation between |kI | and ω.
where P(x) is the polarization field given by and is the local dielectric constant, and I (p) (x) is the indicator function for phase p [cf. (2)]. The vector P(x) 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 = q).
Using the Green's function formalism, the local electric field can be expressed in terms of the following integral equation [11,25]: where the second-rank tensor Green's function G (q) (r) associated with the reference phase q is given by [81] D (q) is a constant second-rank tensor that arises when one excludes an infinitesimal region around the position of the singularity x ′ = x in the Green's function and H (q) (r) represents the contribution outside of the "exclusion" region: wherer ≡ r/|r| is a unit vector directed to r, and H (1) ν (x) is the Hankel function of the first kind of order ν. The Fourier transform of Eq. (29) is particularly simple and concise:G Note that Eq. (31) is independent of the shape of the exclusion region, which stands in contrast to the shapedependent Fourier transform of H (q) ij (r); see SM for details.
Use of (26) and (28) leads to an integral equation for the generalized cavity intensity field F(x): where the integral subscript ǫ indicates that the integral is to be carried out by omitting the exclusion region and then allowing it to uniformly shrink to zero.
Using the definitions (26) and (33), we obtain a linear constitutive relation between P(x) and F(x): where 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 x in R d is chosen to be aligned with the polarization vector P(x) at that position, we have where I is the second-rank identity tensor and A * ∈ [0, 1] is the depolarization factor for a spheroid [11] 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. (35) yields the following shape-dependent tensor: where β pq is the dielectric polarizability defined by Moreover, in these three cases, the generalized cavity intensity field F(x; A * ) reduces to respectively, where D(x) 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 will be elaborated below.
We will now show that our formalism yields an exact relation between the polarization field P(x) and the applied field E 0 (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 (32) as where have temporarily dropped the superscript q. Combination of this equation with (34) yields the following integral equation for the polarization field The desired nonlocal relation is obtained from (41) by successive substitutions: where More explicitly, the nonlocal relation (42) can be expressed as where boldface numbers 1, 2 are short-hand notations for position vectors r 1 , r 2 . Ensemble averaging (44) and invoking statistical homogeneity yields the convolution relation where the operator S depends on relative positions, i.e., S (1, 2) = S (1 − 2), and angular brackets denote an ensemble average. Formally, the nonlocal relation (45) is the same as the one given in Torquato [11] for the static problem, but nonlocality was not explicitly invoked there.
As in the static case [11,15] and quasistatic regime [25], the ensemble-averaged operator S (r), which is given explicitly in the SM in terms of the n-point correlation functions and products of the tensor H(r), depends on the shape of the macroscopic ellipsoidal composite specimen (see Fig. 3). This is due to the fact that H decays like r −d for large r and hence S (r) involves conditionally convergent integrals [11]. To avoid such conditional convergence issues, we follow previous strong-contrast formulations by seeking to eliminate the applied field E 0 in (45) in favor of the average cavity field F (r) in order to get a corresponding nonlocal homogenized constitutive relation between P (r) and F (r) or vice versa. Thus, solving for E 0 in (45) and substituting into the ensemble average of (40) yields Inverting this expression leads to the following nonlocal constitutive relation: where L (q) e (r) is a kernel that is derived immediately below and explicitly given by We are not aware of any previous work that derives such an exact nonlocal homogenized constitutive relation (48) from first principles. Note that the support ℓ s of the kernel L (q) e (r) relative to the incident wavelength λ determines the degree of spatial dispersion. When λ is finite, the relation between P (x) and F (x) in (48) is nonlocal in space. In the regime ℓ s ≪ λ, the nonlocal relation (48) can be well approximated by the local relation . Indeed, in the static limit, L (q) e (r) tends to a Dirac delta function δ(r), expression (48) becomes the position-independent local relation derived earlier [11]. The nonlocal constitutive relation in direct space, Eq. (48), can be reduced to a linear product form in Fourier space by taking the Fourier transform of (48): The wavevector-dependent effective tensor L (q) e (k q ) is postulated (see discussion in the SM) to be given by This linear fractional form for L (q) e (k) is consistent with the one derived for the static limit [11] and for the quasistatic regime [25]. Taking the inverse Fourier transform of (52) yields its corresponding direct-space representation (49). Taking the Fourier transform of (47) yields Comparing (51) to (53), and specifically choosing a spherical exclusion-region, as discussed in Eq. (37), yields the desired exact strong-contrast expansions for general macroscopically anisotropic two-phase media: where A (p) n (k q ) is a wavevector-dependent second-rank tensor that is a functional involving the set of correlation functions S (p) 1 , S (p) 2 , · · · , S (p) n (defined in Sec. II A) and products of the second-rank tensor field H (q) (r), which is explicitly given as [see also Eq. (30)] where k q ≡ |k q |. Specifically, for n = 2 and n ≥ 3, these n-point tensors associated with the polarized phase p are respectively given by where ǫ dr ≡ lim ǫ→0 + |r|>ǫ dr and ∆ (p) n is a position-dependent determinant involving correlation functions of the polarized phase p up to the n-point level: For macroscopically isotropic media, the effective dielectric tensor is isotropic, i.e., ε e (k q ) = ε e (k q ) I. The corresponding strong-contrast expansion for ε e (k q ) is obtained by taking trace of both sides of Eq. (54): where ε e (k q ) = Tr[ε e (k q )]/d, A ó /d for n ≥ 2 and Tr[ ] denotes the trace operation. Furthermore, for statistically isotropic media, the effective dielectric constant becomes independent of the direction of the wavevector, i.e., ε e (k q ) = ε e (k q ).

Remarks:
1. Importantly, the strong-contrast expansion (54) is a series representation of a linear fractional transformation of the variable ε e (k q ) (left-hand side of the equation), rather than the effective dielectric constant tensor itself. The series expansion in powers of the polarizability β pq of this particular rational function of ε e (k q ) has important consequences for the predictive power of approximations derived from the expansion, as detailed in Sec. IV B.
2. The fact that the exact expansion (54), extracted from our nonlocal relation (51), 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 wavenumbers for which our extended homogenization theory applies, i.e., 0 ≤ |k q |ℓ 1.
3. Note that Eq. (54) represents two different series expansions: one for q = 1 and p = 2 and the other for q = 2 and p = 1.
4. The exact expansions represented by (54) are independent of the reference phase q and hence independent of the wavevector k q .
5. For d = 2, the strong-contrast expansion applies for TE polarization only. This implies that the electric field and wavevector are parallel to the plane or transverse to an axis of symmetry in a 3D system whose cross-sections are identical.
6. Formally, the original strong-contrast expansions that apply in the quasistatic regime [25] can be obtained from the nonlocal strong-contrast expansions (54) by simply replacing the exponential functions that appear in the expressions for the secondrank tensors A (p) n (k q ), defined by (57), by unity.
7. For statistically isotropic media, the effective phase speed c e (k q ) and attenuation coefficient γ e (k q ) are determined by the scalar effective dielectric constant ε e (k q ): where n e (k q ) and κ e (k q ) are the effective refractive and extinction indices, respectively. The quantity exp(−2πγ e /c e ) is the factor by which the incidentwave amplitude is attenuated for a period 2π/ω.

B. Convergence Properties and Accuracy of Truncated Series
The form of the strong-contrast expansion parameter β pq in (54) is a direct consequence of the choice of a spherical region excluded from the volume integrals in (54) due to singularities in the Green's functions [25]. It is bounded by which implies that the strong-contrast expansion (54) can converge rapidly, even for infinite contrast ratio ε p /ε q → ∞. 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 [11]. 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 (ε p − ε q )/ε q , and thus is unbounded when ε p /ε q → ∞. Importantly, in the purely static case, the expansion (54) becomes identical to one derived by Sen and Torquato [15] and its truncation after second-order terms [i.e., setting A (p) n = 0 for all n ≥ 3] yields the generalized Hashin-Shtrikman bounds [82] derived by Willis [83] 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 [84]. In the case of an isotropic effective dielectric constant ε e , the optimal Hashin-Shtrikman upper and lower bounds for any phase-contrast ratio ε 2 /ε 1 are exactly realized by the multiscale "coated-spheres" model, which is depicted in Fig. 4 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 [11,85] observed that the strong-contrast expansions (54) in the static limit can be regarded to be ones that perturb FIG. 4. Schematic of the optimal multiscale "coated-spheres" model that realizes the isotropic Hashin-Shtrikman bounds on εe [82]. 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 εe for prescribed volume fraction and contrast ratio, when ε2/ε1 > 1 and ε2/ε1 < 1, respectively. It has recently been proved that these highly degenerate optimal Hashin-Shtrikman multiscale distributions of spheres are hyperuniform [50,79].
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 ε 2 /ε 1 ≫ 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 ε 2 /ε 1 ≪ 1, the matrix phase can be a cellular network [49]. Finally, for moderate values of the contrast ratio ε 2 /ε 1 , even more general microstructures (e.g., those without well-defined inclusions) can be accurately be treated. Importantly, we show that for the dynamic problem under consideration, the first few terms of the expansion (54) yield accurate approximations of ε e (k q ) for a similar wide class of two-phase media (see Sec. VIII). Analogous approximations were derived and applied for the quasistatic regime [25,28].
We now show how lower-order truncations of the series (54) can well approximate higher-order functionals (i.e., higher-order diagrams) of the exact series to all or-ders in terms of lower-order diagrams. Such truncations 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 (54) for ε e yields the rational function in β pq : Expanding (63) in powers of the scalar polarizability β pq yields the series where the first several functionals B (p) n (k q ) are explicitly given in terms of A where T n stands for n successive inner products of a second-rank tensor T . Let us now compare the exact expansion (64) to the one that results when expanding the truncation of the exact expression (54) for [ε e + (d − 1)ε q ][ε e − ε q ] −1 at the 2-point level: where the nth-order functional C (p) n (k q ) for any n is given in terms of the volume fraction φ p and A (p) 2 (k q ), which has the following diagrammatic representation: Here the solid and wavy lines joining two nodes represent the spatial correlation via χ V (r) and a wavevectordependent Green's function H (q) (r) e −ikq·r between the nodes, respectively. The black node indicates a volume integral and carries a factor of dε q . The first several functionals C (p) n (k q ) are explicitly given as Thus, comparing (64) to (66), we see that truncation of the expansion (54) for [ε e + (d − 1)ε q ] · (ε e − ε 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 φ p and the second-order diagram (67). This can be thought of as an approximate but accurate resummed representation of the exact expansion (64). Note that the approximate expansion (66) is exact through second order in β pq . Clearly, truncation of (54) at the three-point level (see Appendix C) will yield even better approximations of the higher-order functionals.

V. STRONG-CONTRAST APPROXIMATION FORMULAS
Here we describe lower-order truncations of the strongcontrast expansions that are expected to yield accurate closed-form formulas for ε e (k q ) that apply over a broad range of wavelengths (k q ℓ 1), volume fractions and contrast ratios for a wide class of microstructures.

A. Macroscopically Anisotropic Media
For the ensuing treatment, it is convenient to rewrite the expansions (54), valid for macroscopically anisotropic media in R d , in the following manner: where the M th-order remainder term is defined as Truncating the exact nonlocal expansion (68) at the twoand three-point levels, i.e., setting R 2 (k q ) = 0 and R 3 (k q ) = 0, respectively, yields Compared to the quasistatic approximation [25], these nonlocal approximations substantially extend the range of applicable wavenumber, namely, 0 ≤ |k q |ℓ 1.

B. 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 ε e (k q ) = Tr [ε e (k q )]/d but depends on the direction of the wavevector k q .

Strong-Contrast Approximation at the Two-Point Level
Solving Eq. (59) for the effective dielectric constant ε e (k q ) yields the strong-contrast approximation for macroscopically isotropic media: where β pq is defined in Eq.
is what we call the nonlocal attenuation function of a composite for reasons described below. The direct-and Fourier-space representations of F (Q) are given as The exponential exp(−iQ · r) in Eq. (73) arises from the phase difference associated with the incident waves at positions separated by r. In the quasistatic regime, this phase factor is negligible, and Eq. (73) reduces to the local attenuation function F (Q) (derived in Ref. [25] and summarized in the SM) because it is barely different from unity over the correlation length associated with the autocovariance function χ V (r). The strong-contrast approximation (72) was postulated in Ref. [66] on physical grounds. By contrast, the present work derives it as a consequence of our exact nonlocal formalism (Sec. IV). For statistically isotropic media, the effective dielectric constant as well as the attenuation function are independent of the direction of the incident wavevector k q , and thus they can be considered as functions of wavenumber, i.e., ε e (k q ) = ε e (k q ) and F (k q ) = F (k q ). Then, the real and imaginary parts of Eq. (73) can be simplified as where Eq. (76) is valid for d = 2, 3. Following conventional usage, we say that a composite attenuates waves at a given wavenumber if the imaginary part of the effective dielectric constant is positive. Recall that attenuation in the present study occurs only because of multiplescattering effects (not absorption). While it is the imaginary part of F (Q) that determines directly the degree of attenuation or, equivalently, Im[ε e ], we see from (76) 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.

Modified Strong-Contrast Approximation at the Two-Point Level
Here we extend the validity of the strong-contrast approximation (72) so that it is accurate at larger wavenumbers and hence better captures spatial dispersion. This is done by an appropriate rescaling of the wavenumber in the reference phase, k q , which we show is tantamount to approximately accounting for higher-order contributions in the remainder term R 2 (k q ). Given that the strongcontrast expansion for isotropic media perturbs around the Hashin-Shtrikman structures (see Fig. 4) in the static limit, it is natural to use the scaling ε HS /ε q k q , where ε HS is the Hashin-Shtrikman estimate, i.e., which gives the Hashin-Shtrikman lower bound and upper bound if ε p > ε q and ε p < ε q , respectively. This yields the following scaled strong-contrast approximation for statistically isotropic media: (78) We now show that the scaled approximation (78) indeed provides good estimates of leading-order corrections of R 2 (k q ) in powers of k q . To do so, we employ the concept of the averaged (effective) Green function¨G (q) (q) ∂ of an inhomogeneous medium which in principle accounts for the all multiple scattering events where k e (ω) ≡ ε e (ω)ω/c = ε e (ω) /ε q k q is the effective wavenumber at a frequency ω and ε e (ω) is the exact effective dynamic dielectric constant, assuming a well-defined homogenization description [17,86]. Since exact complete microstructural information is, in principle, accounted for with the effective Green's function (79), 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 function, i.e., When the functional form of A (p) 2 (Q) or, equivalently, F (Q) is available, it is possible to solve Eq. (81) for ε e (k q ) in a self-consistent manner. Instead, we show that by assuming k e (ω) ≈ ε HS /ε q k q , which results in Eq. (78), we obtain good estimates of the leading-order corrections of R 2 (k q ) in powers of k q . Thus, the scaled approximation (78) provides better estimates of the higher-order three-point approximation (given in Appendix C) than the unmodified strong-contrast approximation (65). This can be easily confirmed in the quasistatic regime from the small-k q expansions of A (p) 2 (k q ) and A (p) 3 (k q ) given in Ref. [25]. 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 FDTD simulations; see Figs. 7 and 8, and Sec. VIII in the SM.

VI. RESULTS FOR THE NONLOCAL ATTENUATION FUNCTION
We report some general behaviors of the nonlocal attenuation function F (Q) [cf. (73) or (74)] for nonhyperuniform and hyperuniform media for long and intermediate wavelengths. We also provide plots of both the real and imaginary parts of F (Q) for the four models of disordered two-phase media considered in this work, which depends on wavenumber Q.
The function F (Q) depends on the microstructure via the spectral densityχ V (Q). Thus, assuming that the latter has the power-law scalingχ V (Q) ∼ Q α as Q → 0, the asymptotic behavior of F (Q) in the long-wavelength limit (Q → 0) follows as where we have used Eqs. (73) and (75). Recall that the exponent α lies in the open interval (0, ∞) for hyperuniform systems (see Sec. II D). For nonhyperuniform systems studied here, we take α = 0. The reader is referred to the 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, (72) and (78), are determined by the asymptotic behaviors F (Q) indicated in Eq. (82), i.e., Thus, hyperuniform media are less lossy than their nonhyperuniform counterparts in the quasistatic regime.
In the case of stealthy hyperuniform media [i.e., χ V (Q) = 0 for 0 ≤ Q < Q U ], the imaginary part of F (Q), defined by (75), is identically zero (transparent or lossless) for any space dimension d for a range of wavenumbers; specifically, (The SM describes how the local attenuation function F (Q) derived in Ref. [25] generally differs from its nonlocal counterpart.) The transparency interval in which Im[ε e (k q )] = 0 predicted by the two strong-contrast approximations [Eqs. (72) and (78)] are thus given by where ε HS is given in Eq. (77). When ε p > ε q , since ε HS > ε 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 (72) agrees with the one obtained from previous simulation results for stealthy hyperuniform "point" scatterers [56], not finite-sized particles considered here. Figure 5 shows F (Q) for the four distinct models of disordered particulate media in R 3 ; overlapping spheres, equilibrium packings, stealthy hyperuniform packings, and hyperuniform polydisperse packings. (Its 2D counterpart is provided in the 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 Im[F (Q)]) are considerably different from one model to another. For example, stealthy hyperuniform media are transparent up to a finite wavelength, and hyperuniform polydisperse packings exhibits much less attenuation than nonhyperuniform systems.

VII. SIMULATION PROCEDURE TO COMPUTE EFFECTIVE DYNAMIC DIELECTRIC CONSTANT
In Ref. [66], we established preliminary comparisons of strong-contrast approximation (72) and numerical simulations via the extended version of the Fast-Fourier-Transform-based technique [87,88]. Due to convergence issues, however, in this paper, we employ a more reliable numerical technique, i.e., the finite-difference timedomain (FDTD) method [89], using an open-source software package [90]. 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).
The general simulation setup is schematically illustrated in Fig. 6. It is carried out in two steps: 1. Obtain the steady-state spatial distributions of electric field E y (r) and electric displacement field D y (r) at a given frequency ω (or the corresponding wavenumber k 1 ). Specifically, the planar source generates Gaussian pulses of electric fields whose wavenumber k 1 spans between min[k 1 ] and max[k 1 ]. Using the aforementioned MEEP package [90], we compute time evolution of electric field E y (r, t) and electric displacement field D y (r, t) for a period of time 6π √ ε 1 /{c min[k 1 ]}. We then compute the temporal Fourier transforms of these fields inside packings. The values of the simulation parameters (indicated in Fig. 6) for the 2D and 3D ordered and disordered models studied in this article are summarized in the SM.
2. Post-process E y (r) and D y (r) to estimate the effective dielectric constant ε e (k 1 , ω) so that it satisfies the following equations: ε e (k 1 ) = ε e (k 1 ) , ε e (k 1 ) =D y (k e ) E y (k e , ω) , where k e represents the effective wavenumber at macroscopic length scales, and where a chosen region V inside the composite. This task is carried by numerically finding a minimizer of |ε e − ε e | 2 with an initial guess ε e = ε HS via the Broyden-Fletcher-Goldfarb-Shanno (BFGS) nonlinear optimization algorithm [91], where ε HS is the Hashin-Shtrikman estimate given by Eq. (77).

VIII. COMPARISON OF SIMULATIONS OF εe(k1) TO VARIOUS APPROXIMATIONS FORMULAS
In this section, we compare our simulations of the effective dynamic dielectric constant ε e (k 1 ) 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 wavenumbers because they are characterized by nontrivial spatial correlations at intermediate length scales.

A. 2D and 3D Periodic Media
We first carry out our FDTD simulations for the effective dynamic dielectric constant ε e (k 1 ) of 2D and 3D periodic packings (square and simple-cubic lattice packings), which necessarily depends on the direction of the incident wave k 1 . While these periodic packings are macroscopically isotropic, due to their cubic symmetry, they are statistically anisotropic (see SM for details). For simplicity, we only consider the case where k 1 is aligned with one of the minimal lattice vectors, i.e., Γ-X direction in the first Brillouin zone. Such periodic models enable us to validate our simulations because ε e (k 1 ) also can be accurately extracted from the lowest two photonic bands that are calculated via MPB, an opensource software package [92]. The results from the bandstructure 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 wavenumber associated with the edge of the first Brillouin zone (i.e., Im[ε e ] = 0 for 0 ≤ |k 1 | π), and resonance-like attenuation due to Bragg diffraction within the photonic bandgap [93] (i.e., a peak in Im[ε e ] or, equivalently, a sharp transition in Re[ε e ] [94]). Thus, our numerical homogenization scheme is valid down to intermediate wavelengths (see the SM for comparison of the band-structure and FDTD computations).
Importantly, while our strong-contrast approximations Eqs. (72) and (78) account for directionality of the incident waves, the MGA and QCA are independent of the direction of k 1 . In Fig. 7, 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 (72) and (78) 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 wavenumber and resonancelike attenuation in the band gaps. Each strong-contrast approximation captures both of these salient characteristics. However, it is noteworthy that the scaled strongcontrast approximation [Eq. (78)] agrees very well with the FDTD simulations. For contrast ratios ε 2 /ε 1 < 1, FDTD simulations are also in very good agreement with the predictions of strong-contrast approximations for a wide range of wavenumbers, as detailed in the SM.  72) and (78), to the Maxwell-Garnett [Eqs. (17) and (16)] and QCA (18) approximations for the effective dynamic dielectric constant εe(k1) of periodic packings to our corresponding computer simulation results. We consider (a) 3D simple cubic lattice and (b) 2D square lattice of packing fraction φ2 = 0.25 and contrast ratio ε2/ε1 = 4. Here k1 is the wavenumber in the reference (matrix) phase along the Γ-X direction, and L is the side length of a unit cell.

B. Disordered Nonhyperuniform and Hyperuniform Media
To test the predictive capacity of approximation formulas for ε 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.,χ V (Q) = 0 for 0 ≤ 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 (72), and the scaled counterpart (78). 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 Im[ε e ] for nonhyperuniform systems, it cannot generally capture the correct transparency characteristics of hyperuniform systems, e.g., it incorrectly predicts Im[ε e (k 1 )] = 0 for all wavenumbers, regardless of whether the medium is stealthy hyperuniform or nonstealthy hyperuniform; see Fig. 8 72) and (78), to the MGA(16) and QCA (18) approximations for the effective dynamic dielectric constant εe(k1) of 3D disordered sphere packings to our corresponding computer simulation results. We consider (a) equilibrium packings and (b) stealthy hyperuniform packings [χ V (Q) = 0 for 0 ≤ Qa < 1.5] of sphere radius a, packing fraction φ2 = 0.25, and phase contrast ratio ε2/ε1 = 4. Here k1 is the wavenumber 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 ε e (k 1 ) for both dis-ordered models, even for large wavenumbers (0 ≤ k 1 a ≤ 1); see Fig. 8. Moreover, the predictions of both strongcontrast 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χ V (Q). For example, they properly predict that stealthy hyperuniform media are lossless up to a finite wavenumber, even if at different cut-off values; see Eq. (86). Corresponding 2D results are presented in the 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 strongcontrast approximation (78) for ordered and disordered media in the previous section, we now apply it to the four different disordered models discussed in Sec. III in order to study how ε e (k 1 ) varies with the microstructure. We first study how ε e (k 1 ) varies with k 1 at a fixed contrast ratio ε 2 /ε 1 = 10 for the four models; see Fig. 9. According to Eq. (84), nonhyperuniform and hyperuniform media in the quasistatic regime have the different scalings, i.e., Im[ε e (k 1 )] ∼ k 1 d and Im[ε e (k 1 )] ∼ k 1 d+α , respectively, where α > 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. 9. Moreover, beyond the quasistatic regime, each model exhibits "effective" transparency for a range of wavenumbers that depends on the microstructure. For 2D and 3D models, hyperuniform polydisperse packings tend to be effectively transparent for a wide range of wavenumbers compared to the nonhyperuniform ones, while the stealthy hyperuniform systems are perfectly transparent for the widest range of wavenumbers, 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 Re[ε e (k 1 )] with k 1 ) [26] because our strongcontrast approximation is consistent with the Kramers-Kronig relations (see Appendix B). Moreover, we see that anomalous dispersion (i.e., a decrease in Re[ε e (k 1 )] with k 1 ) occurs at wavenumbers 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.
We now examine how the imaginary part Im[ε e ] varies with the contrast ratio ε 2 /ε 1 for the disordered models for a given large wavenumber k 1 inside the transparency interval for 2D and 3D stealthy hyperuniform systems. These results are summarized in Fig. 10. The disparity in the attenuation characteristics across microstructures widens significantly as the contrast ratio increases. Clearly, overlapping spheres are the lossiest systems. Hy-  peruniform polydisperse packings can be nearly as lossless as stealthy hyperuniform ones. We also study the effect of packing fraction φ 2 on the effective phase speed c e (k 1 ) and effective attenuation coefficient γ e (k 1 ), as defined by Eqs. (60) and (61), respectively. For concreteness, we focus on 3D stealthy hyperuniform packings. We first generate such packings at a packing fraction φ 2 = 0.4 and Q U a = 1.5, as described in Sec. III D. Without changing particle positions, we then shrink particle radii to attain a packing fraction φ 2 = 0.25, whose stealthy regime is now Q U a ≈ 1.33. The coefficients c e (k 1 ) and γ e (k 1 ) for these packings with ε 2 /ε 1 = 4 are estimated from the scaled approximation (78); see Fig. 11. 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 (wavenumber ranges where the effective attenuation coefficients vanish) is larger for the packing with the higher stealthy cut-off value Q U a = 1.5a  (φ 2 = 0.4), as predicted by Eq. (86).

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 have laid the theoretical foundation that has enabled us to substantially extend previous work in both its generality and applicability. First, we derived exact homogenized constitutive relations for the effective dynamic dielectric constant tensor ε e (k q ) that are nonlocal in space from first principles. Second, our strong-contrast representation of ε e (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 wavenumbers for which our extended homogenization theory applies, i.e., 0 ≤ |k q |ℓ 1 (where ℓ is a characteristic heterogeneity length scale). Third, we extracted from the exact expansions accurate nonlocal closed-form approximate formulas for ε e (k q ), relations (72) and (78), which are resummed representations of the exact expansions that incorporate microstructural information through the spectral densityχ V (Q), 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, cellular networks as well as media without well-defined inclusions (Sec. IV B). Our approximations account for multiple scattering across a range of wavenumbers Fourth, we carried out precise fullwaveform simulations for various 2D and 3D models of ordered and disordered media to validate the accuracy of our nonlocal microstructure-dependent approximations for wavenumbers well beyond the quasistatic regime.
Having established the accuracy of the scaled strongcontrast approximation (78), we then applied it to four models of 2D and 3D disordered media (both nonhyperuniform and hyperuniform) to investigate the effect of mi-crostructure on the effective wave characteristics. Among other findings, we showed that disordered hyperuniform media are generally less lossy than their nonhyperuniform counterparts. We also found that our scaled formula (78) accurately predicts that disordered stealthy hyperuniform media possess a transparency wavenumber interval [0, 0.5 Q U (ε HS /ε q ) −1/2 ) [cf. (86)], where most nonhyperuniform disordered media are opaque. Note that, using multiple-scattering simulations, Leseur, Pierrat, and Carminati [56] were the first to show that stealthy hyperuniform systems should exhibit a transparency interval, but for "point" scatterers, not finite-sized scatterers considered here. Interestingly, their transparency-interval prediction coincides with the one predicted by our less accurate strong-contrast formula [cf. (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 [28]. 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 wavenumber. Moreover, using the spectral densities of the type found by Chen and Torquato [28] for stealthy hyperuniform packings (characterized by a peak value at Q = Q U with intensities that rapidly decay to zero for larger wavenumbers) and formula (78), one can design materials with refractive indices that abruptly change over a narrow range of wavenumbers. Of course, one could also explore the design space of effective wave properties of nonhyperuniform disordered composite media for potential applications.

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. (37) to the strong-contrast expansion with a spherical exclusion-region given in Eq. (59), one can obtain the counterparts of Eq. (59) with disk-like and needle-like exclusion-regions.