Materials with asymmetric microstructure can constitutively couple macroscopic ﬁelds from dif- ferent physics. Examples include piezoelectric materials that couple mechanical and electric ﬁelds and Willis materials that anomalously couple dynamic and elastic ﬁelds. Recently, it was shown anomalous coupling between the elastodynamic and electric ﬁeld emerges when piezoelectricity is incorporated into Willis materials. Here, we investigate such one-dimensional piezoelectric Willis elements using heuristic homogenization, long-wavelength asymptotic analysis and numerical exper- iments. We show that an anomalous electromomentum coupling is necessary to explain the scattering properties of the asymmetric piezoelectric media by a macroscopic description that respects reci- procity and energy conservation. Our ﬁndings elucidate the origins of the electromomentum coupling and provide insight for the future design of this new class of coupled-ﬁeld metamaterials. 1 Introduction Metamaterials are engineered composites materials with extraordinary effective properties that offer func- tionalities beyond those of standard materials [9, 16, 46, 47, 50]. Prominent examples in elastodynamics and acoustics include metamaterials with negative mass density and negative stiffness, which are in sharp contrast with the positive moduli of their constituents [6, 11, 12, 19]. A separate class of metamateri- als, known as Willis materials, exhibit unique properties on the macroscale that their constituents do not posses, namely moduli that couple the momentum and strain ﬁelds as well as the velocity to stress ﬁelds [20, 21, 23, 30, 41]. Willis materials therefore offer additional degrees of freedom to manipulate me- chanical waves similar to the way that bianisotropic materials allow greater control over electromagnetic waves [1, 4, 24, 25, 29, 45]. *Corresponding author: meshmuel@technion.ac.il 1 The emergence of these non-intuitive effective moduli is the result of a homogenization scheme for elastic composites that was initially developed by Willis [51, 52, 54, 55, 56] to understand the elastody- namic response of heterogeneous elastic media, and revisited by many others in the context of metamate- rial research [10, 24, 28, 29, 32–35, 37, 40, 43, 49]. Of relevance to this work are the conclusions of Sieck et al. [45] and Muhlestein et al. [30], who demonstrated that the constitutive relations of Willis form are required to describe a physically admissible macroscopic response, i.e., one that respects the principles of causality and passivity. This necessity is analogous to the need to include bianisotropic coupling to obtain a physically valid macroscopic description of the dynamic behavior of electromagnetic composites [2, 3]. By generalizing the approach of Willis [54] to account for piezoelectric constituent materials, Pernas- Salomón and Shmuel [38] have found that the dynamic response of such materials is described with con- stitutive relations that include effective moduli that couple the velocity and momentum of the medium to the electric ﬁelds. These generalized Willis couplings are therefore referred to as electromomentum couplings. From the viewpoint of applications, the signiﬁcance of this result is twofold. First, it extends the design space of elastodynamic metamaterials beyond Willis materials to include coupling to the elec- tromagnetic domain. Second, it provides a new avenue for mechanical wave manipulation using external electric sources. The homogenization procedures of Willis [52, 53, 54], and in turn Pernas-Salomón and Shmuel [38], provide effective constitutive relations that are nonlocal in time and space. These relations become local in space when the wavelength is much larger than the microstructure, or when they are used to describe a single subwavelength element rather than a bulk, in which case they are termed the Milton-Briane-Willis equations [25–27]. Analysis of the effective parameters shows that the local part of the Willis coupling coefﬁcients originates from broken inversion asymmetry [29, 45]. Recent experiments that demonstrate the local contributions to the Willis coupling were carried out in these settings, i.e., using asymmetric subwavelength structures [7, 8, 17, 20, 21, 30]. The objective of the present work is to provide an analogous evidence of a local electromomentum coupling by investigating the scattering response of an asymmetric piezoelectric element embedded in an elastic background material. We consider Gedanken experiments where a heterogeneous piezoelectric element separates two ho- mogeneous elastic slabs and scatter one-dimensional strain waves incident from the slabs. We calculate analytically the experimentally measurable reﬂection and transmission coefﬁcients, and develop a re- trieval method to extract the effective properties from these coefﬁcients by adapting the schemes in Refs. [13, 30] to piezoelectric media. The heuristic homogenization method applied here is applicable only in the long-wavelength regime [46]. The approach is employed in conjunction with two different homoge- nized models: one that includes the electromomentum coupling, and one that neglects it. Borrowing the terminology in Refs. [3, 45], we refer to the former model as the effective model and to the latter model as the equivalent model. 2 We ﬁnd that the properties that are retrieved using the equivalent model depend on the excitation setup, and violate reciprocity and energy conservation [3, 45], even in the long-wavelength, low fre- quency limit. By contrast, we show that the properties that are retrieved using the effective model satisfy these physical laws and are independent of the excitation setup. This discrepancy of the equivalent model stems from its inability to capture the following: the reﬂection of incident plane waves has a directional phase angle that depends on the electric ﬁeld . The electromomentum coupling captures precisely this phenomenon, and the inclusion of this coupling in the homogenized model renders it physically admissi- ble. This demonstration is analogous to the analysis of Muhlestein et al. [30] for acoustic media where the authors demonstrated, both analytically and experimentally, that scattering from an asymmetric acoustic element can only be represented with physically meaningful material properties when using the effective medium description. We derive explicit low frequency expressions for effective properties using Taylor series. While the leading term in the expansions of the standard properties is a real constant, the leading term in the expansions of the Willis coupling and the electromomentum coupling is a pure imaginary linear function of frequency and purely imaginary. The expansions further show that the couplings depend on the order of the composition, and that certain compositions exhibit zero Willis coupling and ﬁnite electromomentum coupling. The directional phase angle of such elements vanishes in the short circuit conﬁguration and emerges in open circuit conditions, a feature that can be harnessed to control the directional phase change by the opening and shortening of the circuit. The analysis is presented as follows. First, we provide a short summary of the relevant governing equations in Sec. 2. Section 3 provides a simple analysis of plane waves in piezoelectric media with and without the electromomentum coupling, through which we highlight its effect on the characteristic elastic impedance and phase velocity. The conceptual experiments that demonstrate this effect are described in Sec. 4, together with the calculations of the resultant scattering and retrieved properties of the two models. Section 5 provides the low frequency expansions of the effective properties. Numerical examples are then presented in Sec. 6, and the work concludes with a summary of the work with ﬁnal remarks in Sec. 7. 2 Field equations and constitutive relations Pernas-Salomón and Shmuel [38] developed a comprehensive dynamic homogenization theory for piezo- electric composites. Here we review the equations relevant to the current work in the local, one-dimensional, source-free settings1. In the absence of body forces and free charge, a one-dimensional piezoelectric in- homogeneous medium satisﬁes the balance equations σ,x−p,t= 0, D,x= 0, (1) 1For complete details on the theory in the general settings, see Refs. [38, 39]. 3 for the stress σ, momentum density p, and electric displacement ﬁeld D, wherexis the spatial coordinate. Faraday’s equation for the electric ﬁeld Eis satisﬁed by setting E=−φ,x, whereφis the electric potential. At each point, the ﬁelds of Eq. (1) are related to E(orφ,x) and the derivatives of the axial displacement uthrough the constitutive equations [5]  σ D p = C B 0 B−A0 0 0ρ  u,x φ,x u,t , (2) whereρ,A,B, andC, are the mass density, and dielectric, piezoelectric, and elastic moduli, respectively2, which vary in space since the medium is inhomogeneous. We clarify that we restrict attention to a reciprocal and lossless medium, such that the matrix in Eq. (2) is real, symmetric and independent of the boundary conditions. This form of the constitutive relations is referred to as the stress-charge form, whereAandCare measured under the conditions of constant strain and electric ﬁeld, respectively. With the objective of replacing the inhomogeneous medium by a homogeneous one that macroscop- ically behaves the same (in some appropriate sense3), Pernas-Salomón and Shmuel [38] obtained the following constitutive equations for the effective ﬁelds  σ D p = ˜C ˜B ˜S ˜B−˜A˜W ˜S ˜W ˜ρ  u,x φ,x u,t , (3) where we have used the same notation for the effective ﬁelds and the microscopic ﬁelds, while the ef- fective properties are distinguished by overhead tilde. Notably ˜Sand ˜Ware the emergent Willis and electromomentum coupling constants, respectively, which appear only in the macroscopic description. We recall that Eq. (3) is the reduced form of the general relations found in Ref. [38], when considering one-dimensional media free of sources, in the long-wavelength limit. We note that the homogenization process that lead to Eq. (3) uses ensemble averaging to deﬁne the effective ﬁelds [54, 55]; in the long- wavelength limit, this averaging can be approximated by volume averaging. Using simple conceptual experiments, we will demonstrate that one must account for ˜Win order to obtain a physically admissi- ble effective description of piezoelectric scatterers under dynamic loading scenarios. Such description should also respect energy conservation and reciprocity, and be independent of the boundary conditions, as the heterogeneous medium. 2The dielectric and piezoelectric moduli are often denoted in the literature by εandeinstead of the notation AandBused here [5, 42]. 3See Refs. [32, 54, 55] for details. 4 3 Plane wave analysis and the main idea Before we begin with the analysis of the scattering experiments, it is instructive to analyze one-dimensional plane waves in piezoelectric Willis media, through which we can highlight the role of the cross-couplings on the elastic impedance and phase velocity. To this end, we consider a homogenized medium that is gov- erned by Eq. (3), and examine the equation of motion after we substitute the constitutive relations for σ andpto yield ˜Cu,xx+˜Bφ,xx−˜Wφ,xt−˜ρu,tt= 0. (4) Note that ˜Su,txand˜Su,xtcancel out, in agreement with the observations of Muhlestein et al. [31] for the case when the Willis tensor is aligned with the direction of propagation4As is the case for standard piezoelectric media, further derivations depend on the form of the electric ﬁelds. We distinguish between the two cases: ( i) short circuit for which E= 0, and ( ii) open circuit for which D= 0. The piezoelectric analysis here extends that of Refs. [30, 31, 45] for purely acoustic Willis materials, and is related to the analysis by Kutsenko et al. [18] on the Bloch impedance of periodic piezoelectric media. 3.1 Plane waves in a short circuit conﬁguration The condition E= 0≡−φ,xsimpliﬁes Eq. (4) to the standard one-dimensional wave equation for an elastic medium with plane wave modulus ˜Cand density ˜ρ, ˜Cu,xx−˜ρu,tt= 0, (5) which is satisﬁed by plane wave solutions of the form u(x,t) =Uei(±˜kEx−ωt). This deﬁnes the phase velocity ˜v2 E:=˜C/˜ρand the wavenumber ˜kE=ω/˜cEin the absence of electric ﬁeld. In terms of these quantities, we write the characteristic elastic impedance Z=σ/u,tin a short circuit conﬁguration. It immediately follows from the substitution of the plane wave form into the constitutive equation for σthat Z±=±˜ZE+˜S, where ˜Z2 E:= ˜ρ˜C, and the sign of the superscript for Zdesignates the wave direction. The implication is that when ˜Sexists, the impedance is asymmetric with respect to the propagation direction. This motivates us to follow Muhlestein et al. [30, 31] and deﬁne the asymmetry factor in the absence of electric ﬁeld, χE=˜S/˜ZE, in order to rewrite the impedance as Z±= (±1 +χE)˜ZE. (6) Equation (6) shows that the asymmetry factor χEin the absence of an E-ﬁeld is the same as the asymme- 4The contribution of the Willis coupling to the equations is retained in any other case where its polarization is not parallel or anti-parallel to the propagation direction [Eq. (2.8) therein]. 5 try factor in Willis media without the electromomentum coupling. It has been established that the local Willis coupling of a reciprocal and passive medium is purely imaginary, and its leading term in a low frequency expansion is linear [22, 29, 33]. Accordingly, we can relate the asymmetry of Zto a purely real coupling coefﬁcient5ˆSthrough the relation ˜S=iωˆS, that is asymptotically tractable, in terms of which we rewrite the impedance as [30] Z±=±˜ZE+iωˆS. (7) Eq. (7) highlights the fact that the characteristic elastic impedance is complex and changes phase angle with direction in media with local Willis coupling. Muhlestein et al. [30] linked this impedance to the scattering behavior of an acoustic asymmetric element, and related its imaginary part to the asymmetry in the reﬂection coefﬁcients when plane waves are incident from opposing directions. This asymmetric response led to their observation that one must include Willis coupling in the dynamic effective properties to obtain a physically valid effective description of the asymmetric element, since ˜ZEmust be real and direction independent in order for ˜Cand˜ρto satisfy reciprocity and energy conservation. The remaining insights needed to understand the electromomentum coupling are obtained from the comparison with the open circuit conﬁguration, D= 0, considered next. 3.2 Plane waves in an open circuit conﬁguration It follows from Gauss’s law [Eq. (1) 2] that the electric displacement ﬁeld is zero in the absence of free charge when the electrical components are in open circuit boundary conditions. We can then use the derivatives of Eq. (1) 2and the constitutive relation for Dto obtain expressions for φ,xx, andφ,xt, to rewrite the equation of motion (4) as /parenleftBigg ˜C+˜B2 ˜A/parenrightBigg u,xx−/parenleftBigg ˜ρ+˜W2 ˜A/parenrightBigg u,tt= 0. (8) Solutions to Eq. (8) are plane waves of the form u(x,t) =Uei(±˜kDx−ωt)with the wavenumber ˜kD:= ω/˜vDand the phase velocity ˜v2 D:=˜CD/˜ρD, where ˜CD:=˜C+˜B2/˜Aand˜ρD:= ˜ρ+˜W2/˜A. This simple analysis reveals an appealing analogy between the piezoelectric effect and the electromomentum effect: the latter alters the apparent mass density in the same way that the former stiffens the apparent elastic modulus [5]. Therefore, evidence for ˜Wcan be found in a dynamic open circuit experiment that will prove that the effective phase velocity is given by ˜CD/˜ρD, not ˜CD/˜ρ. More insights are obtained from the analysis of the characteristic impedance and asymmetry factor. To this end, we substitute the above form foruinto the constitutive equation for σ, eliminateφ,xusing the constitutive equation for D= 0, 5This actually motivates a different form for the constitutive relations, in which the stress depends on the acceleration and the momentum depends on the strain rate, see Refs. [26, 29, 39, 45]. 6 and obtain σ=i˜kD˜Cu+i˜kD˜B2 ˜Au−iω˜W˜B ˜Au−iω˜Su. (9) Using the deﬁnition of ˜kDwe ﬁnd that the characteristic plane wave impedance in open circuit is Z±:=σ −iωu= (±1 +χD)˜ZD,˜Z2 D:= ˜ρD˜CD, χ D:=1 ZD/parenleftBigg ˜S+˜B ˜A˜W/parenrightBigg , (10) where ˜ZDmust be real for the constitutive properties to respect reciprocity and energy conservation. We again observe that the characteristic impedance is asymmetric in the propagation direction, however the asymmetry factor for the open circuit condition, χD, depends on both the conventional Willis coefﬁcient andthe elecromomentum coupling ˜W. Previous works have shown that ˜W—like ˜S—is purely imaginary in reciprocal and conservative local media, and its leading term in a low frequency expansion is a linear function of frequency [38, 39]. Therefore, we can relate the asymmetry in Zwith a purely real coupling coefﬁcient6that is asymptotically tractable, say ˆW, through the relation ˜W=iωˆW, in terms of which we rewrite the impedance as Z±=±˜ZD+iω/parenleftBigg ˆS+˜B ˜AˆW/parenrightBigg . (11) Equation (11) indicates that the electromomentum coefﬁcient alters the directional dependence by alter- ing the phase angle in addition to the ˆScontribution. This suggests a simple test for the existence of ˜W: Examine the directional phase change of waves that are reﬂected from an asymmetric piezoelectric ele- ment in scattering experiments, both in short and open circuit conﬁgurations; if they do not match, it will provide evidence that ˜Wis non-zero and must be considered in order to determine physically meaningful effective properties of the medium. This is precisely what we demonstrate in the sections that follow by extending the scheme of Muhlestein et al. [30] to piezoelectric media. 4 Scattering from a piezoelectric Willis element Consider a piezoelectric stack connecting two semi-inﬁnite slabs made of an elastic isotropic material with plane wave modulus C0and densityρ0. The piezoelectric stack comprises three layers, 1,2and3, of lengthsl1,l2andl3, respectively, as shown in Fig. 1. Thin electrodes are located at the interfaces between the stack and the elastic slabs. Having setting up the Gedanken experimental apparatus, we conduct a scattering analysis, from which we analytically determine the transmission and reﬂection coefﬁcients of 6This motivates a different form for the constitutive relations, in which the electric displacement ﬁeld depends on the acceleration and the momentum depends on the time rate of the electric ﬁeld, see Ref. [39]. 7 Incident waveTransmitted waveReflected wavex0x1x2x3ρ0,C0ρ0,C0Piezoelectric element Reflected wavex0x1x2x3Short circuit ρ0,C0ρ0,C0Piezoelectric elementIncident waveTransmitted waveOpen circuit FIG. 1: Schematics of the conceptual experiments on a piezoelectric element covered by two electrodes at x0and x3. The element is made of three layers which connect two semi-inﬁnite elastic slabs made of an isotropic material with plane wave modulus C0and densityρ0. When the switch is open (top ﬁgure) and closed (bottom ﬁgure), the setup corresponds to open- and short circuit electric boundary conditions, respectively. For each setup we analyze an experiment where the incident waves are incoming from the left waveguide (top ﬁgure), and an experiment where the incident waves are incoming from the right waveguide (bottom ﬁgure). The phase of the reﬂected waves depends on the direction of the incoming waves and the electric circuit conditions. This dependency is captured by the electromomentum coupling. the piezoelectric stack. We will use these coefﬁcients in an inverse program to determine the effective properties of a ﬁctitious homogenized medium that reproduces the same scattering that the heterogeneous stack generates [13, 30, 48]. The retrieved effective properties of the homogenized medium depend on the model chosen for its constitutive relations. It is important to note that while the analysis that follows considered a single, heterogeneous scatterer, which can also be understood as a unit cell from the perspective of polarizability [2, 21, 41, 45, 48], the approach employed here is representative of experimental methods used to characterized material properties and can be extended to consider the response of many unit cells [30]. In this section we describe a scheme to retrieve the effective properties based on two homogenized models: one that includes the electromomentum coupling [Eq. (3)], and one equivalent model that neglects it. We will show that while the latter model delivers (by construction) the same scattering coefﬁcients as our effective model (hence the term equivalent ), the resulting properties violate reciprocity and energy conservation and therefore lack physical meaning. Further, the properties retrieved using the equivalent model depend on the excitation setup, even in the low frequency limit. 4.1 Retrieved properties from short circuit conditions Consider a setup where the two electrodes at the ends of the piezoelectric stack are connected in a short circuit conﬁguration (switch is closed in Fig. 1), such that the average (or effective) electric ﬁeld is null. 8 Two experiments are considered in this conﬁguration: in experiment L (resp. R), the stack is excited a time-harmonic longitudinal wave that is incoming from the left slab (resp. right). The incident wave is partly reﬂected at the interface with the stack and partly transmitted to the opposite slab. We assume a time dependency e−iωtand write the resultant displacement ﬁeld in the elastic slabs as u(x,t) =  /parenleftBig c+ Leik0x+c− Le−ik0x/parenrightBig e−iωt, x≤x0, /parenleftBig c+ Reik0x+c− Re−ik0x/parenrightBig e−iωt, x≥x3,k0=ω/radicalBiggρ0 C0(12) wherec+ L(resp.c− R) is the prescribed amplitude of the incident wave in experiment L (resp. R), and c− R(resp.c+ L) is zero. These experimentally measurable coefﬁcients can be calculated analytically by introducing the state vector h(x)7that is continuous in x, h(x) = (u,φ,σ,D )T. (13) The values of h(x)at the endsxn−1andxnof thenthpiezoelectric layer are related via [18] h(xn) =Tem nh(xn−1),Tem n= coskDnln 0 sinkDnln/ωZDnBn AnsinkDnln/ωZDn Bn An(coskDnln−1) 1Bn AnsinkDnln/ωZDn/parenleftBig Bn An/parenrightBig2sinkDnln/ωZDn−ln An −ωZDnsinkDnln0 cos kDnlnBn An(coskDnln−1) 0 0 0 1 , (14) wherekDn=ω/radicalBig ρn/CDn,ZDn=√ρnCDnandCDn= (Cn+B2 n/An). Continuity of h(x)implies that h(x3) =TLh(x0),TL=Tem 3Tem 2Tem 1, (15) in experiment L, while in experiment R, the composition order of Tem nis reversed. We focus ﬁrst on experiment L. In this case, it follows from Eq. (15) that φ(x3)−φ(x0) =TL 21u(x0) +TL 23σ(x0) +TL 24D, (16) where TL ijis theij-component of the matrix TL, and we have used the fact that electric displacement is constant inside the stack. Since the electrodes are in short circuit such that the potential difference is 7Henceforth, with an abuse of notation, we omit the dependency on t, and recall that it is always in the form e−iωt. 9 zero, i.e.,φ(x3)−φ(x0) = 0 , then Eq. (16) together with Eq. (15) implies  u(x3) σ(x3) =Ts(x3,x0) u(x0) σ(x0) ,Ts= TL 11−TL 14 TL 24TL 21TL 13−TL 14 TL 24TL 23 TL 31−TL 34 TL 24TL 21TL 33−TL 34 TL 24TL 23 . (17) Equation (17) relates the displacement and stress at the two ends of the stack, and hence the amplitudes of the incoming and outgoing waves. By applying the continuity conditions for the displacement and stress ﬁelds at x0andx3and using Eqs. (12) and (17), we relate the amplitudes of the waves through the following coefﬁcients transfer matrix  c+ R c− R =Ks(R,L) c+ L c− L , (18) where Ks(R,L) = Q0(x3)−1Ts(x3,x0)Q0(x0),Q0(x) = eik0xe−ik0x iC0k0eik0x−iC0k0e−ik0x . (19) Equation (18) returns the output of experiment L when c− R= 0with the reﬂection and transmission am- plitudesrL=c− L/c+ LandtL=c+ R/c+ L, respectively. Similarly, Eq. (18) returns the output of experiment R usingc+ L= 0with the reﬂection and transmission amplitudes rR=c+ R/c− RandtR=c− L/c− R, respectively. The matrix Kscan be expressed in terms of these amplitudes [36], namely, Ks(R,L) = tL−rRt−1 RrLrRt−1 R −rLt−1 Rt−1 R . (20) Equations (17)-(20) together imply the L-experimental results tL= 2e−ik0l/parenleftBig Ts 22+Ts 11−iωZ 0Ts 12−(iωZ 0)−1Ts 21/parenrightBig−1, (21a) rL=/parenleftBig Ts 22−Ts 11−iωZ 0Ts 12+ (iωZ 0)−1Ts 21/parenrightBigtL 2. (21b) The analysis for experiment R mirrors the above procedure, showing that rRandtRequal, respectively, the right hand sides of Eqs. (21b) and (21a), after replacing Ts ij/parenleftBig TL/parenrightBig byTs ij/parenleftBig TR/parenrightBig . Notably, we ﬁnd tL≡tR, as expected for a reciprocal medium, while the phase of rLis may differ from the phase of rR. This difference depends on the spatial asymmetry of the stack, and is a manifestation of the directional dependency of the effective characteristic impedance, as discussed in Sec. 3.1. The objective now is to determine the properties of a ﬁctitious homogeneous medium that exhibits the 10 same scattering coefﬁcients as the stack, i.e., will reproduce the same Ks(R,L)—under the same bound- ary conditions. As per this short circuit conﬁguration, the electric ﬁeld in the ﬁctitious homogeneous medium should vanishes, like the average electric ﬁeld in the stack. The constitutive relations for σand pof both the effective and equivalent models reduce to σ=˜Cu,x−iω˜Su, (22a) p=˜S†u,x−iω˜ρu. (22b) (Later on, we will use superscript eqto distinguish the retrieved properties of the equivalent model that neglects ˜Wfrom the retrieved properties of the effective that includes ˜W). In this case, the constitutive equations for σandp[Eq. (22)] involve only the effective mechanical properties and thus reduce to the Willis form, in agreement with the analysis provided in Sec. 3.1. We note that unlike in Eq. (3), Eq. (22b) does not assume a priori that the conditions imposed by reciprocity are satisﬁed [29, 39], which would imply that the coupling between σandu,tis equal to the coupling between pandu,x. To maintain the additional generality, we denote by ˜S†the coefﬁcient that couplespwithu,x, and show by analysis that ˜S=˜S†. Before we proceed to determine ˜C,˜ρ,˜Sand˜S†, we emphasize that the remaining constitutive equation for the electric displacement ﬁeld depends on the form of the chosen model, i.e., whether or not it accounts for the electromomentum coupling. This will be addressed in Sec. 4.2. To determine the properties ˜C,˜ρ,˜Sand˜S†of a ﬁctitious homogenized medium with the same scat- tering response as the stack, we relate the displacement and stress ﬁelds at the ends of the homogenized medium via (Appendix A) u(x3) σ(x3) =Tm(x3,x0) u(x0) σ(x0) , (23) where Tm(x3,x0) =elMis the transfer matrix of the homogenized medium, deﬁned by M= iω˜S˜C−1 ˜C−1 ω2˜S˜S†˜C−1−ω2˜ρ−iω˜S†˜C−1 . (24) We require that the homogenized medium reproduces the same scattering properties, i.e., the same Ks. This implies that the variables in Eq. (24) are determined from [see Eqs. (19)-(20)] Q−1 0(x3)TmQ0(x0) = tL−rRt−1 RrLrRt−1 R −rLt−1 Rt−1 R , (25) which yield equations for Tmin terms of the reﬂection and transmission coefﬁcients rL,rR,tLandtR. 11 RecallingtL≡tR=:t, simpliﬁes Tmto Tm=1 2t ϑeik0l+e−ik0l+rR−rL1 iωZ0/bracketleftBig ϑeik0l−e−ik0l−rR−rL/bracketrightBig iωZ 0/bracketleftBig ϑeik0l−e−ik0l+rR+rL/bracketrightBig ϑeik0l+e−ik0l−rR+rL , (26) whereϑ=t2−rLrR, andZ0=√ρ0C0is the characteristic impedance of the elastic slabs. By applying Sylvester’s formula to the logarithm of Tm, Eq. (26) admits the following simpliﬁcation [14, 43, 44] M=ψ lsinψ/parenleftbigg Tm−1 2ItrTm/parenrightbigg ,cosψ=1 2trTm=ϑeik0l+e−ik0l 2t, (27) wheree±iψare the eigenvalues of Tm. We deduce from Eq. (27) that M22=−M11, and the effective elastic properties are ˜S=M11 iωM 12=Z0(rR−rL) (t2−rLrR)eik0l−e−ik0l−rR−rL, (28a) ˜S†=−M22 iωM 12≡˜S, (28b) ˜C=M−1 12=2iωZ 0ltψ−1sinψ (t2−rLrR)eik0l−e−ik0l−rR−rL, (28c) ˜ρ=detM ω2M12=−iZ0ψ 2tωlsinψ /parenleftBig (t2−rLrR)eik0l−e−ik0l/parenrightBig2−4rRrL (t2−rLrR)eik0l−e−ik0l−rR−rL , (28d) where, as mentioned, in this setup the elastic equivalent properties are equal to the elastic effective prop- erties. Equation (28) expresses the effective properties directly in terms of the experimentally measurable scattering coefﬁcients [which were analytically calculated using Eq. (21)]. Thus, by construction, the re- trieved models deliver the same scattering response that the piezoelectric stack produces. The above experiments are consistent with the conditions that are required to evaluate the properties of a standard piezoelectric material in the stress-charge form, in the sense that the stiffness was evaluated at a pre- scribed (effective) electric ﬁeld. To summarize the results of our Gedanken experiments so far: we ﬁnd frequency-dependent effective properties, even though the assumed properties of each constituent are frequency-independent. Notably, the Willis couplings ˜Sand˜S†are equal one to another and depend on the spatial asymmetry of the stack throughrR−rL. If the stack has certain symmetry properties, then the reﬂection of incoming waves from the right is the same as from the left and the Willis couplings vanish, as we will explicitly show later using low frequency approximations for the dynamic effective properties. 12 4.1.1 First inconsistency in retrieving the piezoelectric coefﬁcient when neglecting ˜W A different perspective on the retrieval process is to consider the retrieved properties as those that relate the effective ﬁelds in the piezoelectric element. In the low frequency, long-wavelength limit, the effective ﬁelds are deﬁned as the volume average of the microscopic ﬁelds. Evaluating the volume averages over the constitutive relations for σandpdeliver equations for ˜ρ,˜Cand˜Sthat are consistent with Eq. (28a). The third equation for D, namely, /angbracketleftD/angbracketright=˜B/angbracketleftu,x/angbracketright−iω˜W/angbracketleftu/angbracketright, (29) is not sufﬁcient to determine ˜Band˜W. If we ignore the electromomentum effect, i.e., we assume that the effective electric displacement is independent of the effective velocity ﬁeld, then this experiment is sufﬁcient to determine the equivalent piezoelectric coefﬁcient ˜Beqfrom the average ﬁelds. Speciﬁcally, an additional equation ˜Beqis obtained from the requirement that the equivalent properties of the homog- enized model deliver the same average ﬁelds8that are generated in the piezoelectric stack. Thus, when ˜Wis neglected in Eq. (3), the equivalent piezoelectric coefﬁcient is determined from the relationship /angbracketleftD/angbracketright=˜Beq/angbracketleftu,x/angbracketright,so˜Beq:=/angbracketleftD/angbracketright /angbracketleftu,x/angbracketright, (30) where/angbracketleft/angbracketrightdenotes volume averaging over x0≤x≤x3. Now, if ˜Beqis to satisfy energy conservation and reciprocity, in the long-wavelength limit it should be a real-valued quantity that is independent of the direction of the wave. Stating explicitly the implication of the latter condition to our analysis, it requires that ˜Beq L:=/angbracketleftD/angbracketrightL//angbracketleftu,x/angbracketrightLwill be equal to ˜Beq R:=/angbracketleftD/angbracketrightR//angbracketleftu,x/angbracketrightR, where we recall that subscripts L and R correspond to ﬁelds that are generated in the left ( c+ L= 1,c− R= 0) and right ( c+ L= 0,c− R= 1) experiments. However, using Eq. (16) to determine /angbracketleftD/angbracketrightαand Eq. (12) to express /angbracketleftu,x/angbracketrightαin terms of the scattering coefﬁcients, we ﬁnd that ˜Beq α=l/bracketleftBig Tα 21+iωZ 0Tα 23+rαeik0l(Tα 21−iωZ 0Tα 23)/bracketrightBig (1 + (rα−t)eik0l)Tα 24, α= L,R, (31) implying that ˜Beq L/negationslash=˜Beq Rand moreover Im[˜Beq α]/negationslash= 0. These nonphysical artifacts are generated when neglecting the electromomentum effect and persist in the low frequency limit, as we will demonstrate in Sec. 6. Before we proceed to analyze the open circuit conﬁguration, we summarize the results from the short circuit conﬁguration. We have calculated the reﬂection and transmission coefﬁcients for incidence from both sides of the heterogeneous scatterer. In terms of these coefﬁcients we have estimated the effective 8Again, we recall that the use of volume averages instead of ensemble averages is valid only in the long-wavelength, low frequency regime. 13 short circuit dynamic stiffness, ˜C, the short circuit dynamic density ˜ρ, and the Willis coupling, ˜S. We have demonstrated that neglecting electromomentum coupling leads to estimates for the piezoelectric coefﬁcient that violate reciprocity and energy conservation (if the scatterer is asymmetric). The equations that the short circuit conﬁguration provide do not allow the extraction of the electro- momentum coefﬁcient and the dynamic effective dielectric permittivity; these will be determined using scattering experiments on the heterogeneous scatterer in the open circuit conﬁguration. The open circuit conﬁguration will also demonstrate that the retrieved properties using the equivalent description are in- consistent as they depend on the electric boundary conditions, even in the long-wavelength limit. This observation further reinforces the necessity of including electromomentum coupling in heterogeneous piezoelectric media with subwavelength-scale asymmetry. 4.2 Retrieved properties from an open circuit experiment Consider next the same scattering problems as in Sec. 4.1, only now the stack is in an open circuit conﬁgu- ration, i.e., the switch is open in Fig. 1. We repeat the same analysis described in Sec. 4.1 now accounting forD= 0. The objective is to retrieve the homogenized material properties in terms of the reﬂection and transmission coefﬁcients in open circuit conditions, now denoted by rLo,rRo,to. Accordingly, we modify Eq. (15) when considering a null D-ﬁeld, and obtain [cf. Eq. (17)]  u(x3) σ(x3) =To u(x0) σ(x0) ,To= TL 11TL 13 TL 31TL 33 . (32) We relate the transfer matrix Toto the reﬂection and transmission coefﬁcients according to Eqs. (19)-(21) and ﬁnd that for experiment L the coefﬁcients rLoandtLosatisfy Eqs. (21b) and (21a) when Toreplaces Ts. The analysis for experiment R mirrors the procedure described above, showing that rRoandtRoequal, respectively, to the right-hand side of Eqs. (21b) and (21a), where To/parenleftBig TL/parenrightBig is replaced by To/parenleftBig TR/parenrightBig . It is clear that the scattering coefﬁcients in this experiment are different from the coefﬁcients in the short circuit experiment, in view of the difference between ToandTs. Again, we ﬁnd that tLo≡tRo, as it should, while rLoandrRodiffer by a phase. As in Sec. 4.1, we require the properties of the ﬁctitious homogeneous medium to exhibit the same scattering coefﬁcients as the stack when subjected to the same electrical boundary conditions. We will develop the equations using the effective model, and recall that the equations for the equivalent model are obtained by setting ˜W= 0. Since we consider here an open circuit conﬁguration, the D-ﬁeld in the ﬁctitious homogeneous medium is null, and hence its constitutive equation for Dprovides the relation φ,x=˜A−1/parenleftBig˜Bu,x−iω˜Wu/parenrightBig . (33) 14 Subsequently, we employ this relation to replace φ,xin the constitutive relations for σandp, which now read σ=/parenleftBigg ˜C+˜B2 ˜A/parenrightBigg u,x−iω/parenleftBigg ˜S+˜B ˜A˜W/parenrightBigg u, (34a) p=/parenleftBigg ˜S†+˜B ˜A˜W†/parenrightBigg u,x−iω/parenleftBigg ˜ρ+˜W˜W† ˜A/parenrightBigg u. (34b) Finally, we repeat Eqs. (23)-(25) using the transfer matrix that is based on Eq. (34), and denote its logarithm by Mo. The resultant relationships for the effective dynamic properties are ˜S+˜B ˜A˜W=Mo 11 iωMo 12=Z0(rRo−rLo) (t2 o−rLorRo)eik0l−e−ik0l−rRo−rLo, (35a) ˜S†+˜B ˜A˜W†=−Mo 22 iωMo 12≡˜S+˜B ˜A˜W, (35b) ˜C+˜B2 ˜A=Mo−1 12=2iωZ 0ltoψ−1 osinψo (t2 o−rLorRo)eik0l−e−ik0l−rRo−rLo, (35c) ˜ρ+˜W˜W† ˜A=detMo ω2M12=−iρ0ψo 2tok0lsinψo /parenleftBig (t2 o−rLorRo)eik0l−e−ik0l/parenrightBig2−4rRorLo (t2 o−rLorRo)eik0l−e−ik0l−rRo−rLo , (35d) where cosψo=(t2 o−rLorRo)eik0l+e−ik0l 2to. (36) Equation (35) exposes the electromomentum effect on the scattering problem as it is represented by the extracted effective properties. Equation (35a) shows that the electromomentum coupling contributes to the asymmetry in the reﬂection, in addition to the asymmetry associated with the Willis coefﬁcient, ˜S. Further, Eq. (35d) shows that the presence of electromomentum coupling modiﬁes the apparent effective dynamic mass density. These phenomena are in agreement with the analysis in Sec. 3. The above analysis was based on the effective model that includes ˜W. The procedure to retrieve the homogenized properties using the equivalent model is the same up to setting ˜W= 0. Accordingly, the equivalent properties ˜Seq,˜S†eqand˜ρeqequal to the right-hand side of Eqs. (35a), (35b) and (35d), respectively. This results with second inconsistency in the model that neglects the electromomentum coupling: the equivalent properties depend on the electrical boundary conditions, since Eq. (28) and Eq. (35) provide different values for the foregoing properties. If we denote equivalent properties that are extracted from short- and open circuit conditions by subscript sando, respectively, then ˜Seq s/negationslash= ˜Seq o,˜S†eq s/negationslash=˜S†eq oand˜ρeq s/negationslash= ˜ρeq o. This type of dependency of the retrieved properties of inappropriate models on the experimental setup was noted before in Refs. [2, 45]. 15 Having pointed out the second inconsistency in the equivalent model which reinforces the need to include electromomentum coupling, we proceed to determine the effective electromechanical properties in the fully coupled medium. We ﬁrst observe that Eq. (35b) implies ˜W=˜W†, (37) where we used the fact that ˜S≡˜S†from Eq. (28b). While there are three additional equations in Eq. (35), the remaining three coefﬁcients ( ˜A,˜Band˜W) cannot be determined uniquely, since they enter the equations as products. Therefore, one more independent equation is required, which is obtained from the condition on the spatial average of the constitutive equation for D(cf. Sec. 4.1.1): 0 = ˜B/angbracketleftu,x/angbracketright−˜A/angbracketleftφ,x/angbracketright−iω˜W/angbracketleftu/angbracketright. (38) [Equivalently, as the third equation it is possible to use the corresponding equation from the short circuit experiment, namely, Eq. (29). In fact, we use Eq. (29) as a consistency check for the effective properties to be independent of the experimental setup.] Having at hand the solutions for the ﬁeld variables uandφ, we use Eq. (38) together with Eq. (35) to obtain expressions for the rest of the effective properties. We omit these lengthy expressions here, and analyze their low frequency expansions in Sec. 5. Before doing so, we emphasize that there are two assumptions in the procedure described above. ( i) We assume that the effective properties relate volume averages instead of ensemble averages; ( ii) We assume that the retrieved quantities—which depend on the surface data—indeed characterize the effective properties. These assumptions are valid only in the low frequency, long-wavelength ( ˜kl/lessmuch1) regime [30, 46]. In this regime, the effective properties that result from the process are independent of the experimental setup, i.e., independent of the electrical boundary conditions. By contrast, the equivalent properties depend on these conditions, even in the low frequency, long-wavelength limit. We have validated the consistency of the effective model by comparing the effective properties computed here with the effective properties that are calculated from another electrical setup. In this third setup, the piezoelectric element is subjected to a time harmonic voltage drop, while its boundaries are kept clamped (see analysis in Appendix B). In the next section, we provide the low frequency approximations of the effective properties using Taylor expansions. As per the discussion above, the leading terms in these expansions for the effective properties are independent of the experimental setup. 16 5 Low frequency approximations 5.1 Peano expansion To derive explicit expressions in the low frequency, long-wavelength limit, it proves useful to ﬁrst analyze the reﬂection and transmission coefﬁcients. To this end, we use the Peano expansion [15] for the transfer matrix of the layers Tα=I4×4+l/angbracketleftM/angbracketright+O/parenleftBig ω2/parenrightBig ,Mn= 0 0C−1 DnBn AnCDn 0 0Bn AnCDn−Cn AnCDn −ρnω20 0 0 0 0 0 0 , (39) where we recall that αdenotes the direction of the incident waves, and note that the terms of order ω2 depend on this direction. Using this expression in Eq. (21) provides the following power series expansions for the coefﬁcients in the short circuit experiment: rR≈rL=ik0l 2/bracketleftbigg ρ−1 0/angbracketleftρ/angbracketright−C0/parenleftbigg/angbracketleftBig C−1 D/angbracketrightBig +γ/angbracketleftbiggB ACD/angbracketrightbigg/parenrightbigg/bracketrightbigg +O/parenleftBig ω2/parenrightBig , γ =/angbracketleftbiggB ACD/angbracketrightbigg/angbracketleftbiggC ACD/angbracketrightbigg−1 (40a) rL−rR=ω22/summationdisplay p=13/summationdisplay q=p+1lplq CDpCDq/bracketleftBigg ∆qp/parenleftBig Z2 D/parenrightBig +γBpBq ApAq∆qp/parenleftbiggA BZ2 D/parenrightbigg/bracketrightBigg +O/parenleftBig ω3/parenrightBig , (40b) t= 1 +ik0l 2/bracketleftbigg ρ−1 0/angbracketleftρ/angbracketright+C0/parenleftbigg/angbracketleftBig C−1 D/angbracketrightBig +γ/angbracketleftbiggB ACD/angbracketrightbigg/parenrightbigg −2/bracketrightbigg +O/parenleftBig ω2/parenrightBig , (40c) where ∆qpf:=fp−fqis the difference between the value of fin layerpand layerq. Carrying out the same procedure for the coefﬁcients of the open circuit experiment shows that the coefﬁcients have the same form as in Eq. (40), only now the terms that involve products with γdo not appear. Equation (40b) shows that the leading term in the difference between the right- and left reﬂection coefﬁcients is of order ω2, and depends on the asymmetry in the electromechanical properties of the scatterer. 5.2 Explicit low frequency expressions Using the expression for the transmission and reﬂection coefﬁcients, we derive the Taylor expansion9of the effective properties [Eqs. (28), (B.10) and (B.13)] about ω= 0. After some rearrangement of terms 9We used Wolfram Mathematica Series function. 17 we obtain ˜ρ=/angbracketleftρ/angbracketright+O/parenleftBig ω2/parenrightBig , (41a) ˜C=/bracketleftbigg/angbracketleftBig C−1 D/angbracketrightBig +γ/angbracketleftbiggB ACD/angbracketrightbigg/bracketrightbigg−1 +O/parenleftBig ω2/parenrightBig ,whereγ=/angbracketleftbiggB ACD/angbracketrightbigg/angbracketleftbiggC ACD/angbracketrightbigg−1 , (41b) ˜A=/bracketleftBigg/angbracketleftBig A−1/angbracketrightBig −/angbracketleftBiggB2 A2CD/angbracketrightBigg +/angbracketleftBig C−1 D/angbracketrightBig−1/angbracketleftbiggB ACD/angbracketrightbigg2/bracketrightBigg−1 +O/parenleftBig ω2/parenrightBig , (41c) ˜B=/angbracketleftbiggB ACD/angbracketrightbigg /angbracketleftBig A−1/angbracketrightBig/angbracketleftBig C−1 D/angbracketrightBig −l−22/summationdisplay p=13/summationdisplay q=p+1lplq CDpCDq/parenleftbigg ∆qp/parenleftbiggB A/parenrightbigg/parenrightbigg2 −1 +O/parenleftBig ω2/parenrightBig , (41d) ˜S=iω 2l˜C2/summationdisplay p=13/summationdisplay q=p+1lplq CDpCDq/bracketleftBigg ∆qp/parenleftBig Z2 D/parenrightBig +γBpBq ApAq∆qp/parenleftbiggA BZ2 D/parenrightbigg/bracketrightBigg +O/parenleftBig ω2/parenrightBig , (41e) ˜W=iω 2l2˜B /angbracketleftbiggB ACD/angbracketrightbigg2/summationdisplay p=13/summationdisplay q=p+1lplq CDpCDq∆qp/parenleftbiggB A/parenrightbiggq/summationdisplay r=pρrlr(2−δpr−δqr) +O/parenleftBig ω2/parenrightBig , (41f) whereδpr= 1ifp=r, and zero otherwise. We observe that the leading term of all the effective properties except ˜Sand˜Wis independent of ωand the ordering of the layers. These terms agree with their expected values from static homogenization in the limiting decoupled case. Thus, the mass density reduces to the arithmetic mean, while when B= 0the stiffness and the permittivity reduce to the harmonic mean. In contrast, the leading term in the expansions of the Willis coupling and the electromomentum cou- pling is pure imaginary and linear in ω. Furthermore it depends on the ordering of the layers, such that ˜S and˜Wﬂip sign when layers 1 and 3 are interchanged. Since these couplings capture the reﬂection asym- metry, they inherit the dependency on the asymmetry in the electromechanical properties of the scatterer that Eq. (40b) exhibits. We also note that Eq. (40b) can be manipulated to a form that depends explicitly onZEandρBof each layer, with which ˜Scan be written as ˜S=iω 2lχ/braceleftBigg l1l2/bracketleftBigg α3∆21(Z2 E) +β3∆21(Bρ) +A2B1l3 ρ3(B1ρ1C3ρ3−B3ρ3C2ρ2) +A1B2l3 ρ3(B3ρ3C1ρ1−B2ρ2C3ρ3)/bracketrightBigg +l1l3/bracketleftBigg α2∆31(Z2 E) +β2∆31(Bρ) +A3B1l2 ρ2(B1ρ1C2ρ2−B2ρ2C3ρ3) +A1B3l2 ρ2(B2ρ2C1ρ1−B3ρ3C2ρ2)/bracketrightBigg +l2l3/bracketleftBigg α1∆32(Z2 E) +β1∆32(Bρ) +A3B2l1 ρ1(B2ρ2C1ρ1−B1ρ1C3ρ3) +A2B3l1 ρ1(B1ρ1C2ρ2−B3ρ3C1ρ1)/bracketrightBigg/bracerightBigg +O/parenleftBig ω2/parenrightBig , (42) with some property-dependent constants {χ,αj,βj,j= 1,2,3}that are given in Appendix C. Equation (42) explicitly shows that ˜Sdepends on the asymmetry in the elastic impedance Z2 E=ρCand a kind of “piezoelectric impedance”: ρB. This expression nicely reduces to an expression that depends only on 18 the elastic impedance when Bis zero. 5.2.1 Reduction to two layers When a scattering element with two layers is used to construct a periodic bulk material, it will always have inversion symmetry. Therefore, to generate ˜Sand˜WPernas-Salomón and Shmuel [38] broke spatial symmetry in their one-dimensional example using a periodic cell that is made of three layers. To enable a comparison with previous work, we have considered here a three-layer stack in the present work. (Note that our calculations show that in the low frequency regime, the effective properties ˜C,˜ρ,˜Band˜Afor three-layer elements and their periodic counterpart are the same.) However, it is clear that a two-layer scatterer in a background medium is sufﬁcient to break inversion symmetry, and that expressions for bilayer elements are obtained as a particular case of the above results by setting l3= 0. In that case, the expressions for ˜C,˜ρ,˜A, and ˜Bremain the same, only now the averaging is over two layers. The expressions for ˜Sand˜W, on the other hand, simplify to the following expressions ˜S=iω 2l˜Cl1l2[A1A2/angbracketleftA−1/angbracketright∆21(Z2 E) +B1B2/angbracketleftB−1/angbracketright∆21(Bρ)] A1A2CD1CD2/angbracketleftbiggC ACD/angbracketrightbigg +O/parenleftBig ω2/parenrightBig , (43a) ˜W=iω 2l˜Bl1l2/angbracketleftρ/angbracketright∆21/parenleftBig B A/parenrightBig CD1CD2/angbracketleftbiggB ACD/angbracketrightbigg+O/parenleftBig ω2/parenrightBig . (43b) Eq. (43) reproduces the same functional dependency of ˜Sand˜Was the three-layer expressions, in a sim- pler form. Eq. (43a) shows that the Willis coupling depends on the asymmetry in the elastic impedance, ρC, as it is well known from the elastic case, and on the asymmetry in the piezoelectric-like impedance, ρB. Eq. (43b) shows that the electromomentum coupling depends on the asymmetry in the ratio between the piezoelectric coefﬁcient and the permittivity, B/A . Notably, the sign of ˜Sand˜Wchanges when the ordering of the layers is reversed with respect to the coordinate system, since ∆12(·) =−∆12(·). The fact that ˜Sand˜Wdepend on asymmetries of different physical quantities has a very important implication: certain compositions of the piezoelectric element may exhibit zero Willis coupling and ﬁnite electromomentum coupling. This implies that we turn on and off the directional phase asymmetry by the opening and closing the circuit, respectively. One such design rule is ρ1C1=ρ2C2,ρ1B1=ρ2B2and B1A1−1/negationslash=B2A2−1. We note that we have implemented this design rule (i.e., ρCandρBare spatially constant) for a trilayer periodic medium in the rigorous periodic homogenization scheme of Ref. [38]. Interestingly, we found that it also works in the low frequency regime of the periodic case, namely, we received a zero Willis coupling and a nonzero the electromomentum coupling. 19 5.2.2 Generalization to an arbitrary number of layers and the continuous limit The low frequency expressions of Eq. (41) for ˜ρ,˜C,˜A,˜B,˜Sand˜Wcan be generalized to the case of an arbitrary number of n>3layers. The only change needed is to replace the upper limits, 2and3, of the sums overpandqin Eq. (41) by n−1andn, respectively. The continuous limit, i.e., when the properties of the element vary smoothly ( fp→f(x)), corresponds to n→∞ , in which case the sums are replaced by integrals and the average of a continuous property is simply /angbracketleftf/angbracketright=l−1/integraltext1 0fdx. In summary, in the continuous limit ˜ρ,˜Cand˜Aare exactly as shown in Eq. (41), while the expressions for ˜B,˜Sand˜Ware replaced by ˜B=/angbracketleftbiggB ACD/angbracketrightbigg /angbracketleftBig A−1/angbracketrightBig/angbracketleftBig C−1 D/angbracketrightBig −l−2/integraldisplayl 0dx CD(x)/integraldisplayl xdy CD(y)/parenleftBiggB(x) A(x)−B(y) A(y)/parenrightBigg2 −1 +O/parenleftBig ω2/parenrightBig ,(44a) ˜S=−iωl˜C/angbracketleftbigg/angbracketleftbigg ρD,1 2CD/parenleftbigg 1 +γB A/parenrightbigg/angbracketrightbigg/angbracketrightbigg +O/parenleftBig ω2/parenrightBig , (44b) ˜W=iω l2˜B /angbracketleftbiggB ACD/angbracketrightbigg/integraldisplayl 0dx CD(x)/integraldisplayl xdy CD(y)/parenleftBiggB(x) A(x)−B(y) A(y)/parenrightBigg/integraldisplayy xρ(z)dz+O/parenleftBig ω2/parenrightBig , (44c) where the double average /angbracketleft/angbracketleftf, g/angbracketright/angbracketrightof two material functions fandgis deﬁned as /angbracketleft/angbracketleftf, g/angbracketright/angbracketright=l−2/integraldisplayl 0f(x) dx/integraldisplayx 0g(y) dy−l−2/integraldisplayl 0g(x) dx/integraldisplayx 0f(y) dy. (45) Note that the double average is asymmetric in the two material parameters: /angbracketleft/angbracketleftf, g/angbracketright/angbracketright=−/angbracketleft/angbracketleftg, f/angbracketright/angbracketright. Also, the identity Z2 D=ρDCDhas been used to simplify ˜Sin Eq. (44b). 5.2.3 The source of the directional imaginary part of ˜Beq Using low frequency expansions, we can also explicitly identify the source of the directional imaginary part of ˜Beq. To this end, we expand Eq. (31) in power series about ω= 0, and obtain the low frequency approximation ˜Beq L=˜B−˜W˜Bγ−1Z−1 0,˜Beq R=˜B+˜W˜Bγ−1Z−1 0, (46) where ˜Band ˜Ware given in Eqs. (41d) and (41f), respectively. Thus, ˜Beqis a sum of a real quantity (˜B) and an imaginary quantity ( ˜W), where the latter multiplies a sign that depends on the incident wave direction. Thus, for the equivalent model to recover the scattering response of the element, it must absorb the directional phase change that ˜Wmodels into ˜Beq, thereby leading to a piezoelectric coefﬁcient that violates reciprocity and energy conservation. (For similar absorbance of a Willis coupling into an equivalent property see Eq. (20) in Ref. [37]. Indeed, the resultant equivalent properties there—although not referred to as such there—violate causality.) 20 Phase C(GPa)ρ(kg/m3)B(C/m2)A(nF/m ) PZT4 115 7500 15.1 5.6 BaTiO 3 165 6020 3.64 0.97 PVDF 12 1780 -0.027 0.067 PMMA 3.3 1188 0 0.023 Fictitious phase 143.27 6020 18.81 0.97 Table 1: Physical properties of the materials that are used in the computations as the constituents of the element. 6 Numerical experiments In this Section, we apply the procedure in Sec. 4 to numerically calculate the effective properties of differ- ent piezoelectric elements connecting two elastic waveguides made of PMMA. The physical properties of the PMMA, as well as the piezoelectric layers to be considered are given in Table 1. The following calculations are presented as functions of the ordinary frequency f=ω/2πover the range 0−20kHz, in which we have found that the low frequency expansions and the original expressions practically coincide. (The difference is less than 0.1%.) The ﬁrst element that we analyze is made of 2.62 mm PZT4 (layer 1), 0.19 mm BaTiO 3(layer 2) and 0.19 mm PVDF (layer 3). We also analyze the properties that are retrieved when the element is inverted (or rotated 180◦), i.e., when layers 1 and 3 are interchanged. Figure 2a shows the pure imaginary Willis coupling. Values retrieved from the element ( ˜S⊿, black squares) and its inversion ( ˜S/triangleleft, blue squares) satisfy ˜S⊿=−˜S/triangleleft. Figure 2b shows the pure imaginary electromomentum coupling. Like ˜S,˜Wdepends on the orientation (or polarization) of the element: values retrieved from the element ( ˜W⊿, black circles) and its inversion ( ˜W/triangleleft, blue circles) satisfy ˜W⊿=−˜W/triangleleft. Figure 2c presents the values of the standard effective properties, where each property is normalized by its arithmetic mean (normalized properties are denoted by overhead hat). Speciﬁcally, the brown, red, green and black curves correspond to ˆρ,ˆC,ˆAandˆB, respectively. The computation delivers real constants that are independent of the direction of the incident wave, orientation of the element and the circuit conditions, as they should be. We compute next the equivalent properties and focus on the retrieved piezoelectric coefﬁcient ˜Beqas per the short circuit experiment [Eq. (31)]. In contrast with the purely real, direction-independent value of˜B,˜Beqis complex-valued and direction dependent. (We recall that values extracted from the ﬁelds that are generated by leftward- and rightward incoming waves are denoted by ˜Beq Land˜Beq R, respectively.) Fig. 2d depicts the real and imaginary parts of ˜Beq L(˜Beq R) by the grey (blue) curve and green (red) diamonds, respectively. The direction-dependent imaginary part of ˜Bimplies that the equivalent model violates energy conservation and reciprocity, hence losing its physical meaning, even though the model recovers the scattering response of the element in the short circuit experiment. This occurs because ˜Beqabsorbs precisely the directional phase change that the electromomentum couplings models. 21 (a) □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □S∼⊳ S∼⊲x 0 5 10 15 200 -0.50.5 f[kHz]S~ [Mkg/m2s] (b) ○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○ ○W∼⊲W∼⊳x0 5 10 15 200 -1010 f[kHz]W~ [μCs/m3] (c) CABρ 0 5 10 15 2000.511.5 f[kHz]Normalized values (d) ◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇ ◇B~ Req B~ Leqx0 5 10 15 200 -0.50.51.12 f[kHz]B~eq[C/m2] FIG. 2: Retrieved properties of the PZT4-BaTiO-PVDF element and its inversion (see insets) as functions of the ordinary frequency. (a) The pure imaginary Willis coupling. Values retrieved from the element ( ˜S⊿, black squares) and its inversion ( ˜S/triangleleft, blue squares) satisfy ˜S⊿=−˜S/triangleleft. (b) The pure imaginary electromomentum coupling. Values retrieved from the element ( ˜W⊿, black circles) and its inversion ( ˜W/triangleleft, blue circles) satisfy ˜W⊿=−˜W/triangleleft. (c) The standard effective properties, each one is normalized by its arithmetic mean. Brown, red, green and black curves correspond to ˆρ,ˆC,ˆAandˆB, respectively. (d) The equivalent piezoelectric coefﬁcient ˜Beqthat is retrieved from the short circuit experiment. Values extracted from the ﬁelds that are generated by leftward- and rightward incident waves are denoted by ˜Beq sLand˜Beq sR, respectively. While Re ˜Beq sL(blue curve) is equal to Re ˜Beq sR(gray curve), we observe that Im ˜Beq sL(green diamond marks) is nonzero and different from Im ˜Beq sR(red diamond marks), thereby violating energy conservation and reciprocity. 22 (a) Short circuitOpen circuit0 5 10 15 2000.0050.01 f[MHz]Arg(r R)-Arg(r L) (b) ○○○○○○○○○○○○○○○○○○○○□□ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □□○○○○○○○○○○○○○○○○○○○○○□□ □ □ □ □ □ □ □ □ □ □ □ □□○ □ ○ □S∼⊳W∼⊳ S∼⊲x xW∼⊲ 0 5 10 15 200 -40.540.5 f[kHz]S~ ,W~ [μCs/m3] (c) △△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △△△△ △ △ △ △ △ △ △ △ △ △ △ △△○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○△ △○ ○S∼ oeq⊳ S∼ oeq⊲ S∼ seq⊳ S∼ seq⊲x x0 5 10 15 200 -0.30.3 f[kHz]S~eq[Mkg/m2s] (d) ◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇ B~ Req B~ Leqx 0 5 10 15 20-0.500.5 f[kHz]Beq~ [C/m2]0 10 2016.75ReBeq~ FIG. 3: Retrieved properties of the bilayer. (a) The difference in the phase of the right and left reﬂection coefﬁcients in the open circuit (solid blue line) and short circuit (dashed black line). (b) The pure imaginary ˜Sand˜W. Same legend as in Fig. 2. (b) The pure imaginary ˜Seqas retrieved for the element (green markers) and its inversion (red markers) when the circuit is open ( ˜Seq o, circle marks) and shorted ( ˜Seq s, triangle marks). The fact that ˜Seq s/negationslash=˜Seq o implies that the equivalent model is inconsistent. (c) The equivalent piezoelectric coefﬁcient ˜Beqthat is retrieved from the short circuit experiment. Same legend as in Fig. 2. Again, the direction-dependent imaginary part of ˜Beq is a violation of reciprocity and energy conservation. As mentioned previously, there is another inconsistency in the equivalent model: the Willis coupling that the equivalent model retrieves from the short circuit experiment is different from the value that equivalent model retrieves for the open circuit experiment. For this piezoelectric element, the difference between the two results is minor and not presented here. To design an element in which this difference is prominent, we make use of the conclusion in Sec. 5, namely, that ˜Sis zero when ρCandρBare constant in the element, and that ˜Wis nonzero when B/A is asymmetric in the element. Thus, the second element that we consider is a bilayer made of 1.5mm PZT4 and a 1.5mm of ﬁctitious material, the properties of which are set according to the above rule (see Tab. 1 for the corresponding values). We present ﬁrst in Fig. 3a the difference in the phase of the right and left reﬂection coefﬁcients in the open circuit (solid blue line) and short circuit (dashed black line). Indeed, there is no phase difference in the short circuit experiment, and hence both models will retrieve a null Willis coefﬁcient. However, since 23 there is a phase difference in the open circuit experiment, the equivalent model will need a nonzero ˜Seqto capture this, in contrast with the prediction of the equivalent model in the short circuit experiment. These observations are demonstrated in the remaining panels. Speciﬁcally, Fig. 3b shows the effective ˜Sand ˜Wof the element and the inverted element, with the same legend as in Fig. 2. As expected, the former is zero and the latter is nonzero. Fig. 3c shows the values of ˜Seqas retrieved when the circuit is open ( ˜Seq o, circle marks) and shorted ( ˜Seq s, triangle marks), where the green and red colors denote values retrieved for the element (superscript ⊿) and the inverted element (superscript /triangleleft). As expected, we observe that ˜Seq o is different from ˜Seq s: the latter is identically zero and the former is not, and hence the equivalent model is inconsistent. Lastly, we show in 3d the calculation of ˜Beqas extracted from the short circuit experiment [Eq. (31)], and use the same legend as in Fig. 2. Here again, the equivalent piezoelectric coefﬁcient has an imaginary part that depends on the direction of the incident wave. 7 Summary By adapting the homogenization theory of Willis [51, 52, 54, 55, 56] to piezoelectric composites, Pernas- Salomón and Shmuel [38] have found that a macroscopic cross-coupling between the momentum and electric ﬁeld emerges if the media has broken inversion symmetry. The electromomentum coupling con- stitutes an additional degree of freedom in the design space of metamaterials, which reﬂects an unusual interplay between the kinetic energy and electric energy at the macroscale. The objective of this work was to provide a heuristic demonstration of this coupling, analyze its effect, and elucidate its physical origins. Towards this end, we have analyzed ﬁrst the effect of the electromomentum coupling on plane waves in homogenized media under short circuit and open circuit conditions. This analysis shows that only un- der open circuit conditions, the coupling modiﬁes the phase velocity and introduces a directional phase angle to the characteristic elastic impedance. The former feature is similar to the effect that the piezo- electric coupling has on the phase velocity, and the latter feature is similar to the effect that the Willis coupling has on the characteristic impedance. The analysis suggests that these features can be quantiﬁed from the scattering properties of a one- dimensional asymmetric piezoelectric scatterer. To do so, we have calculated analytically its scattering coefﬁcients in response to rightward- and leftward incident waves, under the two circuit conditions. We then used an inverse program to determine the properties of a homogenized medium according to two models: an effective model that includes both Willis and electromomentum coupling, and an equivalent model that neglects electromomentum coupling [13, 30, 45]. This heuristic homogenization approach is applicable only in the long-wavelength, low frequency regime [46]. The conclusions that we have found and demonstrated numerically correspond to this limit, are summarized next. We highlight the main observation from this analysis: the reﬂection of rightward- and leftward inci- 24 dent waves is different by a phase, where this phase depends on the electric circuit conditions, i.e., if the circuit is shorted or open. It is precisely the extra phase in the open circuit that only the electromomen- tum coupling captures. Together with this coupling, the properties that are retrieved from the scattering experiments using the effective model satisfy reciprocity and energy conservation, and are independent of the excitation setup. By contrast, the properties of the equivalent model that are retrieved from the short circuit setup violate reciprocity and energy conservation, and are different from the values that are retrieved from the open circuit setup, even in the low frequency regime. We note that such violations of physical laws by equivalent models that neglect the electromomentum coupling were also found in studies of inﬁnite piezoelectric media, see Ref. [28]. We have obtained explicit expressions for the low frequency approximations of the effective prop- erties using Taylor expansions. The leading term in the expansions of the effective mass, stiffness and permittivity is a real constant that agrees with their expected value from static homogenization in the limiting decoupled case. By contrast, the leading term in the expansions of the Willis coupling and the electromomentum coupling is linear and purely imaginary, hence responsible for the phase change. The expansion of the former clearly shows that the Willis coupling depends on the asymmetry in the impedance—as known in the elastic case— andthe asymmetry in a kind of a piezoelectric impedance, that is the product of the piezoelectric coefﬁcient by the mass density. The leading term in the expansion of the electromomentum coupling depends on the asymmetry in ratio betwen the piezoelectric constant and the permittivity. Hence, it is possible to design elements that have zero Willis coupling and ﬁnite electromomentum coupling. The directional phase angle of such element vanishes in the short circuit setup—since it depends only on the Willis coupling—and persists in the open circuit setup since it de- pends on the electromomentum coupling. This observation suggests a convenient way to control the directional phase change by the opening and shorting of the circuit. We expect the demonstration of the electromomentum effect in this simple setting to promote fur- ther studies of its implications in more general settings, and its engineering for applications that require control over elastic waves. Acknowledgments This research was supported by the Israel Science Foundation, funded by the Israel Academy of Sci- ences and Humanities (Grant no. 2061/20), the United States-Israel Binational Science Foundation (Grant no. 2014358), and Ministry of Science and Technology (grant no. 880011). A.N.N. and M.R.H. acknowl- edge support from NSF EFRI award no. 1641078 and M.R.H. acknowledges additional support ONR YIP award no. N00014-18-1-2335. 25 Appendix A. The transfer matrix for the mechanical ﬁelds of the homogenized medium in short circuit. Consider a ﬁnite slab of homogeneous medium governed by Eq. (3), admitting plane waves of the form u(x,t) =Ue±i˜kEx−iωt. In a short circuit conﬁguration such that the electric ﬁeld vanishes, its stress and momentum reduce to the Willis form σ=˜Cu,x−iω˜Su, (A.1) p=˜S†u,x−iω˜ρu. (A.2) We extract the strain from Eq. (A.1): u,x=˜C−1σ+iω˜S˜C−1u, (A.3) substitute it into Eq. (A.2), and then use the resultant expression for pin the governing equation σ,x= −iωpto obtain σ,x=−iω˜S†˜C−1σ+ω2/parenleftBig˜S˜S†˜C−1−˜ρ/parenrightBig u. (A.4) Equations A.3) and (A.4) are written in matrix form as  u,x σ,x =M u σ where M= iω˜S˜C−1 ˜C−1 ω2/parenleftBig˜S˜S†˜C−1−˜ρ/parenrightBig −iω˜S†˜C−1 . (A.5) Considering plane wave solutions, Eq. (A.5) leads to the following relation between the displacement and stress at two different points in the slab:  u(xb) σ(xb) =Tm(xb,xa) u(xa) σ(xa) , x 0≤xa,xb≤x3, (A.6) where Tm(xb,xa) =e(xb−xa)M(A.7) is the transfer matrix relating the displacement and stress ﬁelds at the space coordinates xaandxb. 26 Appendix B. Retrieved properties in the clamped ends conﬁguration. This appendix provides an analysis of a third setup that veriﬁes the consistency of the effective properties. We consider a setup with zero average strain by clamping the ends x0andx3, and excite the stack by applying a time harmonic ( e−iωt) voltage drop. (This conﬁguration is consistent with an effective stress- charge form since we evaluate the effective permittivity and piezoelectric coefﬁcients at a prescribed effective strain when an electric ﬁeld generates stress and electric displacement ﬁeld.) We determine the effective properties from the requirement that the ﬁctitious medium will deliver the same average ﬁelds as in the piezoelectric stack. To this end, we ﬁrst calculate the exact ﬁelds in the stack, as follows. From Eq. (15) we obtain σ(x0) =−T⊿ 14 T⊿ 13D, (B.1) which together with Eqs. (14) and (B.1) determine u(x)andσ(x)as u(x) =  D/parenleftBiggB1 A1−T⊿ 14 T⊿ 13/parenrightBiggsinkD1(x−x0) ωZD1, x 0≤x≤x1, u(x1) coskD2(x−x1) +/bracketleftbigg σ(x1) +B2 A2D/bracketrightbiggsinkD2(x−x1) ωZD2, x 1≤x≤x2, u(x2) coskD3(x−x2) +/bracketleftbigg σ(x2) +B3 A3D/bracketrightbiggsinkD3(x−x2) ωZD3, x 2≤x≤x3,(B.2) where T⊿=Tem 3Tem 2Tem 1and u(x1) =/parenleftBigg Tem 1(14)−Tem 1(13)T⊿ 14 T⊿ 13/parenrightBigg D, u (x2) =Tem 2(11)u(x1) +Tem 2(13)σ(x1) +Tem 2(14)D, (B.3) σ(x1) =/parenleftBigg Tem 1(34)−Tem 1(33)T⊿ 14 T⊿ 13/parenrightBigg D, σ (x2) =Tem 2(31)u(x1) +Tem 2(33)σ(x1) +Tem 2(34)D, (B.4) andTem n(ij)is theij-component of the matrix Tem n. The stress as function of xisσ(x) =CDnu,x− (Bn/An)D, forxn−1≤x≤xn,n= 1,2,3. Furthermore, from Eqs. (15) and (B.1) we have that φ(x3)−φ(x0) =/parenleftBigg T⊿ 24−T⊿ 23T⊿ 14 T⊿ 13/parenrightBigg D. (B.5) Having determined u(x)andσ(x)allows us to average over the stack, which together with Eq. (B.5) and 27 the fact that (in absence the of charge) Dis constant in the stack provide l/angbracketleftu/angbracketright= Γ⊿ 1/angbracketleftD/angbracketright, l/angbracketleftσ/angbracketright= Γ⊿ 2/angbracketleftD/angbracketright, l/angbracketleftφ,x/angbracketright= Γ⊿ 3/angbracketleftD/angbracketright, (B.6) where the terms Γ⊿ iare Γ⊿ 1=/parenleftBiggB1 A1−T⊿ 14 T⊿ 13/parenrightBigg1−coskD1l1 ωkD1ZD1 +/parenleftBigg Tem 1(14)−Tem 1(13)T⊿ 14 T⊿ 13/parenrightBiggsinkD2l2 kD2+/parenleftBigg Tem 1(34)−Tem 1(33)T⊿ 14 T⊿ 13+B2 A2/parenrightBigg1−coskD2l2 ωkD2ZD2 +/bracketleftBigg Tem 2(11)/parenleftBigg Tem 1(14)−Tem 1(13)T⊿ 14 T⊿ 13/parenrightBigg +Tem 2(13)/parenleftBigg Tem 1(34)−Tem 1(33)T⊿ 14 T⊿ 13/parenrightBigg +Tem 2(14)/bracketrightBiggsinkD3l3 kD3 +/bracketleftBigg Tem 2(31)/parenleftBigg Tem 1(14)−Tem 1(13)T⊿ 14 T⊿ 13/parenrightBigg +Tem 2(33)/parenleftBigg Tem 1(34)−Tem 1(33)T⊿ 14 T⊿ 13/parenrightBigg +Tem 2(34)+B3 A3/bracketrightBigg1−coskD3l3 ωkD3ZD3(B.7) Γ⊿ 2=/parenleftBiggB1 A1−T⊿ 14 T⊿ 13/parenrightBiggsinkD1l1 kD1−B1 A1l1 +/parenleftBigg Tem 1(34)−Tem 1(33)T⊿ 14 T⊿ 13+B2 A2/parenrightBiggsinkD2l2 kD2−ωZD2/parenleftBigg Tem 1(14)−Tem 1(33)T⊿ 14 T⊿ 13/parenrightBigg1−coskD2l2 kD2−B2 A2l2 +/bracketleftBigg Tem 2(31)/parenleftBigg Tem 1(14)−Tem 1(13)T⊿ 14 T⊿ 13/parenrightBigg +Tem 2(33)/parenleftBigg Tem 1(34)−Tem 1(33)T⊿ 14 T⊿ 13/parenrightBigg +Tem 2(34)+B3 A3/bracketrightBiggsinkD3l3 kD3 −ωZD3/bracketleftBigg Tem 2(11)/parenleftBigg Tem 1(14)−Tem 1(13)T⊿ 14 T⊿ 13/parenrightBigg +Tem 2(13)/parenleftBigg Tem 1(34)−Tem 1(33)T⊿ 14 T⊿ 13/parenrightBigg +Tem 2(14)/bracketrightBigg1−coskD3l3 kD3−B3 A3l3, (B.8) Γ⊿ 3=T⊿ 24−T⊿ 23T⊿ 14 T⊿ 13. (B.9) We are now in a position to determine the effective properties from the relations between average ﬁelds. Since we use the equations for the average stress and electric displacement ﬁeld when the average strain is zero, then the retrieved piezoelectric coefﬁcient is the same whether we use ˜Was in model (3), or neglect it as in the equivalent model; by contrast, the equivalent permittivity is different from the effective permittivity if ˜Wis neglected. When ˜Wis included, the effective properties read ˜A⊿=−l+iω˜W⊿Γ⊿ 1 Γ⊿ 3,˜B⊿=Γ⊿ 2+iω˜S⊿Γ⊿ 1 Γ⊿ 3; (B.10) here, we use superscript⊿to distinguish the results for this element from the results for an inverted element, the layers 1 and 3 of which are intechanged. Eq. (B.10) provides a third equation to determine 28 ˜W. Thus, together with Eq. (35a) we obtain ˜W⊿=−l/parenleftBig˜S⊿eq o−˜S⊿/parenrightBig Γ⊿ 2+iω˜S⊿eq oΓ⊿ 1(B.11) where ˜S⊿eq oequals the right-hand side of Eq. (35a). We note that if ˜Wis neglected, then the expressions for the equivalent permittivity from this experiment is determined by setting ˜W= 0in Eq. (B.10) 1. It is straightforward to carry out the same analysis if the element is inverted, i.e., if layers 1and3in the element are interchanged. In this case ˜A/triangleleftand˜B/triangleleftof the inverted element are given by ˜A/triangleleft=−l+iω˜W/triangleleftΓ/triangleleft 1 Γ/triangleleft 3,˜B/triangleleft=Γ/triangleleft 2+iω˜S/triangleleftΓ/triangleleft 1 Γ/triangleleft 3; (B.12) here, Γ/triangleleft 1,Γ/triangleleft 2andΓ/triangleleft 3correspond to the case where layers 1and3of the laminate are swapped, and can be obtained, respectively, from the expressions for Γ⊿ 1,Γ⊿ 2andΓ⊿ 3by interchanging the properties and matrices corresponding to the layers 1and 3, and replacing the matrix T⊿by the matrix T/triangleleft=Tem 1Tem 2Tem 3. Again, by substituting Eq. (B.12) into Eq. (35a) we obtain ˜W/triangleleft=−l/parenleftBig˜S/trianglelefteq o−˜S/triangleleft/parenrightBig Γ/triangleleft 2+iω˜S/trianglelefteq oΓ/triangleleft 1, (B.13) where ˜S⊿eq ois given by the right-hand side of Eq. (35a) for the inverted element. Now, since Eqs. (14) and (15) verify that Γ/triangleleft 1=−Γ⊿ 1,Γ/triangleleft 2= Γ⊿ 2andΓ/triangleleft 3= Γ⊿ 3, and Eqs. (28a) and (35a) imply that ˜S/triangleleft=−˜S⊿and ˜S⊿eq=−˜S/trianglelefteq, and hence the electromomentum coupling changes sign when the element is inverted, i.e., ˜W/triangleleft=−˜W⊿, like the Willis coupling. Together, the former imply that ˜A⊿=˜A/triangleleft,˜B⊿=˜B/triangleleft. This further The above analysis demonstrates another inconsistency in the equivalent model that neglects ˜W: the extracted values for the equivalent piezoelectric coefﬁcients are different in the short circuit experiment and in this experiment, i.e., Eq. (31) and Eq. (B.10) 2, are different. This is evident from the fact that Eq. (31) is complex-valued, while Eqs. (B.10) 2is real-valued. Appendix C. The coefﬁcients in the alternative expression of ˜S. The constants that appear in Eq. (42) are χ=/parenleftBig A1A2CD1CD2C3l3+A1A3CD1CD3C2l2+A2A3CD2CD3C1l1/parenrightBig−1˜C, 29 and α1=A2A3C1l1+A1CD1(A2l3+A3l2), β 1=B1B2B3l1+A1CD1(B3l2+B2l3) α2=A1A3C2l2+A2CD2(A1l3+A3l1), β 2=B1B2B3l2+A2CD2(B3l1+B1l3) α3=A1A2C3l3+A3CD3(A1l2+A2l1), β 3=B1B2B3l3+A3CD3(B2l1+B1l2)