Effective Elastic Wave Characteristics of Composite Media

We derive exact expressions for effective elastodynamic properties of two-phase composites in the long-wavelength (quasistatic) regime via homogenized constitutive relations that are local in space. This is accomplished by extending the"strong-contrast"expansion formalism that was previously applied to the static problem. These strong-contrast expansions explicitly incorporate complete microstructural information of the composite via an infinite set of $n$-point correlation functions. Utilizing the rapid-convergence properties of these series expansions (even for extreme contrast ratios), we extract accurate approximations that depend on the microstructure via the spectral density, which is easy to compute or measure for any composite. We also investigate the predictive power of modifications of such approximation formulas postulated elsewhere [J. Kim and S. Torquato, Proc. Nat. Acad. Sci. ${\bf 117}$, 8764 (2020)] to extend their applicability ${\it beyond ~the ~quasistatic ~regime.}$ The accuracy of these nonlocal microstructure-dependent approximations is validated by comparison to full-waveform simulation results for certain models of dispersions. We apply our formulas to a variety of models of nonhyperuniform and hyperuniform disordered composites. We demonstrate that hyperuniform systems are less lossy than their nonhyperuniform counterparts in the quasistatic regime, and stealthy hyperuniform media can be perfectly transparent for a wide range of wavenumber. Finally, we discuss how to utilize our approximations for engineering composites with prescribed elastic wave characteristics.


Introduction
The theoretical determination of the effective elastic wave characteristics of multiphase composite media is of great importance in geophysics [1,2,3], exploration seismology [4,5], diagnostic sonography [6], crack diagnosis [7,8], architectural acoustics [9] and acoustic metamaterials [10], among many examples. Such effective elastic properties generally depend on the phase properties, phase volume fractions φ i , frequency ω or wavenumber k I of the incident elastic waves, and an infinite set of correlation functions that characterizes the composite microstructure [11,12,13]. There have been numerous theoretical/computational attempts to estimate the effective elastic wave characteristics [11,12,13,14,15,16,17,18]. However, the preponderance of previous closed-form approximation formulas for the effective elastodynamic properties apply only in the quasistatic regime [13,16], i.e., applicable when k I ℓ ≪ 1, where ℓ is a characteristic heterogeneity length scale ‡, and under restrictive conditions. One such closed-form approximation is the Gaunaurd-Überall approximation [19,16], which we employ to compare to simulation data and our nonlocal formulas described below.
Our focus in this paper is the theoretical determination of the effective dynamic stiffness tensor C e (k I , ω) of a two-phase elastic composite in d-dimensional Euclidean space R d , which depends on the frequency ω or wavevector k I of the incident elastic waves beyond the quasistatic regime; see figure 1. From this effective property, one can determine the corresponding effective wave speeds c L,T e and attenuation coefficients γ L,T e . To achieve this goal, we first generalize the strong-contrast expansion formalism that has been employed to treat the static elastic problem [20,21,22] to the elastodynamic problem in the quasistatic regime by establishing homogenized constitutive relations that are local in space. Because of the interplay between longitudinal and transverse waves and the complexity of the fourth-rank tensors that are involved, this task is considerably more challenging than the derivation of its electromagnetic counterparts [23,24]. The terms of the resulting quasistatic strong-contrast expansions are explicitly given in terms of integrals over products of Green's functions and the n-point correlation functions S (i) n (x 1 , · · · , x n ) of the random two-phase medium to infinite order. Here, the quantity S (i) n (x 1 , · · · , x n ) gives the probability of finding n points at positions x 1 , · · · , x n simultaneously in phase i. This implies that multiple scattering to all orders is exactly treated in the long-wavelength or quasistatic regime. It is noteworthy that the strong-contrast expansions are given in terms of expansion parameters that are rational functions of the phase moduli. This endows strong-contrast expansions with rapid convergence properties, even for large phase contrast ratios. This behavior is to be distinguished from standard perturbation treatments that result in so-called "weakcontrast" expansions [22] that slowly converge and only apply for small phase contrast ratios.
Due to the fast-convergence properties of strong-contrast expansions, their lowerorder truncations yield accurate closed-form approximate formulas for the effective dynamic moduli that apply to a wide class of microstructures. Postulated nonlocal variants of these formulas are resummed representations of the strong-contrast ‡ Some multiple-scattering approximations for effective elastic waves are accurate beyond the quasistatic regime; see Ref. [18] and references therein. However, these formulas require complicated scattering coefficients of individual scatterers. expansions that still accurately capture multiple scattering to all orders via the microstructural information embodied in the spectral densityχ V (Q). The quantitỹ χ V (Q) is the Fourier transform of the autocovariance function χ V (r) ≡ S (i) where r ≡ x 2 − x 1 , which can be easy to ascertain for general microstructures theoretically, computationally, or via scattering experiments [25]. We also verify the accuracy of the postulated approximations via full-waveform simulations for certain benchmark models. This validation allows us to use them to predict the effective elastic wave characteristics accurately well beyond the quasistatic regime for a wide class of composite microstructures without computationally expensive full-blown simulations. As discussed in section 2, such a broad microstructure class includes particulate composites consisting of identical or polydisperse particles of arbitrary shapes (ellipsoids, cylinders, polyhedra) that may or not overlap, cellular networks as well as systems without well-defined inclusions. Such broad applicability is a notable advantage of our formulas over other multiple-scattering approximations, such as Keller's approximation [14,26] §. Thus, our postulated formulas can be employed to accelerate the discovery of novel elastodynamic composites by appropriate tailoring of the spectral densities [27,28] and then generating the microstructures satisfying them [28], as elaborated in section 7.
While our strong-contrast formulas for the effective dynamic elastic moduli can be applied to periodic two-phase media, the primary applications are spatially correlated disordered media because they can provide advantages over periodic ones with high crystallographic symmetries [29,30], including perfect isotropy and robustness against defects [31,32]. We are interested in both "garden-variety" models [22] as well as exotic hyperuniform forms [33,34,35] of disordered two-phase media. Hyperuniform two-phase systems are characterized by an anomalous suppression of volume-fraction fluctuations in the infinite-wavelength limit [33,34,35], i.e., the spectral densityχ V (Q) obeys the condition Such hyperuniform two-phase media encompass all periodic systems, many quasiperiodic media, and exotic disordered ones; see Ref. [35] and references therein. Disordered hyperuniform systems are exotic states of matters that lie between crystals and liquids; they behave like crystals in the way they suppress large-scale density fluctuations and yet are like liquids because they are statistically isotropic without any Bragg peaks [33,34,35]. Hyperuniform systems have attracted considerable attention over the last decade because of their close connections to a broad spectrum of topics that arise in physical [29,35,36,37,38,39,40,41,42,43,44,45,46,47,48], mathematical [49,50,51], and materials sciences [28,52,53,54,55] as well as the emerging § Keller's approximation is derived for the simplified case in which only longitudinal waves propagate in a very special system: colloidal suspensions of spherical particles in which the fluid has zero shear modulus. Such systems can be treated with the scalar Helmholtz equation, which is to be contrasted with our treatment of the full elastodynamic equations for macroscopically anisotropic media.
technological importance of the disordered varieties [32,35,45,56,52,57,58,59,60]. We apply our nonlocal strong-contrast formulas to predict the real and imaginary parts of the effective elastic moduli for model microstructures that possess some typical disorder (nonhyperuniform) as well as those with exotic hyperuniform disorder (section 5). We are particularly interested in exploring the elastic 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,61,62,63], i.e., where Q ≡ |Q|. Disordered stealthy hyperuniform materials have been shown to exhibit novel optical, acoustic, mechanical, and transport properties [23,53,54,57,64,65,66,67]. Among other results, we show here that disordered hyperuniform media are generally less lossy than their nonhyperuniform counterparts. We also demonstrate that disordered stealthy hyperuniform particulate composites exhibit novel wave characteristics, including the capability to act as low-pass filters that transmit elastic waves isotropically without loss up to a selected wavenumber. Our results demonstrate that one can design the effective wave characteristics of a disordered composite, hyperuniform or not, by engineering spatial correlations of microstructure at prescribed length scales.
In section 2, we present the strong-contrast formalism to derive corresponding expansions of the effective elastic wave characteristics of macroscopically anisotropic two-phase media in the quasistatic regime. While we assume that both phases are elastically isotropic for simplicity, the effective elastic properties are described by a full fourth-rank tensor (what we mean by macroscopically anisotropic media) due to possibly statistically anisotropic microstructures. In section 3, we extract strong-contrast approximations from the exact expansions. In section 4, we extend the validity of the strong-contrast approximations for the effective dynamic moduli so that they apply well beyond the quasistatic regime. The accuracy of these nonlocal approximations is verified by comparison to full-waveform simulations for certain benchmark models. In section 5, we describe four models of disordered composites that we treat in the paper, two of which are nonhyperuniform and two of which are hyperuniform. In section 6, we investigate the microstructure-dependence of the effective elastic wave characteristics for these models. Finally, we provide concluding remarks in section 7.

Exact Strong-Contrast Expansions
Here we extend the general strong-contrast formalism that was devised for the purely static elastic problem [20,21,22] to the elastodynamic problem in the long-wavelength (quasistatic) regime. We first present a compact derivation of the expansions for the effective stiffness tensor C e (k I , ω) of a macroscopically anisotropic medium (section 2.1) and then specialize them to a macroscopically isotropic medium (section 2.2). Detailed derivations are given in the Supplementary Material (SM) [69].  Figure 1: (a) Schematic of a large ellipsoidal, macroscopically anisotropic two-phase composite medium embedded in an infinite reference phase of mass density ρ I and stiffness tensor C I (gray regions) under an applied elastic waves ǫ 0 (x) =ǫ 0 exp(i(k I · x − ωt)) of frequency ω. The wavelength λ associated with the applied wave can span from the quasistatic regime (2πℓ/λ ≪ 1) down to the intermediate-wavelength regime (2πℓ/λ 1), where ℓ is a characteristic length scale. (b) After homogenization, the same ellipsoid can be regarded to be a specimen of a homogeneous medium with an effective stiffness tensor C e (k I , ω), which depends on ω and k I . As noted in the main text, we omit the ω dependence of C e (k I , ω) because (without loss of generally) we assume a linear dispersion relation between |k I | and ω. In the infinite-volume limit, we show that the effective wave characteristics are independent of the shape of the ellipsoidally-shaped composite.
We will exploit the same useful mathematical properties of the strong-contrast formalism that has been used to treat elastostatics [21,22], as we briefly outline here. We begin with the integral solution of the elastodynamic equations in terms of the fourth-rank tensor Green's function. The singular nature of the Green's function requires us to exclude a region around the singularity, but we recognize that the choice of the shape of this "exclusion" region enables us to generate an infinite family of exact series expansions. By choosing a spherical exclusion region, we are able to derive strong-contrast expansions that rapidly converge, even for large phase contrast ratios. The terms of the resulting strong-contrast expansions are explicitly given in terms of absolutely convergent integrals over products of Green's functions and certain n-point correlation functions through all orders. The rapid convergence of strongcontrast expansions enables us to extract accurate approximation formulas from the exact expansions.

Macroscopically Anisotropic Media
Here, we consider macroscopically anisotropic two-phase media with isotropic phases but whose effective elastic properties are described by a full fourth-rank tensor C e . Macroscopic anisotropy arises with isotropic phases because the microstructure can generally be statistically anisotropic, e.g., layered media and oriented ellipsoids in a matrix. (See Ref. [68] for a description of anisotropic phases for elastostatics.) We follow closely the strong-contrast formalism of Torquato [21,22] but apply it to derive the analogous series expansions for the effective dynamic moduli. We consider a macroscopically large ellipsoidal specimen of two-phase statistically homogeneous but anisotropic composite in R d embedded inside an infinitely large reference phase I with mass density ρ I and stiffness tensor C I ; see figure 1. The microstructure is perfectly general, and its inhomogeneity length scale ℓ is much smaller than the specimen size, i.e., ℓ ≪ L. The shape of this specimen is purposely chosen as 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. Mathematically, we prove this below by showing that the strong-contrast formalism leads to effective properties that involve absolutely convergent integrals.
For a two-phase medium, we define the indicator function for phase i as [22,70] For statistically homogeneous media, its ensemble average is simply the phase volume fraction, i.e., φ i ≡ I (i) (x) so that φ 1 + φ 2 = 1. The local stiffness tensor C (x) of such a medium can be written as where C i denotes the stiffness tensor of phase i (= 1, 2). For simplicity, we take the reference phase to be phase q (equal to 1 or 2). In the ensuing discussion, we make the following three assumptions on the phase properties.
(a) Phase 1 and phase 2 are elastically isotropic in d-dimensional Euclidean space R d , i.e., where K i and G i are bulk and shear moduli of phase i (= 1, 2), respectively. Here, the hydrostatic projection tensor Λ h and shear projection tensor Λ s are constant fourth-rank tensors given by where δ ij is the Kronecker delta symbol. The tensor Λ h projects onto fields that are isotropic everywhere, whereas the tensor Λ s projects onto fields that are trace-free (see the SM [69] for useful identities).
(b) Each phase is dissipationless, namely, the elastic moduli K i and G i for i = 1, 2 are real-valued and frequency-independent.
(c) The mass densities of both phases are identical, i.e., Assumption (a) enables us to decompose the elastic waves in phase i (= 1, 2) into longitudinal and transverse waves with their respective wave speeds c L i and c T i , which are defined by Under this assumption, the Poisson ratio ν i of phase i is expressed as [22] and bounded in the range of −1 ≤ ν i ≤ 1/(d − 1) [22,71]. Assumption (b) means that these speeds are independent of frequency ω, implying the following linear dispersion relations: where k L i and k T i are longitudinal and transverse wavenumbres, respectively. Assumption (c) is achievable for many pairs of solid materials; see discussion in Ref. [67]. We suppose that the applied or incident elastic strain field ǫ 0 (x) is a plane wave of an angular frequency ω, well-defined propagation directionk in the reference phase, and the associated wavelength λ. Our interest is in deriving an exact expression for the effective stiffness tensor C e (ω) or, equivalently, C e k Lq in the quasistatic regime (ℓ ≪ λ), where k Lq is the longitudinal wavenumber (11) in the reference phase. While each phase is dissipationless, as stated in (b), the composite is generally lossy (i.e., C e is complex-valued) due to scattering from the inhomogeneities. Nonetheless, our results can be straightforwardly extended to viscoelastic media (with complex-valued moduli), but this will not be done in the present work.
In this paper, 'wave speed' always refers to the speed associated with the phase of the wave. This term is used instead of 'phase speed' because 'phase' in this paper refers to a constituent material of a composite.
Under the assumptions (a)-(c), the local displacement field u(x) solves the timeharmonic wave equation: where the Einstein summation is implied, c Lq and c T q are given in (9), and P ij (x) is the induced stress polarization field given by and ǫ (x) is the local strain tensor [71]. The symmetric second-rank tensor P (x) is the induced field relative to the reference phase q, and hence is nonzero only in the "polarized" phase p (p = q). Following Torquato [21,22], we use a Green's function formalism to solve the wave equation (12) for u for an arbitrary macroscopically anisotropic two-phase medium: where u 0 (x) is related to the applied strain ǫ 0 , and g (q) (r) is the the second-rank Green's function associated with (12). Taking the symmetric part of the gradient of u(x) gives an integral equation for the local strain tensor ǫ (x): where the fourth-rank Green function G (q) (r) associated with the reference phase q is given by [22] is a constant fourth-rank tensor that arises when one excludes an infinitesimal volume around the singularity of the Green function at x ′ = x, and H (q) (r) is the contribution outside of this "exclusion" region. The fourth-rank tensor H (q) (r) is symmetric under the following index changes, i.e., and its explicit expression is given by where r L ≡ k Lq r, r T ≡ k T q r, I is the fourth-rank identity tensor, and H ν (x) is the Hankel function of the first kind of order ν. The three fourth-rank tensors T i (r) for i = 1, 2, 3 are defined, in component form, as r ≡ r/ |r|, andr i is the ith component ofr. Formulas for the traces of H (q) (r) are provided in the SM [69]. Note that (18) reduces to its static counterpart given in Refs. [20,21,22] up to a multiplicative factor ρ q in the static limit. The constant tensor D (q) depends on the "exclusion-region" shape. Due to the fast-convergence properties of the resulting expansion discussed below and in Refs. [20,22,24], we choose a sphericalexclusion region, for which The integral equation (15) is written in a compact linear operator form as Excluding the contribution from the exclusion region in (23), we define generalized cavity strain tensor: Use of (23), (24), and (13) demonstrates that P and f are directly related as follows: where κ pq and µ pq are the scalar polarizabilities for bulk and shear moduli, respectively, defined by Note that Eqs. (22) and (26) are identical to their static counterparts [21,22,72] up to a multiplicative factor ρ q . We now find a series expansion for the following homogenized constitutive relation where · denotes an ensemble average, and the effective constant tensor L (q) e k Lq is explicitly given as To do so, we solve P (x) in terms of ǫ 0 by iteratively substituting (24) and (25). We then obtain the relation (29) from the aforementioned expansion and an ensemble average of (24) by eliminating ǫ 0 in order to avoid conditional convergence problems. Following the strong-contrast formalism of Torquato [21,22], we obtain an expression of the effective tensor L (q) e in the form of series expansion: where n is a position-dependent determinant involving up to the n-point correlation function associated with the dispersed phase p, i.e., Here, S n (x 1 , · · · , x n ) is the n-point correlation function defined as which gives the probability for simultaneously finding n points at x 1 , x 2 , · · · , x n in phase p [22,70]. Here, it is important to note that the integrals (32) and (33) are absolutely convergent because while H (q) (r) decays as r −d for large r, ∆ n identically vanishes at the boundary of the specimen [22]. Therefore, this proves that the effective elastodynamic properties in the infinite-volume limit are independent of the shape of the ellipsoidal composite shown in Fig. 1. The detailed derivation of the strong-contrast expansion (31) is given in the SM. Importantly, the exact series expansion (31) accounts for complete microstructural information (infinite set of n-point correlation functions) and hence multiple scattering to all orders in the quasistatic regime.

Remarks:
(i) Importantly, the strong-contrast expansion (31) is a series representation of a linear fractional transformation of the effective stiffness tensor C e k Lq (left-hand side).
The series expansion in powers of the polarizabilities κ pq and µ pq of this particular rational function of C e k Lq has important consequences for the predictive power of approximations derived from the expansion. While this desirable feature is briefly discussed below, the reader is referred to Ref. [24] for detailed explanations for the corresponding electromagnetic problem. We note that the strong-contrast formalism is a significant departure from standard perturbative expansions that lead to "weak-contrast" expansions in which the expansion parameters are simple differences in the phase moduli, implying that they converge only for small contrast ratios [21,22].
(ii) The homogenized constitutive relation (29) is local in space [i.e., P (x) at point x depends on f (x) at the same position x] and strictly valid in the long-wavelength regime. In such a regime, the effective elastic moduli are independent of the direction of incident waves, as shown in the expansion (31). For shorter wavelengths, however, the associated relation must be nonlocal in space [i.e., P (x) at point x depends on f (x ′ ) at positions around x], which can result in "wavevector"dependent effective elastic moduli, as was rigorously shown for the analogous electromagnetic wave problem [24].
(iii) Note that the expansion (31) represents two different series: one for q = 1 and p = 2 and the other for q = 2 and p = 1.
(v) In contrast to its static counterpart, the n-point microstructure-dependent tensors B (p) n k Lq given in Eqs. (32) and (33) are functions of a frequency ω of the elastic waves or, equivalently, the longitudinal and transverse wavenumbers k Lq and k T q . Throughout this work, we use k Lq as an independent variable, instead of ω or k T q , for the following three reasons. First, the tensor H (q) (r) as well as the effective stiffness tensors are conveniently written in terms of k Lq . Second, k Lq is directly proportional to ω and k T q [cf. (9)]. Furthermore, k Lq is directly related to a length scale, which is suitable for describing microstructural information rather than the temporal quantity ω.
(vi) The exact expansions (31) are independent of the reference phase q and hence of the associated wavenumber k Lq .
(vii) Note that the strong-contrast formalism for the elastodynamic problem shares similar mathematical structure to the electromagnetic counterpart [23,24]. In both cases, the wave equations can be simplified to the Helmholtz equation [i.e., (∇ 2 + k 2 ) u(x) = 0], which results in integral operator descriptions of their expansions being formally identical. However, there are important fundamental distinctions between the two problems. Among other things, while electromagnetic waves have only transverse propagation modes, elastic waves always have both transverse and longitudinal modes with different wave speeds. The interplay between these two propagation modes makes the theoretical determination of the effective elastodynamic properties generally more complex than its electromagnetic counterpart.

Macroscopically Isotropic Media
Here we assume that the composite is macroscopically isotropic. In this case, the effective stiffness tensor C e can be expressed in the effective bulk and shear moduli (denoted by K e and G e , respectively). Then, the series expansion (31) can be reduced to Utilizing the properties of two tensors Λ h and Λ s (see the SM [69]), one can separate (36) into two expansions by taking the quadruple inner products of Λ h and Λ s with (36). One is associated with the effective bulk modulus, and the other is related to the effective shear modulus: respectively, where C Assuming that the composite is passive (i.e., it does not generate mechanical energy), and the time-harmonic factor of waves is e −iωt , the imaginary parts of the effective elastic moduli must be non-positive, implying that for any non-negative k Lq . In light of these properties, we have The effective elastic wave characteristics, including wave speeds c L,T e and attenuation coefficients γ L,T e , are directly related to the effective moduli as follows: where ρ e = ρ p = ρ q , and the superscripts L and T denote longitudinal and transverse waves, respectively. Note that exp(−2πγ L e /c L e ) and exp(−2πγ T e /c T e ) represent the factors by which the amplitudes of the incident waves are attenuated inside the composite for a period of time 2π/ω. Thus, if γ L e = γ T e = 0 at certain wavenumbers (or frequencies), the composite is perfectly transparent, i.e., elastic waves propagate without any loss.

Remarks:
(i) Any statistically isotropic medium is macroscopically isotropic, but the converse is not true. For example, while cubic lattice packings are statistically anisotropic, they are macroscopically isotropic due to the cubic symmetry (see section 4.2).
(ii) The dynamic strong-contrast expansions represented by (37) and (38) possess fast-convergence properties for a wide class of microstructures, even at extreme phase contrast ratios (see Ref. [24] for detailed explanations). Such convergence properties are attributed to the following two aspects. First, even for extreme contrast ratios K p /K q or G p /G q , the two expansion parameters κ pq and µ pq are rational functions of the phase moduli and bounded by where ν q is the Poisson ratio of the reference phase q. Secondly, as Torquato [20,21] observed, the strong-contrast expansions in the static limit can be regarded to be ones that perturb around the wide class of optimal structures [20,21,73], including the optimal multiscale Hashin-Shtrikman "coated-spheres" assemblages. The reader is referred to Refs. [20,21] for details. It suffices to note here that such optimal two-phase media are characterized by a disconnected dispersed phase that is distributed throughout a connected matrix. These observations imply that the first few terms of the expansions (37) and (38) can yield accurate approximations of the effective properties for a class of particulate composites as well as more general microstructures, even for extreme contrast ratios, provided that the dispersed phase is prevented from forming large clusters compared to the specimen size. Depending on whether the high-stiffness phase percolates or not, 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 [73] as well as media without well-defined inclusions. The reader is referred to Ref. [24] for a more detailed discussion of this issue.

Approximations at the Two-Point Level
Due to the fast-convergence properties of strong-contrast expansions, their truncations at low orders should yield accurate approximations for the effective bulk and shear moduli for the aforementioned wide class of microstructures over a broad range of volume fractions and contrast ratios; see also Ref. [24] for additional details. In what follows, we present such approximations by truncating the strong-contrast expansions after the second-order term. The corresponding approximations at the three-point level are presented in Appendix A. Detailed derivations are provided in Sec. I in the SM [69]. Truncating (37) and (38) at the two-point level and solving the left-hand sides of these truncated series for K e and G e , respectively, yields where C k Lq are defined respectively as and the local attenuation function F (Q) is defined as where Γ(x) is the gamma function, H ν (x) is the Hankel function of the first kind of order ν, χ V (r) ≡ S (p) 2 (r) − φ p 2 is the radial autocovariance function, and the spectral densityχ V (Q) is its Fourier transform. Some important properties of F (Q) are given in Appendix C. Use of these properties of F (Q) immediately shows that in the static limit (ω = 0), the parameters C 2 (0) are identically zero, which is consistent with previous studies [72,22].
Remarkably, F (Q) also appears in the quasistatic strong-contrast approximations for the electromagnetic characteristics [23,24]. This commonality between the two wave problems at the two-point level allowed us to establish cross-property relations for the effective elastic and electromagnetic wave characteristics in Ref. [67].

Improved Approximations at the Two-point Level
In order to extend the series expansions and approximations discussed in section 2.1section 3 beyond the quasistatic regime, one needs to generalize the strong-contrast expansion formalism to theories of elastodynamics that are nonlocal in space (see a recent review [74]) from first principles, as we did for the electrodynamic problem in Ref. [24]. Unlike the dielectric problems, however, such generalizations are nontrivial in the case of the elastodynamic problem because an elastically isotropic medium generally possesses multiple elastic wavenumbers at a given frequency ω.

Nonlocal Strong-Contrast Approximation
Based on the following two observations, we postulate nonlocal strong-contrast approximations for the effective elastodynamic properties at the two-point level that are expected to be accurate beyond the quasistatic regime. First, the local strongcontrast expansions for the elastodynamic and electromagnetic problems are similar in that the local attenuation function F (Q), given by (47), appears in the local strongcontrast approximations of the effective dielectric constant that was rigorously derived in Ref. [23]; see also Ref. [24]. Second, guided by our exact formulation of the nonlocal effective electromagnetic characteristics [24], such generalizations at the twopoint level are tantamount to replacing the wavenumber-dependent local attenuation function F (Q), defined in (47), with the wavevector-dependent nonlocal attenuation function F (Q) defined by [24,67] Unlike F (Q), F (Q) accounts for the contribution from spatial variation of the sinusoidal incident waves exp(−iQ · r) and thus more accurately estimates the scattering effects of waves associated with wavevector Q from the long-to intermediate-wavelength regimes.
(Important properties of F (Q) for a statistically isotropic medium are provided in Appendix C.) From these two observations, it is reasonable to assume that one can extend the range of applicable wavelengths by replacing F (Q) in the local strongcontrast approximations at the two-point level [Eqs. (43) and (44)] with F (Q), which are numerically verified in section 4.2. The resulting approximations are given respectively by where k Lq and k Tq are the longitudinal and transverse wavevectors of the incident waves, respectively, and F (Q) is given in (49). We emphasize that the nonlocal strongcontrast approximations for both elastic and electromagnetic properties share a common microstructure-dependent parameter F (Q), which enabled us to establish cross-property relations linking those properties in Ref. [67]. Note that, as we shown in a recent  2)] for the effective dynamic bulk K e k L 1 and shear G e k L 1 moduli as functions of dimensionless waveunmber k L 1 L of periodic packings to the corresponding simulation results. We consider 3D cubic lattice packing of packing fraction φ 2 = 0.05, contrast ratios K 2 /K 1 = G 2 /G 1 = 2, and Poisson ratio ν 1 = 1/3. Here, k L 1 is the longitudinal wavenumber in the reference phase (phase 1) along the Γ-X direction, and L is the nearest-neighbor distance.
paper [24], the analytic properties of F (|Q|) lead the nonlocal approximations (51) and (52) satisfies Kramers-Kronig relations for elastic waves [75,76]. (These nonlocal approximations were first postulated in Ref. [67] on physical grounds for establishing the cross-property relations.) For a statistically isotropic composite, as shown in (51) and (52), the imaginary part of F (|Q|) directly determines the degree of attenuation, i.e., Im[K e ] and Im[G e ] or, equivalently, γ L e and γ T e defined in (41) and (42). In the quasistatic regime, assuming that the spectral density has the power-law scalingχ V (Q) ∼ Q α , the effective attenuation coefficients γ L,T e k Lq exhibit where nonhyperuniform systems take α = 0, whereas hyperuniform ones take α > 0 (see Appendix C). Thus, hyperuniform media are less lossy than their nonhyperuniform counterparts as the wavenumber tends to zero. Remarkably, the stealthy hyperuniform media are perfectly transparent up to a finite wavenumber: where c T q / c Lq = (1 − 2ν q )/[2(1 − ν q )], and ν q is the Poisson ratio of the reference phase q.

Comparison of Simulations to Various Approximations
Here we compare various approximations formulas for the effective dynamic elastic moduli to computer simulations, which are highly nontrivial calculations. In particular, we utilize our fast-Fourier-transform (FFT) numerical scheme presented elsewhere [67]. This procedure extends the one first devised for the effective static elastic moduli [77] in order to treat elastodynamics. The reader is referred to the SM [69] and Ref. [67] for details.
In order to ensure convergence of the simulation procedure, we choose to study simple cubic lattice packings in which identical spheres of radius a of phase 2 are embedded in the matrix phase (phase 1). While the periodic packings are macroscopically isotropic, due to cubic symmetry, they are statistically anisotropic, implying that effective properties can depend on the direction of the incident wave k L1 . For simplicity, we only consider the case where k L1 is aligned with one of the minimal lattice vectors, i.e., Γ-X direction in the first Brillouin zone. Simple cubic lattice packings also provide stringent tests of the predictive power of the approximations at finite wavenumbers because they exhibit two salient and nontrivial elastic properties due to spatial correlations at intermediate length scales: transparency up to finite wavenumbers associated with the edges of the first Brillouin zone (i.e., Im[K e ] = 0 for 0 ≤ k L1 π and Im[G e ] = 0 for 0 ≤ k T 1 π), and resonance-like attenuation due to Bragg diffraction within the phononic bandgap (i.e., a peak in the imaginary parts or, equivalently, a sharp transition in the real parts).
We perform simulations for the case of simple cubic lattice of spheres in a matrix in which the packing fraction is φ 2 = 0.05, contrast ratios are K 2 /K 1 = G 2 /G 1 = 2, and the Poisson ratio of the reference phase is ν 1 = 1/3. In figure 2, we compare the simulation results to the predictions from the strong-contrast approximations [Eqs. (43) and (44) for local approximations, and Eqs. (51) and (52) for the nonlocal counterparts] as well as the Gaunaurd-Überall approximation (GUA) [(B.1) and (B.2)]. While all approximations agree with the simulations in the quasistatic regime, the GUA and local strong-contrast approximations fail to capture properly two key features: no loss of energy up to finite wavenumbers and resonance-like attenuation in the band gaps. However, the nonlocal strong-contrast approximations capture these two features and agree well with the simulation results, even beyond the quasistatic regime.

Disordered Model Microstructures
Here, we describe the four models of 3D disordered two-phase media that are statistically isotropic to study the microstructure-dependence of effective elastic properties. The models include two nonhyperuniform systems (overlapping spheres and equilibrium packings) and two hyperuniform systems (class I hyperuniform polydisperse packings and stealthy hyperuniform packings). In each mode, spherical particles of phase 2 are distributed throughout a matrix phase (phase 1).
Overlapping spheres Equilibrium packing Hyperuniform polydisperse packing Stealthy hyperuniform; Q U a = 1.5 Figure 4: (Color online) The spectral densityχ V (Q) as a function of dimensionless wavenumber Qa for the four models of 3D disordered media at φ 2 = 0.25: overlapping spheres, equilibrium packings, class I hyperuniform polydisperse packings, and stealthy hyperuniform packings. For hyperuniform polydisperse packings, a is the mean sphere radius. The remaining models consist of identical spheres of radius a. Figure 3 depicts representative images of these four models in two dimensions at the selected volume fraction of the disperse phase φ 2 = 0.25 for the purpose of visualization. Note that the degree of volume-fraction fluctuations decreases from figure 3(a) to (d). We also compute the corresponding spectral densities for these models in three dimensions and plot them in figure 4.

Overlapping Spheres
Overlapping spheres (also called fully-penetrable-sphere model) refer to systems of identical spheres of radius a whose centers are spatially uncorrelated in a matrix phase [22]. At a given number density ρ in d-dimensional Euclidean space R d , the autocovariance function of this model can be written analytically [22]. For d = 3, it explicitly writes where φ 1 = exp(−ρ v 1 (a)) is the volume fraction of the matrix phase (phase 1), v 1 (a) = 4πa 3 /3 is the volume of a sphere of radius a, x ≡ r/(2a), and Θ(x) (equal to 1 for x > 0 and zero otherwise) is Heaviside step function. For d = 3, the particle phase (phase 2) percolates when φ 2 0.29 (Ref. [78]). In this work, we consider this model for φ 2 well below the percolation threshold.

Hyperuniform Polydisperse Packings
Class I hyperuniform sphere packings with a polydispersity in size can be constructed from nonhyperuniform progenitor point patterns via a tessellation-based procedure [55,80]. Specifically, we employ the centers of 3D configurations of equilibrium packings (section 5.2) at a packing fraction 0.45 and particle number N = 1000 as the progenitor point patterns. We begin with the Voronoi tessellation [22] 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 resulting packing fraction is where ρ is the number density of particle centers, V F is the volume of the periodic fundamental cell, and a represents the mean sphere radius. In the small-|Q| regime, the spectral density of the resulting particulate composites exhibit a power-law scalingχ V (Q) ∼ |Q| 4 , which are of class I.

Stealthy Hyperuniform Packings
Stealthy hyperuniform packings of identical spheres, which are also class I, are defined by the spectral density vanishing around the origin; see (2). We obtain the spectral density from their realizations for d = 3 that are numerically generated via the following two steps. Specifically, we first generate stealthy hyperuniform point configurations that include N particles in a fundamental cell F under periodic boundary conditions via the collective-coordinate optimization technique [61,62,63], which amounts to finding numerically the ground-state configurations for the following potential energy; where S(Q) is the structure factor of the particle centers, and a soft-core repulsive term [81] u(r) = (1 − r/σ) 2 , r < σ 0, otherwise .
Different from the usual procedure [61,62,63], the interaction (58) used here also includes a soft-core repulsive energy (60), as was done in Ref. [81]. Thus, the resulting configurations are still disordered and stealthy hyperuniform, and their nearest-neighbor distances are larger than the length scale σ due to the soft-core repulsion u(r). Finally, to create packings, we circumscribe the points by identical spheres of radius a < σ/2 under the constraint that they cannot overlap. The parameters used to generate these packings are summarized in the SM [69].

Spectral Densities for the Four Models
Here, we plot the spectral densityχ V (Q) for the fouar models for d = 3 at a selected particle-phase volume fraction φ 2 = 0.25. From the long-to intermediatewavelength regimes (Qa 4), their spectral densities exhibit notable microstructuredependence. For example, overlapping spheres have the largest degree of volume-fraction fluctuations, followed by equilibrium packings. By contrast, in the small-wavelength regime (Qa ≫ 4), all four curves collapse onto a single curve because these models consist of spheres of similar sizes and thus have similar local structures. and γ L,T e , as a function of dimensionless wavenumber k L 1 a for the four 3D models of disordered composites of spheres of radius a and φ 2 = 0.25. The Poisson ratios of the matrix and dispersed phases are ν 1 = 0.4 and ν 2 = 0.25, respectively, and the phase contrast ratios are K 2 /K 1 = 10, G 2 /G 1 = 28, which correspond to glass beads in an epoxy matrix [82].
Here, k L1 is the longitudinal wavenumber in the reference phase, and c L1 and c T 1 are longitudinal and transverse wave speeds, respectively. The insets in the lower panels are log-log plots of the respective larger panels.

Predictions from Strong-Contrast Approximations
Having established the accuracy of the nonlocal strong-contrast approximations, (51) and (52), for simple cubic lattice packings in section 4.2, we now apply them to predict the effective elastodynamic characteristics of the four different disordered models discussed in section 5. Specifically, we study how the effective elastic moduli [K e (k L1 ), G e (k L1 )], wave speeds c L,T e (k L1 ), and attenuation coefficients γ L,T e (k L1 )] vary with the microstructure. For simplicity, we take the matrix phase to be the reference phase (phase 1) and the dispersed phase to be the polarized phase (phase 2). Figure 5 shows the scaled effective wave characteristics [i.e., c L e / c L1 and γ L e /c L e for longitudinal waves and c T e / c T 1 and γ T e /c T e for transverse waves] vary with k L1 at fixed phase properties K 2 /K 1 = 10, G 2 /G 1 = 28, and ν 1 = 0.4 for the four models. While all models are effectively lossless (i.e., small values of γ e ) for a range of wavenumber around the origin, they become increasingly lossy as the wavenumber increases; see the lower panels of figure 5. In the quasistatic regime, as shown in the insets of figure 5, hyperuniform and nonhyperuniform exhibit qualitatively different attenuation characteristics [cf. (53)]: hyperuniform composites generally tend to be less lossy than their nonhyperuniform counterparts. Remarkably, stealthy hyperuniform media can be perfectly lossless, even well beyond the quasistatic regime; see (54). Such microstructuredependence of the effective attenuation behaviors vividly demonstrates that γ L,T e can be engineered by the spatial correlations of composites.
We now examine how the imaginary part Im[G e ] varies with the contrast ratio G 2 /G 1 for the disordered models for a given large wavenumber k Lq inside the transparency interval (wavenumber ranges where the imaginary parts of the effective bulk and shear moduli vanish) given in (54) for the stealthy hyperuniform packing.
Here, we fix the phase Poisson ratios to be ν 1 = 0.4 and ν 2 = 0.25, as we did for the case shown in figure 5. These results are summarized in figure 6. 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. Unlike the imaginary part, the real part Re[G e ] is virtually independent of model microstructure and thus is not shown in this work. We also do not include the corresponding plot for K e because its behavior is qualitatively similar to that of G e . Since stealthy hyperuniform packings exhibit novel physical properties, such as perfect transparency, we further study the effect of packing fraction φ 2 on their effective elastic moduli K e (k L1 ) and G e (k L1 ) and the effective Poisson ratio ν e (k L1 ). Specifically, we are interested in examining stealthy hyperuniform packings consisting of auxetic particles of ν 2 = −1 and a matrix phase with ν 1 = 1/3 and G 2 /G 1 = 10. Auxetic (negative Poisson ratio) materials laterally dilate (shrink)in response to axial elongation (contraction) [83], and are known to have superior energy-absorbing properties [84]. We first generate such packings at a packing fraction φ 2 = 0.4 and Q U a = 1.5, as described in section 5.4. Without changing particle positions, we then shrink the sphere radii to attain a packing fraction φ 2 = 0.25, whose stealthy regime is now Q U a ≈ 1.33.
In figure 7, we plot the effective bulk and shear moduli as well as the effective  Poisson ratio using approximations (51) and (52), and (10). To quantifythe damping characteristics of such composites, we also include in this figure the corresponding loss tangents defined by for some general effective property X e , which are frequently measured in experiments.
For the bulk and shear moduli, the loss tangents represent the ratios of mechanically attenuated energy to the stored elastic energy [85]. We see that these stealthy dispersions are effectively auxetic, i.e., Re[ν e ] < 0 [see figure 7(c)]. Figure 7(a) reveals that such stealthy auxetic composites have exceptionally large loss tangent values in the intermediate-wavelength regime, compared to typical values ( 10 −1 as in the cases in figure 5), which implies that they are excellent energy absorbers, as expected. The transparency interval (i.e., tan δ Ke = tan δ Ge = 0) is slightly larger for the higher density packing with the higher stealthy cut-off value Q U a = 1.5a, as predicted by (54). The complex Poisson ratio implies that the lateral and axial vibrations are out of phase. The absolute value of tan δ νe is approximately proportional to the difference between the shear and bulk loss factors (i.e., degrees of energy loss due to shear and compression); see Ref. [86].

Conclusions and Discussion
Closed-form approximations of the effective dynamic elastic moduli derived previously only apply at long wavelengths (quasistatic regime) and for very special macroscopically isotropic disordered composite microstructures [16], namely, nonoverlapping spheres or spheroids in a matrix. In this work, we have provided the theoretical underpinnings to substantially extend previous work in both its generality and applicability. First, we derived exact homogenized constitutive relations for the effective dynamic elastic stiffness tensor C e k Lq from first principles that are local in space. Second, our strong-contrast representation of C e k Lq exactly accounts for complete microstructural information (n-point correlation functions for n ≥ 1) for general microstructures and hence multiple scattering to all orders in the quasistatic regime. Third, we extracted from the exact expansions accurate local closed-form approximate formulas for K e k Lq and G e k Lq , relations (51) and (52), 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-stiffness 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 (section 2). Fourth, we extended these local approximations beyond the quasistatic regime by postulating nonlocal formulas based on the similarities between electrodynamic and elastodynamic problems and our rigorous formulation of the nonlocal effective dynamic dielectric properties [24]. We carried out precise fullwaveform elastodynamic simulations for certain 3D benchmark models to validate the accuracy of our nonlocal formulas for wavenumbers well beyond the quasistatic regime, i.e., 0 ≤ k Lq ℓ 1 (where ℓ is a characteristic heterogeneity length scale).
Having verified the accuracy of the postulated strong-contrast approximations (51) and (52) for dispersions, we then applied them to the four disordered model microstructures in three dimensions (both nonhyperuniform and hyperuniform) to investigate the microstructure-dependence of the effective elastic wave characteristics. We demonstrated that disordered hyperuniform media are generally less lossy than their nonhyperuniform counterparts. We also found that our approximations predict that disordered hyperuniform media possess a transparency wavenumber interval (54) around which most nonhyperuniform media exhibit strong attenuation. We note that using finite-element method calculations and supercell techniques, Gkantzounis, Amoah, and Florescu [66] showed that 2D stealthy hyperuniform packings should exhibit a transparency interval for elastic waves, which are qualitatively consistent with our predictions.
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 carrying out computationally expensive full-blown simulations. Thus, our nonlocal formulas can be used to accelerate the discovery of novel elastodynamic composites by appropriate tailoring of the spectral densities and then constructing the corresponding microstructures by using Fourier-space inverse methods [28]. For example, from our findings, it is clear that stealthy disordered particulate media can be utilized as low-pass filters that transmit elastic waves "isotropically" up to a selected wavenumber. Of course, one could also explore the design space of effective elastic wave properties of nonhyperuniform disordered composite media for potential applications.
There are interesting open problems for future exploration. Could the exact local strong-contrast expansions, such as (31) and (36), be generalized to the cases in which the mass densities of both phases are different, i.e., ρ 1 = ρ 2 ? This is a highly nontrivial extension. One possible approach to answer this question is to introduce the concept of the dynamic matrix, which is used to derive dispersion relations for elastic waves in simple harmonic lattices [87] in order to separate the local mass density ρ(x) and displacement field u(x). Another challenging problem is the derivation of strongcontrast expansions of the effective dynamic elastic moduli that are nonlocal in space from the first principles in the manner obtained for the electromagnetic problem [24]. This problem is also quite challenging partly because, unlike the electromagnetic waves, one needs to account for the interplay between longitudinal and transverse propagation modes of elastic waves at a given frequency. Finally, it desirable to formulate fullwaveform elastodynamic simulations for two-phase media that are more efficient than the dynamic FFT scheme used here. regardless of whether the composites are hyperuniform or not. The reader is referred to Ref. [24] for derivations. A two-phase composite is effectively lossless for elastic waves at a given frequency ω if and only if Im[K e k Lq ] = 0 and Im[G e k Lq ] = 0, which are equivalent to Im[F k Lq ] = Im[F k T q ] = 0 when using the nonlocal approximations (51) and (52). One can show that these conditions are satisfied in the transparency interval (54) for stealthy hyperuniform media [cf. (2)].