Homogenization of piezoelectric planarWillis materials undergoing antiplane shear

Homogenization theories provide models that simplify the constitutive description of heterogeneous media while retaining their macroscopic features. These theories have shown how the governing fields can be macroscopically coupled, even if they are microscopically independent. A prominent example is the Willis theory which predicted the strain–momentum coupling in elastodynamic metamaterials. Recently, a theory that is based on the Green’s function method predicted analogous electro–momentum coupling in piezoelectric metamaterials. Here, we develop a simpler scheme for fibrous piezoelectric composites undergoing antiplane shear waves. We employ a source-driven approach that delivers a unique set of effective properties for arbitrary frequency– wavevector pairs. We numerically show how the resultant homogenized model recovers exactly the dispersion of free waves in the composite. We also compute the effective properties in the long-wavelength limit and off the dispersion curves, and show that the resultant model satisfy causality, reciprocity and energy conservation. By contrast, we show how equivalent models that neglect the electromomentum coupling violate these physical laws. © 2021 Elsevier B.V. All rights reserved.

Central to the field of metamaterials in elastodynamics is the theory of Willis [20,21,22,23], which was developed at the level of generality of random media, and relies on ensemble averaging to handle the randomness. Notably, the theory also applies when periodicity is present as a particular case. The theory shows that the effective response of a material point is nonlocal in space (i.e., depends on the state of other points) and time (i.e., depends on previous states). The theory further shows that the effective velocity may induce an effective stress, and that the effective strain may induce an effective momentum. These nonlocal couplings-which are absent in conventional materials-are now termed the Willis couplings. To date, the theory is still being refined and analyzed by Willis [24,25,26] and others [27][28][29][30][31][32][33][34][35][36][37][38][39]. These theoretical studies are accompanied with experiments that measure and exploit the Willis coupling for applicational purposes [28,[40][41][42][43][44][45][46][47][48].
Our results are presented as follows. Section 2 summarizes the governing equations and introduces our definition of effective fields in periodic piezoelectric media. Section 3 develops the source-driven homogenization scheme for antiplane shear waves of Bloch-Floquet type in two-dimensional piezoelectric composites, using the plane wave expansion method. A numerical study is carried out in Section 4, followed by a summary of our main results and insights in Section 5.

Governing equations and effective fields
The state of a piezoelectric medium is determined from the equations of motion in conjunction with the quasielectrostatic approximation of Maxwell's equations. This approximation is valid when at the considered frequencies, the mechanical wavelengths are much shorter than the electromagnetic wavelengths. Accordingly, the equations are ∇ · σ + f =ṗ, ∇ · D = q, ∇ × E = 0, (1) where f is the body force density, σ is the Cauchy stress, p is the linear momentum, D is the electric displacement, q is the charge density and E is the electric field. It is convenient to identically satisfy Eq. (1) 3 using an electric potential φ, such that E = −∇φ. The state variables of the medium are related through the constitutive relations where u is the displacement field and η is an inelastic strain, the purpose of which is discussed in detail in Refs. [25,32,33,49]; the symbols ρ, A, B, C denote the mass density and the dielectric, piezoelectric and elasticity tensors, which are functions of the position x for a heterogeneous medium. In the sequel, we refer to the quantities in the right-and left column vectors as kinematic-and kinetic variables, respectively. The objective of homogenization theories is to replace the description of the heterogeneous medium by a fictitious homogeneous medium that effectively behaves the same in some appropriate sense. This is carried out by first defining average effective fields, and then determining the effective constitutive relations between them. Here, we focus on time-harmonic Bloch-Floquet waves in periodic media, and adopt the definition of Willis [49], Alù [61] and the references therein for the average fields. Thus, we replace field variables of the form by the effective fields whereζ (x) is periodic with the same periodicity of the medium, Ω is the periodic cell, and ω is the frequency. As discussed in Refs. [25,49], this definition-which filters out the fast oscillations at the microscopic scale-delivers effective variables that satisfy exactly macroscopic balance equations that have the same form as Eq. (1). (See also Refs. [30,61] for more arguments that support this definition.) Our goal is to develop a method for calculating a unique set of effective properties that relates the average fields in two-dimensional piezoelectric composites undergoing antiplane shear waves. To do so, we will use the plane wave expansion method in the averaging process, as described next.

The homogenization scheme
The basic idea in the following homogenization scheme using a plane wave expansion is simple: we represent all the field variables in their Bloch form and expand periodic functions in Fourier space; this allows us to extract from the governing equations the relations between the effective fields. This procedure is described next, for the problem of antiplane shear waves. 3 We consider a piezoelectric composite that is periodic in the (x 1 , x 2 ) plane and uniform in the x 3 direction. The composite is driven by the presence of sources in the form which generate antiplane shear waves in the composite; note that since { s 0 } are constants, it follows that ⟨s⟩ ≡ s. As per the discussion in Section 1, we recall that these mathematical sources need not be experimentally realizable. Specifically, here form (5) demands the eigenstrain to be independent of the phase, and the charge to be controlled in time and spacedifficult requirements to achieve in practice. However, as we show in the sequel, this formulation delivers a unique set of effective properties, satisfying reciprocity and energy conservation; this set, in turn, can describe the overall response of the composite in the presence of arbitrary physical sources. Owing to the form of the prescribed sources, the displacement field has only one component in the direction of x 3 , and this component propagates in the (x 1 , x 2 ) plane; we denote this component by u. We represent the constitutive relations of the composite at each point in the following matrix form: where A, B and C are 2 × 2 matrices, 0 n×m is a n × m zero matrix, and here and henceforth the assumed harmonic dependence on time is used to replace time derivatives with −iω. These quantities enter the reduced form of the governing equations for antiplane shear waves: The first step in the homogenization scheme is to rewrite Eq. (6) using the Bloch ansatz [Eq. (3)], which provides We expand the periodic parts in Eq. (8) in Fourier series: where {G} is the infinite set of reciprocal lattice vectors; if we assume a square lattice with a period a, then { G = 2π a n 1 e 1 + 2π a n 2 e 2 , n 1 , Since Eq. (10) holds for any value of x and t, the expression in the square brackets must vanish, which implies We multiply Eq. (10) by e −iG ′′ ·x , and integrate the result over the unit cell. The orthogonality of the Fourier functions eliminates all terms, except those that satisfy G ′′ where (·) GG ′ denotes the Fourier coefficient of (·) along the basis function e i(G−G ′ )·x .
Having the constitutive equations in the desired form, we proceed to manipulate the governing equations in the same manner. Accordingly, we consider the Bloch form of Eq. (7), namely, Again, since Eq. (13) holds for any value of x and t, the expressions in the brackets must vanish, which allows us to obtain after expanding the periodic parts in Fourier series. We multiply Eq. (14) by e −iG ′′ ·x , and integrate the result over the unit-cell. This provides for each G where Next, we assemble all the matrix equations associated with Eq. (15) for each G to a single infinite matrix system in the form For computational purpose, the size of the matrices and column vectors in Eq. (17) is truncated by defining the number of Fourier components (equivalently, the number of G vectors) to N. In this case, f A is a column vector of 2N components which contains all the Fourier coefficients off G andq G , h A is a column vector of 5N components which contains all the Fourier coefficients ofσ 13 where L A is a 5N × 5N matrix that assembles the Fourier components of the mechanical properties, J A is a 5N × 2N matrix that assembles the Fourier components of the differential operators, w A is a column vector of 2N components that assembles the Fourier components of the displacement and the electric potential and m A is a column vector of 5N components that assembles the Fourier components of the inelastic strain. Explicit expressions for Eqs. (17) and (18) are given in Appendix A. We now extract from Eq. (18) the average fields out of h A , and write them as where J 0 and L 0 are the parts of the matrices that multiply the average fields (associated with G = 0) and (•) s denotes the reduced matrix or vector without the G = 0 terms. 4 Our goal is to express the terms w s usingw andm. To this end, Next, we utilize the equations that do not include G = 0 in Eq. (20) to write where Q s is a square matrix which contain the terms which do not multiply the average fields and T is the part of the matrix D T A L A which multiplies the average fields. (This process is similar to the one carried out in Ref. [57] for flexural waves in Euler-Bernoulli beams.) We are now able to express the fluctuating terms of w using its average value and the properties of the composite, namely, 4 Here L 0 is a 5 × 5 matrix, J 0 is a 5 × 2 matrix,w is a vector with 2 components,m is a vector with 5 components, L s is a 5 × (5N

Table 1
The mass density, dielectric, piezoelectric and elastic moduli of PMMA and PZT4. These materials are used in the computations as the matrix and fiber phases, respectively, of a composite with square lattice of period a = 5 mm. We consider two different PZT4 fibers: one with a circular cross-section of a radius R = 1.9 mm, and one with a circular sector that is defined by the area πR 2 /8.
Eq. (23) is the desired form, since it relates the independent kinematic variables J 0w to the dependent kinetic variables h through a matrix that defines the effective properties. These are given through the following identification The resultant effective properties are nonlocal in time and space, since they are functions of ω and κ. In general, ω and κ are independent and prescribed by the impressed sources; only if the sources are set to zero, then (ω, κ) pairs are related through the dispersion relations that are defined by the normal modes of the composite. Notably, the nonzero blocksW andW † are identified with the electromomentum couplings that were first discovered by Pernas-Salomón and Shmuel [33], using a different formulation that is based on Green's functions. LikeS, the electromomentum coupling is of Willis type, in the sense that it couples fields that are not coupled at the microscopic level, and arises at finite frequencies and/or wavevectors. Symbolically, we write the effective constitutive relations as Refs. [29,30,32,36] discuss how the nonlocal nature of the effective operator generates an ambiguity in the definition of the independent effective variables. They suggest a more physical set by adding the time rate of the original independent effective variables. We adapt this notion, and define an alternative form to Eq. (25), namely, whereŜ = iS/ω,Ŵ = iW/ω and so on. Accordingly, the modified Willis coupling actually depends on the acceleration rather than on the velocity. Likewise, the electromomentum coupling actually depends on the time rate of the electric field rather on the electric field itself.

Numerical analysis
In this section, we apply the method that we have developed in Section 3 to numerically evaluate the effective properties-and specifically the electromomentum tensor-of exemplary piezoelectric composites. Specifically, we consider two compositions of PZT4 fibers within a PMMA matrix, the properties of which are given in Table 1. The phases are poled in the x 3 direction, such that in the (x 1 , x 2 ) plane they exhibit transverse isotropy. We consider a square lattice with the period a = 5 mm, and analyze two different cross-sections: one with a circular fiber of a radius R = 1.9 mm ( Fig. 1(a)), and one with a circular sector that is defined by the area πR 2 /8 ( Fig. 1(b)). Before we present our results, we note that our scheme was validated against several benchmark cases, as detailed in Appendix B. In Section 4.1, we show that the effective model defined by our homogenization scheme also exhibits the same dispersion relation as the composite in the case of free waves. As per our discussion in Section 1, this reinforcing feature on its own is not sufficient to validate the model, since the effective properties should also satisfy the mathematical restrictions imposed by causality, energy conservation and reciprocity, as we demonstrate in the sequel. By contrast, we will also show that homogenized models that ignore the electromomentum tensor violate these physical laws, even if these models recover the dispersion relation of the composite.

Free waves
We begin by analyzing the effective properties in the case of free waves, i.e., in the absence of sources, where we recall that the effective properties in such case are not uniquely defined. It is important to note that nonlocality induces components in the constitutive tensor that do not appear in local constitutive equations, even if the medium is isotropic, see Appendix C. Thus, to avoid information overload, we restrict attention to selected components of the effective tensors, and to the section -X − Γ − X of the 1st Brillouin zone ( Fig. 1(c)).
We begin by presenting in Fig. 2 the effective properties of the axisymmetric unit cell as functions of normalized frequencyω = ωa/2πc m , where c m is the velocity of shear wave in the PMMA matrix. We observe that all the standard properties-ρ,Ã,B andC-are real and satisfy the reciprocity requirement [32] ρ here and henceforth uppercase indices run from 1 to 6 according to Voigt notation. (For example, the components B 113 ,C 2323 andŜ † 331 are denotedB 15 ,C 44 andŜ 35 , respectively.) We examine next the coupling of Willis typeŜ andŴ (and their adjoints), the computation of which demonstrates the following features, supporting the validity of our homogenization scheme. (i) The tensorsŜ,Ŝ † ,Ŵ andŴ † are zero at κ = 0, and pure imaginary for κ ̸ = 0, in agreement with the understanding that the imaginary part captures nonlocal interactions, and that the real part should vanish when there is inversion symmetry 5 [33,36]; (ii) the tensors satisfy the symmetrieŝ as required from reciprocity and energy conservation [32]. From the calculation of the effective properties, it is possible to evaluate the dispersion relation that the homogenized medium exhibits. To do so, we substitute the effective constitutive relations [Eq. (26)] into the macroscopic governing equations [as mentioned, they are of the same form of Eq. (1), where the effective fields replace the microscopic ones], and extract the ω (κ) relation after some algebraic manipulation, as detailed in Appendix D. This resultant dispersion relation along the points Γ XMΓ is shown in Fig. 2(g), where we compare the corresponding acoustic and optical branches of the homogeneous medium (green circle marks), with the branches of the composite (black solid curves). Remarkably, the prediction of the homogenized medium recovers exactly the curves of the composite.
We recall that there is an ambiguity in the definition of the Bloch vector: if κ is a solution, then so is the sum of κ with any reciprocal lattice vector. To avoid this ambiguity, κ is restricted to the 1st Brillouin zone. Willis [74] questioned if this ambiguity extends to the definition of the effective properties, when later Srivastava and Nemat-Nasser [75] and Pernas-Salomón and Shmuel [57] showed that homogenizing only over the 1st Brillouin zone leads to nonphysical results. Specifically, Pernas-Salomón and Shmuel [57], who developed a similar method to calculate the effective properties of a periodic Euler-Bernoulli beam, showed that the effective properties using values of κ from the 1st Brillouin zone do not recover branches in the dispersion relation beyond the acoustic branch, 6 unless the values of κ are taken from ascending Brillouin zones. They further showed that these properties lead to unphysical reflections in scattering problems, see also Ref. [75]. This holds here as well: values of κ in the 1st Brillouin zone do not recover the optical branch of the original composite. Accordingly, to recover the second branch we have selected values of κ from the second Brillouin zone. reciprocity requirement [Eq. (27)]. In contrast, the computation of the couplings of Willis typeŜ,Ŵ and their adjoints provides different results than those of the axisymmetric cell, namely, here they do not vanish at κ = 0. Specifically, they are pure real in this limit, and become complex-valued for finite κ. Furthermore, we observe that their real part is even in κ, in consistency with reciprocity [32]. The evaluation of the dispersion relation using the homogenized properties is carried out in Fig. 3(b), and is compared with the dispersion relation of the original composite. Here again, the dispersion relation that the effective medium predicts recovers exactly the dispersion relation of the composite.
As mentioned, since there are no sources, the definition of the effective properties is not unique, and it is possible to define equivalent properties that neglect W (and W † ), while recovering the dispersion relation. However, we show next that equivalent properties without the electromomentum coupling are not physical, since they violate reciprocity and energy conservation. To do so, we note first that owing to the absence of sources, the independent fields are in fact algebraically dependent, as known from studies on standard (i.e., purely elastic) Willis materials [30]. Here, they are related via (Appendix D) for a generalized Willis medium that exhibits the piezoelectric and electromomentum effects. Accordingly, it is possible to absorbŴ into an equivalent dielectric tensor: that define the following equivalent constitutive equation for the electric displacement field: By construction, the equivalent description recovers the dispersion relation without usingŴ, however it violates the energy conservation and reciprocity conditions, 7Ã . This violation can be identified immediately from Eq. (30), by inspecting the symmetries that its comprising components (Ŵ,Ã andB) satisfy, and deducing that the diagonal terms ofÃ eq have a nonzero imaginary part that is odd κ. This is demonstrated Fig. 4, where A " eq (red curves) that is extracted for the composite with the asymmetric fibers is presented as function of the frequency of the dispersion curves along -X − Γ − X. For comparison, we also presentÃ ′′ (green circles); it is identically zero, as is should be.

Local approximation of the effective relations
The most prevalent use of metamaterials is in the long-wavelength regime, where the generated fields have a length much larger than the period of the composite. In such cases, the response of a metamaterial can be approximated as local, which implies that in our scheme we can approximate κ + G by G. We analyze next the local effective properties that result from this approximation. For brevity, we focus on the richer case of asymmetric fibers and present only selected components of the effective properties, where for completeness additional components are given in Appendix E. The selected properties as functions of normalized frequency are displayed in Fig. 5, and exhibit the following notable features. First, they respect causality, passivity and energy conservation, since they satisfy 8 [36,61] ∂ρ ∂ω ≥ 0, and they are all real. They exhibit a resonance at the normalized frequency ≈ 0.74, where beyond this frequency they change sign. Furthermore, we observe that the cross couplings of Willis type also respect causality, passivity and energy conservation since they satisfy [32] ∂B iJ ∂ω ≤ 0, In sharp contrast, a similar derivation for local equivalent properties that neglect the electromomentum coupling will result in a description that violates the physical requirements above. This is demonstrated in Fig. 6, where we show A eq as function of the normalized frequency. Specifically, in panel (a) we show A ′ eq , and observe that the components A ′ eq11 and A ′ eq12 have a negative slope; and in panel b we show A " eq , and observe it is nonzero with a negative slope. 9 Interestingly, a similar demonstration of the need for the inclusion of the electromomentum coupling in order to obtain a physical model was shown by Pernas-Salomón et al. [54]. There, the authors used a heuristic homogenization approach, based on the scattering response of an asymmetric piezoelectric element.

Effective properties outside the dispersion curve
We recall that in the presence of sources, the frequency and the effective wavevector are independent one of the other. We demonstrate next the calculations of the effective properties in such case, i.e., outside the dispersion curves. By way of example, in Fig. 7 we examine the case κ 1 = π/a, κ 2 = 0, which corresponds to a wavelength that is twice the period of the medium, and evaluate the properties as functions of the normalized frequency. We observe that the normalized resonance frequency has shifted to ≈ 0.48. We observe that the effective properties satisfy Eq. (32) and exceptŴ andŜ are all reals as they should be.

Summary
Rigorous homogenization theories are key in understanding the physics of heterogeneous media, and hence in the development of metamaterials. These theories have shown how composites can exhibit nonlocal effective properties that 8 We note that the derivatives in Eq. (32) should be evaluated at a fixed κ, as the partial derivates symbol implies. This requirement does not hold for the total derivatives, which correspond to the slopes in Figs. 2-3. (There, the properties are evaluated as functions of the frequency on the dispersion curves, where it is a function of κ.) 9 The requirement that in the local limit, the dielectric tensor of electromagnetic media satisfies ∂Ã ij /∂ω ≥ 0 can be found in Ref. [63]; the extension of this condition to generalized Willis media was derived in Ref. [32]. (e) Willis; and (f) electromomentum tensors as functions of the normalized frequency, calculated for PZT4 circular sector fibers in a PMMA matrix (Table 1).  (Table 1).
do not exist in their constituents [21,33,49,[59][60][61]. Here, we have developed explicit formulae to the homogenization theory of Pernas-Salomón and Shmuel [33], specialized for antiplane shear waves in fibrous piezoelectric composites. The formulae above are more accessible than the previous ones, since here we rely on the plane wave expansion method in the averaging process, instead of the Green's function method in Ref. [33], the knowledge of which is often unknown. Both formulations follow two principles that were advocated in rigorous electromagnetic and elastodynamics homogenization theories [21,31,36,49,[59][60][61]64]. First, the volume averaging over the microscopic Bloch-Floquet waves is carried out only over their periodic part. Second, the microscopic governing equations and constitutive relations include driving sources. The reason for considering them is mainly mathematical: while these sources enter our formulation in a form that is challenging to access experimentally, if at all possible, they allow us to derive a unique set of effective properties that is applicable for general frequency-wavevector pairs, not only for those that are related by the normal  (Table 1). (a) Mass density; and selected components of the (b) elastic; (c) dielectric; and (d) piezoelectric tensors. Panels e and f (resp. g and h) show the real and imaginary parts of the Willis tensor (resp. the electromomentum tensor). modes of the composite. Importantly, the resultant effective properties are physically meaningful, since they satisfy the mathematical restrictions that reciprocity, causality and energy conservation impose [32].
We have evaluated our formulae for two composites made of PZT4 inclusions in a PMMA matrix. The tensorial nature of the problem has allowed us to demonstrate certain symmetries of the effective properties that cannot be observed in one-dimensional problems. Our numerical examples have shown that the effective medium recovers the exact dispersion of free waves in the original composites. We have also calculated the approximated local effective properties, and the nonlocal effective properties at representative arbitrary frequency-wavevector pairs. Our computations demonstrate how the effective properties indeed respect reciprocity, causality and energy conservation.
When sources are not considered, i.e., for free waves, it is possible to define alternative homogenized descriptions that also recover the dispersion relations of the original composite. In spite of the fact that such descriptions produce this particular measurable quantity properly, they are not admissible, since they generally predict nonreciprocal, nonconservative response for reciprocal, conservative media. Such violations were discussed in detail in Refs. [36,61,64], who studied analogous electromagnetic and acoustic problems. Lending the terminology in these works, we referred to the alternative models that neglect cross-couplings of Willis type as equivalent models. As explained in Refs. [36,61,76], while the secondary properties (impedance and wavevector) of the effective and equivalent sets coincide, the primary properties of the equivalent set lack physical meaning. We have evaluated the properties of the equivalent model for our case study, and showed that indeed they (more specifically, the equivalent dielectric tensor) violate reciprocity, causality and energy conservation.
It is imperative to examine the implications of these homogenization results-derived for an infinite medium-in predicting the scattering response of a finite piezoelectric medium with asymmetric microstructure (cf. the analysis and applications in Refs. [43,57,75,77,78] for acoustic and elastic problems). In one-dimensional scattering problems, Pernas-Salomón et al. [54] showed that the resultant electromomentum coupling is related to directional phase angle that depends on the electric circuit conditions; such a study in two-dimensional settings as considered here is left for future work.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
The description of matrices is given here where the components of the assembly matrix L A are

Appendix B. Validation of the homogenization scheme
In addition to the validating features listed in Section 4, the validity of our effective model is further supported by its predictions in the low frequency, long-wavelength limit. In this limit, the effective mass density of our model is equal to the volume average, as it should [79]. Furthermore, in this limit, there are also several benchmark results in the onedimensional setting. To test our model against these results, we used our two-dimensional homogenization scheme to analyze a layered unit cell whose properties are uniform along x 2 and periodic along x 1 . For this microstructure, we considered a wavevector whose x 2 -component is null and its x 1 -component is infinitely small, but greater than 0. By doing so, we reduce our problem to a one-dimensional problem in the long-wavelength limit of propagation in the periodicity direction. When taking B (x) − → 0, we find thatC 55 (:=C 1313 ) andÃ 11 are equal to their weighted harmonic mean, agreeing with their expected values from static homogenization in the limiting decoupled case. When taking nonzero B (x), we find numerically that Result (B.1) was obtained analytically by Pernas-Salomón et al. [54] using heuristic homogenization of one-dimensional piezoelectric scatterers. Pernas-Salomón et al. [54] also obtained that (i) the local Willis modulus depends on the sum of the asymmetry in the elastic impedance Z 2 E = ρC, as well-known in the purely elastic case [36,44], and the asymmetry in a piezoelectric-like impedance, ρB; ((ii) the local electromomentum modulus depends on the asymmetry in the ratio B/A.
Using a one-dimensional version of the scheme developed here, we have computed the effective properties of a trilayer whose constituents are such ρB, ρC and B/A exhibit inversion symmetry; we received zeroS 53 (:=S 133 ) and zeroW 13 .
When breaking the symmetry of B/A in the x 1 -direction, we received nonzeroW 13 . Subsequently, we examined a case were the inversion symmetry of ρB is broken; indeed, this case resulted with nonzeroS 53 . The agreement in the above benchmarks between the two homogenization approaches reinforces their validity. We numerically observed, however, that in our periodic homogenization, if the mass density is spatially constant, then both the Willis coupling and the electromomentum coupling are null; in the heuristic homogenization, this condition is not necessary.
Finally, we analyzed a second case, by considering a wavevector whose x 1 -component is null and its x 2 -component is infinitely small, but greater than 0. By doing so, we reduce our problem to a one-dimensional problem in the longwavelength limit of propagation along the layers. When taking B(x) − → 0, we find thatC 44 (:=C 2323 ) andÃ 22 are equal to their volume average, agreeing with their expected values from static homogenization in the limiting decoupled case.

Appendix C. On nonlocal components in constitutive equations
As mentioned in the main text, nonlocality induces components in the effective tensor that do not appear in local constitutive equations, even if both media are isotropic. (Such components are somewhat overlooked in the literature, but they do exist.) For example, the constitutive relation of a statistically homogeneous nonlocal linear isotropic elastic solid is [80] σ ij (x) = where ϵ is the strain and for some material moduli λ, µ, α and β that are functions of the translation ⏐ ⏐ x − x ′ ⏐ ⏐ . The moduli α and β multiply terms that vanish in local isotropic solids. Thus, by way of example, the component C 1234 must be null in local isotropic solids, while having a finite value for κ > 0. Analogous results hold for electromagnetic constitutive relations, see, e.g., Eq.
(D.1)  Table 1 for the local approximation (κ = 0) case, as function of the normalized frequency.
Our goal is to find the equation that we will use to recover the dispersion relation for any prescribed κ. To this end, We use Eq. (D.2) to obtain iκ T (C iκ ⟨u⟩ +B † iκ ⟨φ⟩ −S ⟨iωu⟩ ) = −iω (S † iκ ⟨u⟩ +W † iκ ⟨φ⟩ −ρ ⟨iωu⟩ ) . (D.7) We rearrange the equation to get quadratic equation that will allow us to get the dispersion relation where a, b and c are (D.10) Using Eq. (D.10) we can obtain the eigenfrequencies that corresponds to the effective properties. The first two branches that are fully recovered are shown in Fig. 3. The dispersion relation of the composite was calculated from the eigenvalue problem det Appendix E. Additional components of the effective properties for the asymmetric cell Fig. 8 presents additional effective properties for the asymmetric cell in the local approximation. Fig. 9 presents additional effective properties for the asymmetric cell when the 1-component of wavevector is κ 1 = π/a and 2-component is null.