Symmetry Breaking Creates Electro-Momentum Coupling in Piezoelectric Metamaterials

TL;DR

This paper introduces a homogenization method revealing new electro-momentum couplings in piezoelectric metamaterials caused by symmetry breaking, enabling advanced active wave control through external stimuli.

Contribution

It develops a rigorous homogenization framework to identify symmetry-breaking induced couplings in responsive metamaterials, specifically unveiling electro-momentum interactions in piezoelectric systems.

Findings

Discovered Willis-like couplings in piezoelectric metamaterials.

Developed a homogenization method for responsive metamaterials.

Demonstrated potential for active wave control via external stimuli.

Abstract

Metamaterials posses microstructure designed to acquire properties not found in nature. An epitome in acoustics and solid mechanics is Willis coupling, which refers to the particle velocity-stress coupling, and of great significance since it controls mechanical waves. We here reveal new couplings, analogous to Willis coupling, when considering patterns of materials that mechanically interact with magnetic, electric or thermal fields. To this end, we develop a rigorous homogenization method for the effective properties of such responsive metamaterials. As an example, we apply the scheme to piezoelectric materials, and unveil coupling of the velocity and electric fields. Hence, Willis-like couplings in responsive metamaterials open new avenues for active wave control by modulation of external stimuli.

Figures

Tables

Table 1. Table 1: Physical properties of the phases comprising the periodic piezoelectric laminate.

Phase	$C\left(\mathrm{GPa}\right)$	$\rho \left(\mathrm{kg}/m^3\right)$	$B\left(C/m^2\right)$	$A$$\left(\text{nF/m}\right)$
PZT4	$115$	7500	15.1	5.6
${\text{BaTiO}}_{\text{3}}$	$165$	6020	3.64	0.97
PVDF	12	1780	-0.027	$0.067$
${\text{Al}}_2{\text{O}}_3$	300	3720	0	$0.079$
PMMA	3.3	1188	0	$0.023$

Equations

$$\begin{array}{r}\nabla \cdot \mathbf{\sigma }+\mathbf{f}-\dot{\mathbf{p}}=0,\end{array}$$

$$\nabla \cdot \mathbf{D}=q,\nabla \times \mathbf{E}=0,$$

$$\left(\begin{array}[]{c}\bm{\sigma}\\ \mathbf{D}\\ \mathbf{p}\end{array}\right)=\left(\begin{array}[]{ccc}\mathbf{C}&\mathbf{B}^{\mathsf{T}}&0\\ \mathbf{B}&-\mathbf{A}&0\\ 0&0&\rho\end{array}\right)\left(\begin{array}[]{c}\nabla\mathbf{u}-\bm{\eta}\\ \nabla\phi\\ \dot{\mathbf{u}}\end{array}\right),$$

$$A_{ij}=A_{ji},B_{ijk}=B_{jik},C_{ijkl}=C_{jikl}=C_{jilk}=C_{klij},{\sigma }_{ij}={\sigma }_{ji}.$$

$$\begin{array}{r}\nabla \cdot ⟨\mathbf{\sigma }⟩+\mathbf{f}-\dot{⟨\mathbf{p}⟩}=0,\end{array}\nabla \cdot ⟨\mathbf{D}⟩=q,$$

$$\displaystyle\left(\begin{array}[]{c}\left\langle\bm{\sigma}\right\rangle\\ \left\langle\mathbf{D}\right\rangle\\ \left\langle\mathbf{p}\right\rangle\end{array}\right)=\left(\begin{array}[]{ccc}\tilde{\mathbf{C}}&\tilde{\mathbf{B}}^{\mathsf{T}}&\tilde{\mathbf{S}}\\ \tilde{\mathbf{B}}&-\tilde{\mathbf{A}}&\tilde{\mathbf{W}}\\ \tilde{\mathbf{S}}^{\dagger}&\tilde{\mathbf{W}}^{\dagger}&\tilde{\bm{\rho}}\end{array}\right)\left(\begin{array}[]{c}\left\langle\nabla\mathbf{u}\right\rangle-\bm{\eta}\\ \left\langle\nabla\phi\right\rangle\\ \left\langle\dot{\mathbf{u}}\right\rangle\end{array}\right),$$

$${\mathsf{D}}^{\mathsf{T}}\mathsf{h}=-\mathsf{f},$$

$$\mathsf{h}=\mathsf{L}\left(\mathsf{b}-\mathsf{m}\right),$$

$$\displaystyle\begin{aligned} &\mathsf{B}=\left(\begin{array}[]{cc}\nabla&0\\ 0&\nabla\\ s&0\end{array}\right),\mathsf{m}=\left(\begin{array}[]{c}\bm{\eta}\\ 0\\ 0\end{array}\right),\mathsf{D}=\left(\begin{array}[]{cc}\nabla\cdot&0\\ 0&\nabla\cdot\\ -s&0\end{array}\right),\\ &\mathsf{w}=\left(\begin{array}[]{c}\mathbf{u}\\ \phi\end{array}\right),\mathsf{f}=\left(\begin{array}[]{c}\mathbf{f}\\ -q\end{array}\right),\mathsf{b}=\mathsf{B}\mathsf{w},\end{aligned}$$

$$\mathsf{D}^{\mathsf{T}}\mathsf{L}\mathsf{B}\mathsf{G}=-\left(\begin{array}[]{cc}\delta\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\mathbf{I}&0\\ 0&\delta\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\end{array}\right),$$

$$\mathsf{G}=0\mathrm{over}\partial {\Omega }_{\mathsf{w}},{\left(\mathsf{BG}\right)}^{\mathsf{T}}\mathsf{LN}=0\mathrm{over}\partial {\Omega }_{\mathsf{t}},$$

$$\mathsf{N}^{\mathsf{T}}=\left(\begin{array}[]{ccc}\otimes\mathbf{n}&0&0\\ 0&\otimes\mathbf{n}&0\end{array}\right).$$

$$\begin{array}{rl}\mathsf{w}(\mathbf{x},y)& ={\int }_{\Omega }{\mathsf{G}}^{\mathsf{T}}\mathsf{f}\mathrm{d}\Omega +{\int }_{\partial {\Omega }_{\mathsf{t}}}{\mathsf{G}}^{\mathsf{T}}{\mathsf{N}}^{\mathsf{T}}\mathsf{L}\left(\mathsf{Bw}-\mathsf{m}\right)\mathrm{d}a\\ & +{\int }_{\Omega }{\left(\mathsf{BG}\right)}^{\mathsf{T}}\mathsf{Lm}\mathrm{d}\Omega -{\int }_{\partial {\Omega }_{\mathsf{w}}}{\left(\mathsf{BG}\right)}^{\mathsf{T}}\mathsf{LNw}\mathrm{d}a,\end{array}$$

$$\begin{array}{rl}\mathsf{w}(\mathbf{x},y)& ={\int }_{\Omega }{\mathsf{G}}^{\mathsf{T}}\mathsf{f}\mathrm{d}\Omega +{\int }_{\partial {\Omega }_{\mathsf{t}}}{\mathsf{G}}^{\mathsf{T}}{\mathsf{N}}^{\mathsf{T}}\mathsf{L}\left(\mathsf{Bw}-\mathsf{m}\right)\mathrm{d}a\\ & -{\int }_{\Omega }{\left(\mathsf{BG}\right)}^{\mathsf{T}}\mathsf{L}\left(⟨\mathsf{Bw}⟩-\mathsf{m}\right)\mathrm{d}\Omega +⟨\mathsf{w}⟩.\end{array}$$

$$\mathsf{w}=\left\{{\mathsf{G}}^{\mathsf{T}}{⟨\mathsf{G}⟩}^{-\mathsf{T}}⟨{\left(\mathsf{BG}\right)}^{\mathsf{T}}\mathsf{L}⟩-{\left(\mathsf{BG}\right)}^{\mathsf{T}}\mathsf{L}\right\}\left(⟨\mathsf{Bw}⟩-\mathsf{m}\right)+⟨\mathsf{w}⟩.$$

$$⟨\mathsf{h}⟩\left(\mathbf{x}\right)=\tilde{\mathsf{L}}\left(⟨\mathsf{Bw}⟩-\mathsf{m}\right),$$

$$\tilde{\mathsf{L}}=⟨\mathsf{L}⟩-⟨\mathsf{LB}{\left(\mathsf{BG}\right)}^{\mathsf{T}}\mathsf{L}⟩+⟨\mathsf{LB}{\mathsf{G}}^{\mathsf{T}}⟩{⟨\mathsf{G}⟩}^{-\mathsf{T}}⟨{\left(\mathsf{BG}\right)}^{\mathsf{T}}\mathsf{L}⟩.$$

$$⟨\zeta ⟩=\frac{1}{l}{\int }_{{\Omega }_{\mathrm{p}}}{\zeta }_0\left(x-y\right)\mathrm{d}y$$

$$D_{0,x}^{\left(n\right)}=0,$$

$$D_0^{\left(n\right)}=B_0^{\left(n\right)}{u}_{0,x}^{\left(n\right)}-A_0^{\left(n\right)}\text{ϕ0,x(n)}$$

$$\text{ϕ0,xx(n)}=\frac{B_0^{\left(n\right)}}{A_0^{\left(n\right)}}{u}_{0,xx}^{\left(n\right)}.$$

$${\sigma }_{0,x}^{\left(n\right)}-s^2{\rho }_0^{\left(n\right)}{u}_0^{\left(n\right)}=-f_0^{\left(n\right)},$$

$${\sigma }_0^{\left(n\right)}=C_0^{\left(n\right)}\left(u_{0,x}^{\left(n\right)}-{\eta }_0^{\left(n\right)}\right)+B_0^{\left(n\right)}{\varphi }_{0,x}^{\left(n\right)}$$

$$\left(C_0+\frac{B_0^2}{A_0}\right)\left(u_{0,xx}-{\eta }_{,x}\right)-s^2{\rho }_0{u}_0=-f_0,$$

$$\displaystyle\begin{aligned} G_{0}\left(x,x^{\prime}\right)&=\left\{\begin{array}[]{l}Vu^{+}\left(x\right)u^{-}(x^{\prime}),\quad x<x^{\prime},\\[10.0pt] Vu^{+}(x^{\prime})u^{-}(x),\ \quad x^{\prime}<x,\end{array}\right.\\ \left\langle G\right\rangle\left(x,x^{\prime}\right)&=\frac{1}{l}\int_{\Omega_{\mathrm{p}}}G_{0}\left(x-y,x^{\prime}-y\right)\mathrm{d}y,\end{aligned}$$

$$V^{-1}=\left(C_0+\frac{B_0^2}{A_0}\right)\left(u_{,x}^{+}{u}^{-}-u^{+}{u}_{,x}^{-}\right),u^{\pm }=u_{\mathrm{p}}^{\pm }\left(x\right)e^{\pm i{k}_{B}x},$$

$$\begin{array}{rl}u\left(x,y\right)=& \left(G⟨G⟩{}^{-1}⟨G_{,x^{\prime }}\check{C}⟩-G{}_{,x^{\prime }}\check{C}\right)\left(⟨u_{,x^{\prime }}⟩-\eta \right)\\ & +\left(G⟨G⟩{}^{-1}⟨G\rho ⟩-G\rho \right)s^2⟨u⟩+⟨u⟩,\end{array}$$

$$\begin{array}{rl}\tilde{S}\left(\xi \right)& ={\tilde{S}}^{†}\left(-\xi \right)=-\mathrm{conj}{\tilde{S}}^{†}\left(\xi \right)\\ \tilde{W}\left(\xi \right)& ={\tilde{W}}^{†}\left(-\xi \right)=-\mathrm{conj}{\tilde{W}}^{†}\left(\xi \right),\end{array}$$

$$\tilde{\mathsf{L}}\left(\xi \right)={\int }_{{\Omega }_{\mathrm{p}}}\tilde{\mathsf{L}}\left(x\right)e^{i\xi x}\mathrm{d}x.$$

$$\left(\begin{array}[]{c}\nabla\cdot\bm{\sigma}-s^{2}\rho\mathbf{u}\\ \nabla\cdot\mathbf{D}\end{array}\right)+\left(\begin{array}[]{c}\mathbf{f}\\ -q\end{array}\right)=\left(\begin{array}[]{c}\mathbf{0}\\ 0\end{array}\right).$$

$$\left(\begin{array}[]{c}\mathbf{G}_{11}^{\mathsf{T}}\cdot\left(\nabla\cdot\bm{\sigma}\right)+\mathbf{G}_{21}\nabla\cdot\mathbf{D}\\ \mathbf{G}_{12}\cdot\left(\nabla\cdot\bm{\sigma}\right)+\mathsf{G}_{22}\nabla\cdot\mathbf{D}\end{array}\right)-\left(\begin{array}[]{c}s^{2}\rho\mathbf{G}_{11}^{\mathsf{T}}\cdot\mathbf{u}\\ s^{2}\rho\mathbf{G}_{12}\cdot\mathbf{u}\end{array}\right)+\left(\begin{array}[]{c}\mathbf{G}_{11}^{\mathsf{T}}\cdot\mathbf{f}-\mathbf{G}_{21}q\\ \mathbf{G}_{12}\cdot\mathbf{f}-G_{22}q\end{array}\right)=\left(\begin{array}[]{c}\mathbf{0}\\ 0\end{array}\right),$$

Full text

Keywords: dynamic homogenization, piezoelectric composite, Bloch Floquet waves, metamaterials, effective properties, constitutive relations, wave propagation, Willis coupling

Faculty of Mechanical Engineering, Technion–Israel Institute of Technology, Haifa 32000, Israel

1 Introduction

The effective or macroscopic properties of materials are modeled by the coupling parameters between physical fields in germane constitutive equations [1]; extraordinary properties are engineered by cleverly designing the microstructure of artificial materials. Such metamaterials were developed in optics, acoustics and mechanics for various objectives [2, 3, 4, 5].

One of the grand challenges in metamaterial design is to obtain control over traversing waves [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. Mechanical waves are governed by the balance of momentum; at the microscale, the momentum is coupled only to the material point velocity by the mass density. Willis discovered that in elastic materials with a specific microstructure, the macroscopic momentum is coupled also to the strain by the now termed Willis coupling [17]. This coupling thus offers a designable degree of freedom to manipulate waves.

A series of theoretical studies were carried out to characterize Willis coupling and understand its physical origins [18, 19, 20, 21, 22, 23, 24, 25, 26]. Guided by accompanying predictions that Willis coupling is connected with unusual phenomena such as asymmetric reflections and unidirectional transmission, recent experimental realizations of Willis metamaterials that demonstrate these phenomena were reported [27, 28, 29, 30, 31, 32, 33, 34, 35]. To date, these investigations were limited to metamaterials that are deformable only by mechanical forces.

Here, we consider constituents that additionally deform by non-mechanical stimuli, namely, piezoelectric materials responding to electric fields [36, 37]. We construct asymmetric patterns of such responsive materials, and show that their macroscopic momentum can additionally be coupled by design to the stimulus. We call this macroscopic property the electro-momentum coupling (Fig. 1). Akin to the intrinsic piezoelectric and engineered Willis couplings, this coupling appears in (meta)materials with no inversion symmetry. Beyond its theoretical significance, mechanical metamaterials designed with this property can actively manipulate waves by modulation of the external stimulus, contrary to typical metamaterials whose functionality is fixed, cf. Refs.  [38, 39, 40, 41]. To this end, it is required to carry out complementary studies on the connection between the electro-momentum coupling, scattering properties, and medium composition [42], similarly to the process that was required in employing Willis coupling for metamaterial design [28, 25, 34]; the present work opens the route for these studies.

The macroscopic properties of metamaterials are analytically calculated using homogenization or effective medium theories [43, 44, 45, 46, 47, 48, 49, 50, 51, 52]. Guided by the effective elastodynamic theory of Willis [53, 54, 55, 56, 57, 58, 59], we develop a homogenization method for piezoelectric metamaterials, whose application unveils the electro-momentum coupling. Our method is based on three elements. Firstly, it employs a unified framework we developed to account for the microscopic interactions between the mechanical and non-mechanical fields. Secondly, it uses an averaging scheme whose resultant effective fields identically satisfy the macroscopic governing equations. To this end, we have adapted the ensemble averaging approach of Willis [57] to the current setting. Lastly, it incorporates driving forces that render the effective properties unique, as firstly advocated by Fietz and Shvets [44], and later in Refs. [57, 23].

Before proceeding, a short discussion on the applicability of Willis homogenization scheme, and by transitivity our scheme, is in order. First, we note that the scheme is independent of any assumptions, and delivers effective properties that identically satisfy the field equations and boundary conditions, therefore considered exact [18, 26]; in fact, asymptotic homogenization schemes were shown to be approximations of Willis homogenization [20, 26]. Thus, for infinite periodic medium, the scheme reproduces precisely the its corresponding band diagram. Still, for the homogenized fields to serve as good approximations of the microscopic ones, certain homogenizability conditions should be satisfied [18]. These show that beyond the long-wavelength low-frequency limit there is an additional range in which dynamic homogenization is meaningful [50].* We further note that the new coupling terms emerging from this scheme were found essential for obtaining an effective description that satisfies fundamental principles such as causality [23], similarly to the* need for the analogous bianisotropic coupling in electromagnetics [60, 45].

The paper is structured as follows. In Sec. 2 we summarize the equations governing elastic waves in heterogeneous linearly responsive media, i.e., media that respond mechanically to non-mechanical loads. In Sec. 3 we develop our dynamic homogenization scheme for such media, and highlight the essential features of the resultant macroscopic properties. We apply our theory to longitudinal waves in piezoelectric layers assembled periodically with broken inversion symmetry in Sec. 4. A summary of our work concludes the paper in Sec. 5.

Before we proceed, we stress out that while this introduction was concerned with piezoelectric materials which respond to electric stimuli, our results immediately apply to piezomagnetic materials, owing to the mathematical similarity between their microscopic equations. We also expect our framework to benefit other types of elastic media that respond to non-mechanical stimuli, such as thermoelastic composites. These extensions are discussed in Appendix A

2 Equations of piezoelectric composites

Consider a linear material with mass density $\rho$ and stiffness tensor $\mathbf{C}$ occupying the volume $\Omega$, subjected to prescribed body force density $\mathbf{f}$ and inelastic strain $\bm{\eta}$. Physically, $\bm{\eta}$ can result from a plastic process or phase transformation, where the mathematical motivation to account for it will be explained in the sequel [44, 57]. The response of the material is governed by the balance of linear momentum $\mathbf{p}$

[TABLE]

where $\bm{\sigma}$ is the Cauchy stress second-order tensor and the superposed dot denotes time derivative. The material mechanically interacts with the electric displacement $\mathbf{D}$, and electric field $\mathbf{E}$, which satisfy

[TABLE]

where $q$ is a prescribed free charge density, acting as a source similarly to $\mathbf{f}$. Note that Eq. (2)2 is identically satisfied by setting $\mathbf{E}=-\nabla\phi$.

The constitutive relations between the fields are written in symbolic matrix form [61]

[TABLE]

where $\mathbf{u}$ is the displacement field, $\mathbf{A}$ and $\mathbf{B}$ are the dielectric and piezoelectric (tensorial) properties, and the transpose of $\mathbf{B}$ is defined by $B^{\mathsf{T}}_{ijk}=B_{kij}$. We assume the standard tensor symmetries, which in components read

[TABLE]

For later use, we denote the matrix (resp. column vector) in the right (resp. left) hand side of Eq. (3) by $\mathsf{L}$ (resp. $\mathsf{h}$).

We clarify that the elements in the symbolic matrices here and in what follows are differential operators and tensors of different order, and their product should be interpreted accordingly. For example, the product $\mathsf{L}_{11}\mathsf{m}_{1}$ represents double contraction, which in tensor notation is $\mathbf{C}:\bm{\eta}$, and in index notation is $C_{ijkl}\eta_{kl}$, while $\mathsf{L}_{12}\mathsf{b}_{2}$ is the single contraction $\mathbf{B}^{\mathsf{T}}\cdot\nabla\phi$, or $B^{\mathsf{T}}_{ijk}\phi_{,k}$. We further note that the symbolic matrix structure can be cast into standard matrix representation, using Voigt notation. Accordingly, $\mathbf{C},\mathbf{B}^{\mathsf{T}}$ and $\mathbf{B}$ are representable by $6\times 6,6\times 3$ and $3\times 6$ matrices, respectively, such that $\mathsf{L}$ is a symmetric $12\times 12$ matrix. (The zeros are $9\times 3$ and $3\times 9$ null matrices, and $\rho\mathbf{I}$ is a $3\times 3$ matrix.) The symmetric tensor $\bm{\sigma}$ is mapped to a $6\times 1$ column vector (and so is the symmetric part of $\nabla\mathbf{u}$), such that $\mathsf{h}$ is a $12\times 1$ column vector, and so on.

The prescribed boundary conditions are $\mathbf{u}=\mathbf{u}_{0}$ and $\phi=\phi_{0}$ over $\partial\Omega_{\mathsf{w}}\subset\partial\Omega$, and across the remaining boundary $\partial\Omega_{\mathsf{t}}=\partial\Omega\backslash\partial\Omega_{\mathsf{w}}$ are $\bm{\sigma}\cdot\mathbf{n}=\mathbf{t}_{0}$ and $\mathbf{D}\cdot\mathbf{n}=-w_{e}$, where $\mathbf{n}$ is the outward normal, $\mathbf{t}_{0}$ is the traction, and $w_{e}$ is the surface charge density.

When the medium is randomly heterogeneous, its properties $\rho,\mathbf{A},\mathbf{B}$ and $\mathbf{C}$ are functions of the position $\mathbf{x}$ and a parameter $y$ of a sample space $\Lambda$ with certain probability measure. Importantly, a periodic medium—the prevalent case of interest for metamaterials—can be analyzed as random, by considering different realizations of the composite generated by periodizing representative volume elements whose corner is a uniformly distributed random variable, and identifying $y$ with this variable. [57, 59].

Next, observe that ensemble averaging of Eqs. (1) and (2) over $y\in\Lambda$, denoted by $\left\langle\cdot\right\rangle$, provides

[TABLE]

where $\left\langle\mathbf{f}\right\rangle=\mathbf{f}$ and $\left\langle q\right\rangle=q$ since $\mathbf{f}$ and $q$ are sure (prescribed). Eq. (5) suggests the use of $\left\langle\bm{\sigma}\right\rangle,\left\langle\mathbf{D}\right\rangle$ and $\left\langle\mathbf{p}\right\rangle$ as effective fields that identically satisfy the governing equations. The effective properties are thus the quantities that relate these effective fields with $\left\langle\nabla\mathbf{u}\right\rangle,\left\langle\dot{\mathbf{u}}\right\rangle$ and $\left\langle\nabla\phi\right\rangle$, to form effective constitutive relations. Together with Eq. (5), they establish a meaningful description of the material when the ensemble averaged fields fluctuate slowly enough relatively to the scale of the microstructure; for a rigorous description of the applicability conditions for homogenization see Ref. [18]. In periodic media undergoing Bloch-Floquet waves, these ensemble averages reduce to volume averages over the periodic part of each field111The equivalence between ensemble averaging and volume averaging of the periodic part will be exemplified in Sec.  4 in the scalar case, without the loss of generality. [57, 18]. Qualitatively, the effective fields have the form of the curve in the right sketch of Fig. 2, after the fluctuations of the curve in the left sketch have been averaged out. The outstanding problem is to calculate the effective properties. Before we derive them, we can now provide a formal statement of our main result: homogenization shows that the effective constitutive relations are in the form

[TABLE]

with the electro-momentum coupling $\tilde{\mathbf{W}}^{\dagger}$ (Fig. 2), where $\left(\cdot\right)^{\dagger}$ denotes the adjoint operator with respect to the spatial variable. Thus, our homogenization process exposes effective couplings between $\left\langle\mathbf{D}\right\rangle$ and the average velocity, and between $\left\langle\nabla\phi\right\rangle$ and the average momentum. We denote the matrix of the effective properties in Eq. (15) by $\tilde{\mathsf{L}}$. The terms $\left\langle\nabla\mathbf{u}\right\rangle-\bm{\eta}$ and $\left\langle\dot{\mathbf{u}}\right\rangle$ that $\tilde{\mathsf{L}}$ operates on are independent, owing to $\bm{\eta}$, thus rendering $\tilde{\mathsf{L}}$ unique; otherwise, the fields $\left\langle\nabla\mathbf{u}\right\rangle$ and $\left\langle\dot{\mathbf{u}}\right\rangle$ are derived from the same potential, resulting with non-unique $\tilde{\mathsf{L}}$ [57], see, e.g., the example in Ref.  [62]. The calculation of $\tilde{\mathsf{L}}$ is detailed next.

3 Derivation of the effective properties

We adapt the ingenious approach of Willis [57] to responsive metamaterials as follows. Firstly, we cast the problem into matrix equations whose entries are tensors and other operators. Accordingly, Eqs. (1), (2)1 and (3) take the form

[TABLE]

where

[TABLE]

and we employed the Laplace transform to replace time derivatives by products with $s$; to reduce notation, we retain the same symbols for the temporal and transformed fields, bearing in mind they now represent functions of $s$.

We define the Green’s function $\mathsf{G}\left(\mathbf{x},\mathbf{x}^{\prime}\right)$ via

[TABLE]

where $\mathbf{I}$ is the second order identity tensor and $\delta$ is the Dirac delta; the entries $\mathsf{G}_{11}$, $\mathsf{G}_{12}$ and $\mathsf{G}_{21}$, and $\mathsf{G}_{22}$ are second order tensor-, vector-, and scalar-valued functions, respectively, with the homogeneous boundary conditions

[TABLE]

where

[TABLE]

We derive next a useful expression for $\mathsf{w}$, by right multiplying the transpose of Eq. (18) by $\mathsf{w}\left(\mathbf{x}\right)$, subtracting it from the left product of Eq. (16) by $\mathsf{G}^{\mathsf{T}}\left(\mathbf{x},\mathbf{x}^{\prime}\right)$, and integrating the difference over the volume $\mathbf{x}^{\prime}\in\Omega$. The result is

[TABLE]

The development of Eq. (20) employs the divergence theorem, boundary conditions for $\mathsf{G}$ and symmetries of $\mathbf{A},\mathbf{B}$ and $\mathbf{C}$, and its detailed derivation using index and tensor notations is provided in Appendix B. Next, we manipulate the last integral in Eq. (20) as follows. First, we formally extend the integral domain from $\partial\Omega_{\mathsf{w}}$ to the whole boundary $\partial\Omega$ using the fact that the integrand vanishes over $\partial\Omega_{\mathsf{t}}$, owing to homogeneous boundary conditions that $\mathsf{G}$ satisfies over $\partial\Omega_{\mathsf{t}}$. Next, we can replace $\mathsf{w}$ with $\left\langle\mathsf{w}\right\rangle$, since it is sure on $\partial\Omega_{\mathsf{w}}$ and the integrand vanishes where it is not. By transforming the surface integral it back to the volume via the divergence theorem, and with the aid of Eq. (18), we can rewrite $\mathsf{w}$ as

[TABLE]

Using the ensemble average of Eq. (21) and the fact that $\mathsf{m},\mathsf{f},\mathbf{t}_{0}$ and $w_{e}$ are sure, we have that (see Appendix B for more details)

[TABLE]

Finally, we substitute Eq. (22) into Eq. (17) and ensemble average the result to obtain the effective constitutive relations

[TABLE]

with

[TABLE]

Eq. (24) generalizes the result of Willis [57] to metamaterials that interact with non-mechanical fields. The components of $\tilde{\mathsf{L}}$ define the effective properties in Eq. (15), which are non-local operators in space and time. The non-zero adjoint terms $\tilde{\mathsf{L}}_{23}$ and $\tilde{\mathsf{L}}_{32}$ reveal the coupling between $\left\langle\mathbf{D}\right\rangle$ and $\left\langle\dot{\mathbf{u}}\right\rangle$, and between $\left\langle\mathbf{p}\right\rangle$ and $\left\langle\nabla\phi\right\rangle$, respectively, denoted by $\tilde{\mathbf{W}}$. In the case that $\mathbf{D}$ and $\nabla\phi$ are vectors, the coupling $\tilde{\mathbf{W}}$ is a second order tensor. Owing to the symmetries of $\mathbf{A},\mathbf{B}$ and $\mathbf{C}$, the operator $\mathsf{L}$ is symmetric, so that $\mathsf{G}$ self-adjoint with the usual symmetries of Green functions, and hence $\tilde{\mathsf{L}}$ is self-adjoint as well, justifying the notation $\left(\cdot\right)^{\dagger}$ for $\tilde{\mathsf{L}}_{31}$ and $\tilde{\mathsf{L}}_{32}$.

4 Application to piezoelectric

layers

We exemplify the emergence of the electro-momentum coupling by calculating the effective properties of an infinite repetition of commercially available piezoelectric layers. Specifically, we will study two different periodic cells, namely, (i) all-piezoelectric cell made of PZT4-BaTiO3-PVDF layers, henceforth called composition 1, and (ii) one piezoelectric BaTiO3 layer between elastic Al2O3 layer and another elastic PMMA layer, henceforth called composition 2. The material properties of the comprising phases are given in Tab. 1, where the values for the piezoelectric materials correspond to the coefficients in the direction of the poling.

The periodic cell is denoted $\Omega_{\mathrm{p}}$, and its period is denoted $l$; in the calculations that follow we fix $l=3\,$mm. The layers are oriented such that the poling direction is along the direction of lamination, say $x$. The composite is driven by a body force density $f$ acting in the $x$ direction, along which axial inelastic strain $\eta$ is present; there are no free or surface charge sources in the problem ($w_{e}=q=0$). As a result, longitudinal waves propagate in the $x$ direction, such that the problem is one-dimensional and the pertinent fields can be treated as scalars. The objective is to obtain the macroscopic description of this problem by means of our homogenization scheme.

To this end, we analyze the periodic medium as random by treating the position of the period as a uniformly distributed random variable, with uniform probability density $l^{-1}$ over $\Omega_{\mathrm{p}}$. Accordingly, any $l$-periodic function $\zeta_{y}\left(x\right)$ in realization $y\in\Omega_{\mathrm{p}}$ is $\zeta_{0}\left(x-y\right)$; its ensemble average

[TABLE]

is independent of $x$, and equals the spatial average in any realization.

We examine first the governing equations in realization $y=0$. In the absence of free charge, Gauss law within any layer reads

[TABLE]

where superscript $n=a,b$ and $c$ denotes values in the first, middle, and third layer, respectively. Note that $D$ is continuous across the layers, and since there are no electrodes and surface charge, we have that $D=0$ everywhere. Further, Eq.  (26) together with the constitutive relation

[TABLE]

implies that

[TABLE]

The second governing equation is the equation of motion, namely,

[TABLE]

and we recall that the Laplace transform has been used, as in Sec.  3. Substituting in the constitutive relation for the stress

[TABLE]

and Eq.  (28) yields

[TABLE]

where the superscript notation was suppressed for brevity. Since the homogeneous equation derived from Eq. (31) has periodic coefficients, its Green function is constructed using Bloch-Floquet solutions. The Green function and its ensemble average are thus

[TABLE]

where $\left\langle\mathsf{G}\right\rangle$ is only a function of $x-x^{\prime}$, $V$ is

[TABLE]

and $u_{\mathrm{p}}^{\pm}$ are $l$-periodic functions whose standard (and tedious) calculation is detailed in Appendix C.

As the simplicity of the problem enabled a solution via a single Green function, a simpler equivalent to Eq. (22) for $u$ follows, namely,

[TABLE]

where $\check{C}=C+\dfrac{B^{2}}{A}$. Observing that in the prescribed settings $D=0$ in any realization, we obtain the remaining fields $\phi\left(x,y\right)$ from Eq. (3) in terms of $u\left(x,y\right),A\left(x,y\right)$ and $B\left(x,y\right)$.

Finally, the equivalents of Eqs. (23)-(24) are derived directly by substituting $u\left(x,y\right)$ and $\phi\left(x,y\right)$ into Eq. (3), ensemble averaging, and identifying the terms that multiply $\left\langle u_{,x^{\prime}}\right\rangle-\eta,\left\langle\phi_{,x^{\prime}}\right\rangle$ and $\left\langle su\right\rangle$ as the effective properties according to Eq. (15). The price for using a single Green function in the absence of charge is one degree of freedom in calculating $\tilde{\mathsf{L}}$, which we eliminate by enforcing $\tilde{B}^{\mathsf{T}}=\tilde{B}$. Together with this choice, the resultant equations determine the components of $\tilde{\mathsf{L}}$. The explicit expressions are provided in Appendix D, and imply that

[TABLE]

where $\mathrm{conj}$ denotes complex conjugate, and $\xi$ is the Fourier transform variable, i.e.,

[TABLE]

These relations also imply that Re$\,\tilde{W}$ is an odd function of $\xi$, while Im$\,\tilde{W}$ is an even function, similarly to the functional form of $\tilde{S}$ [23]. Eq. (35) consolidates the notation used by Willis [56, 57, 58, 59], who interprets $\tilde{S}$ and $\tilde{S}^{\dagger}$ as formal adjoints with respect to the spatial variable, and other notations in the literature [63, 64, 23], which use $-\mathrm{conj}\,\tilde{S}$ instead of $\tilde{S}^{\dagger}$.

The numerical results in Fig. 3 display the Fourier transform $\tilde{W}\left(\xi\right)$ and its adjoint versus frequency $\frac{\omega}{2\pi}$, where $s=-i\omega$, across the first pass band. The calculations are done for the microstructure $l^{\left(a\right)}=1\,$mm (thickness of the first layer in the periodic cell) $,l^{\left(b\right)}=1.4\,$mm (middle layer), and $l^{\left(c\right)}=0.6\,$mm (last layer). Panel 3a shows the long-wavelength limit $\xi l=0$, which contains the homogenization limit [65, 50], where in Panel 3b we examine $\xi l=2$. For brevity, we omit the plots for the remaining effective properties, noting they are frequency-dependent, real at $\xi l=0$ (except $\tilde{S}$ and $\tilde{S}^{\dagger}$, which are pure imaginary), and generally complex for $\xi l>0$.

Evidently, the contrast between the piezoelectric coefficients of the layers greatly affects the magnitude of $\tilde{W}$. For example, the maximal value of Im$\,\tilde{W}$ at $\xi l=0$ for composition 1 is 21$\>\mu$Cs$\mathrm{m}^{-3}$, while for composition 2 it is over 900$\>\mu$Cs$\mathrm{m}^{-3}$.

Both compositions share the following notable features. Firstly, the electro-momentum coupling $\tilde{W}$ vanishes in the quasi-static limit $\omega\rightarrow 0$, as it should. In the long-wavelength limit and $\omega>0$, it attains non-zero pure imaginary values, as Eq. (35) implies. Above the long-wavelength limit, the numerical results confirm that the electro-momentum is complex such that Im$\tilde{\,W}=$Im$\,\tilde{W}^{\dagger}$ and Re$\tilde{\,W}=-$Re$\,\tilde{W}^{\dagger}$. Our observations conform with the insights of Sieck et al. [23] on the microstructure-induced Willis coupling: the imaginary part originates from broken inversion symmetry, hence appears even in the long-wavelength limit, while mesoscale effects of multiple scattering create its real part, hence appear at $\xi l>0$. This is further demonstrated in Panel  3c, where we evaluate $\tilde{W}$ and $\tilde{W}^{\dagger}$ at $\xi l=2$ in a system that is symmetric under inversion by setting $l^{\left(b\right)}=0$ in composition 1. Indeed, the coupling is pure real in this setting, and satisfies satisfies $\tilde{W}\left(\xi\right)=-\tilde{W}^{\dagger}\left(\xi\right)$.

Fig. 4a directly evaluates the dependency of $\tilde{W}$ on the microstructure, by plotting it against $l^{(b)}$ while setting $l^{(a)}=l^{(c)}$, at the representative frequency 0.1 MHz. The couplings vanish for $l^{(b)}=0$ and 3 mm, as it should for microstructures with inversion symmetry. Notably, the change from $l^{(b)}=0$ to $l^{(b)}=0^{+}$ is discontinuous for composition 2, as it reflects a shift from an elastic composition to a piezoelectric one. Interestingly, the coupling of composition 2 vanishes also at $l^{(b)}=0.79\,$mm, where the coupling of composition 1 exhibits singularity in the vicinity of $l^{(b)}=0.87\,$mm. These values are functions of the frequency; for example, at 0.05 MHz we calculated lower values, namely, $l^{(b)}=0.77\,$mm and $l^{(b)}=0.84\,$mm, respectively.

Lastly, in Fig. 4b we evaluate the frequency versus the momentum ratios between $\left|\tilde{W}^{\dagger}\left\langle\nabla\phi\right\rangle\right|,\left|\tilde{S}^{\dagger}\left\langle\nabla u\right\rangle\right|$ and $\left|\tilde{\rho}\left\langle\dot{u}\right\rangle\right|$ at $\xi l=0$, for composition 1 with $l^{(b)}=0.87\,$mm (black curves), and composition 2 (blue curves) with the microstructure studied in Fig. 3. Composition 1 exhibits a singularity at 1 MHz, conforming with the singularity observed in Fig. 4a. In the vicinity of 1 MHz, the growth of $\left|\tilde{W}^{\dagger}\left\langle\nabla\phi\right\rangle\right|$ is much faster than $\left|\tilde{S}^{\dagger}\left\langle\nabla u\right\rangle\right|$. Composition 2 demonstrates lower ratios, as expected from a metamaterial with an arbitrarily chosen microstructure.

5 Summary and discussion

We have developed an exact source-driven homogenization scheme for responsive metamaterials, based on the ensemble averaging approach of Willis and applied it in a numerical example considering periodic repetitions of commercially available piezoelectric layers. The scheme reveals that the effective non-mechanical and momentum-velocity fields can be coupled by properly designing the microstructure and composition of the medium. We conjecture that the corresponding couplings are necessary for obtaining an effective description that respect fundamental principles such as causality, similarly to the need for Willis coupling in elastodynamics* [23] and bianisotropic coupling in electromagnetics [60, 45].* The new couplings reflect energy conversion between electrical and mechanical energy in a distinct way than it occurs at the microscale. Therefore, they capture a new mechanism that can be employed for rectifying mechanical waves by modulation of external stimuli. We expect that future studies will show that extraordinary wave response such as asymmetric reflections and unidirectional transmission are modeled by the new couplings [42], analogously to Willis coupling   [28, 32, 33, 34]. These kind of results are to be guided by theoretical analyses based upon proper homogenization schemes [25]. The present work opens the route for such studies on the properties of the electro-momentum coupling and its experimental realizations for metamaterial design.

acknowledgements

We thank anonymous reviewers for their constructive comments that helped us improve this paper. We acknowledge the support of the Israel Science Foundation (Grant 1912/15), United States-Israel Binational Science Foundation (Grant 2014358), and MOST (Grant 880011). Science and Technology (Grant no. 880011).

Appendix A. Implications on other stimuli-responsive

media

The structure we developed and its resultant effective description immediately apply for piezomagnetic media, by observing the following mathematical connections. The fields $\mathbf{D}$ and $\mathbf{E}$ are equivalent to the magnetic induction and magnetic field (usually denoted by $\mathbf{B}$ and $\mathbf{H}$, respectively), since they both satisfy identical differential equations. The latter two fields are constitutively related by the second-order permeability tensor (usually denoted by $\mu$), and their coupling with the stress and strain is captured by the piezomagnetic third-order tensor [66]. Thereby, homogenization of piezomagnetic composites fits exactly into our scheme, which predicts new macroscopic second-order tensors that couple the magnetic induction and the velocity, and the momentum and the magnetic field.

Implications of our scheme to thermoelasticity are less immediate and require a separate treatment. However, at the very basic level, we can draw analogies between $\mathbf{D}$ and the increase of entropy per unit volume with respect to a reference state (denoted $\vartheta$) , and between $\mathbf{E}$ and the change in temperature with respect to some base temperature (denoted $\theta$). The fields $\vartheta$ and $\theta$ are scalar fields that are constitutively coupled through the constant of specific heat. The microscopic cross-coupling between $\vartheta$ and $\theta$ to the stress and strain is captured by the thermal expansion second-order tensor, see Eq.  (2.24) in [66]. We now proceed to the differential equations that govern thermoelasticity. While $\mathbf{D}$ and $\mathbf{B}$ are subjected to the same differential equation, there is no such spatial constraint on $\vartheta$; in some approximated formulations of thermoelasticity, the equations of heat conduction and energy balance can be combined to form with the equation of motion two coupled field equations for the temperature and displacement, see Eqs.  (2.10-23) or Eqs.  (9.6-7)-(9.6-8) in [67]. We conjecture that homogenizing this system will expose new macroscopic couplings between $\vartheta$ and the velocity, and between the momentum and $\theta$.

Appendix B. Detailed derivation of $\tilde{\mathsf{L}}$

in index and tensor notation

Eq. (16) in tensor notation is

[TABLE]

The left product of Eq. (B.1) with \mathsf{G}^{\mathsf{T}}\text{\left(\mathbf{x},\mathbf{x}^{\prime}\right)} is

[TABLE]

or, in index notation,

[TABLE]

The quantities $\mathsf{G}_{11}$, $\mathsf{G}_{12}$ and $\mathsf{G}_{21}$, and $\mathsf{G}_{22}$ are second order tensor-, vector-, and scalar-valued functions, respectively, and hence the symbolic $2\times 2$ matrix $\mathsf{G}$ is representable by a $4\times 4$ matrix. ($\mathsf{G}_{11}$, $\mathsf{G}_{12}$ and $\mathsf{G}_{21}$, and $\mathsf{G}_{22}$ are represented by $3\times 3,3\times 1,1\times 3$ and $1\times 1$ blocks, respectively.) The equation that defines $\mathsf{G}\left(\mathbf{x},\mathbf{x}^{\prime}\right)$ is

[TABLE]

where, in tensor notation, the elements of $\delta\mathsf{I}$ are

[TABLE]

and the elements of $\mathsf{D}^{\mathsf{T}}\mathsf{L}\mathsf{B}\mathsf{G}$ are

[TABLE]

or, in index notation,

[TABLE]

We emphasize that the convention contraction is different when $\mathsf{G}$ is involved. For example, $\mathbf{C}:\nabla\mathbf{G}_{11}$ is $C_{ijkl}G_{11kp,l}$ and $\mathbf{B}^{\mathsf{T}}\cdot\nabla\mathbf{G}_{21}$ is $B_{ijm}^{\mathsf{T}}G_{21p,m}$. The standard convention (contraction with the first two indices of the right tensor) can be recovered observing that

[TABLE]

and thus $\mathbf{C}:\nabla\mathbf{G}_{11}=\nabla\mathbf{G}_{11}^{\mathsf{T}}:\mathbf{C}.$ Similarly,

[TABLE]

and thus $\mathbf{B}^{\mathsf{T}}\cdot\nabla\mathbf{G}_{21}=\nabla\mathbf{G}_{21}\cdot\mathbf{B}$. Using such identities, the standard convention for divergence operation applies, i.e., acting on the last free index.

Right multiplying the transpose of Eq. (B.4) with $\mathsf{w}\left(\mathbf{x}\right)$ provides

[TABLE]

Subtracting Eq. (B.18) from (B.3) yields

[TABLE]

Eq. (B.23) is simplified using the following relations

[TABLE]

which, upon integration over the volume $\mathbf{x}\in\Omega$ and application of divergence theorem, then reads

[TABLE]

or, in tensor notation,

[TABLE]

where the tensor product $\otimes$ between vectors $\mathbf{a}$ and $\mathbf{b}$ is defined by the action on a third vector $\mathbf{c}$, namely [68],

[TABLE]

implying that $\left(\mathbf{a}\otimes\mathbf{b}\right)_{ij}=a_{i}b_{j}$ and $\left(\mathbf{T}\otimes\mathbf{b}\right)_{ijk}=T_{ij}b_{k}$ for second-order tensors $\mathbf{T}$. Accordingly,

[TABLE]

which, by defining

[TABLE]

allows us to write Eq. (B.51) in the symbolic matrix form

[TABLE]

We analyze next the last integral. First, we formally extend its domain from $\partial\Omega_{\mathsf{w}}$ to the whole boundary $\partial\Omega$ using the fact that the integrand vanishes over $\partial\Omega_{\mathsf{t}}$, owing to homogeneous boundary conditions

[TABLE]

that $\mathsf{G}$ satisfies over $\partial\Omega_{\mathsf{t}}$. Next, we can replace $\mathsf{w}$ with $\left\langle\mathsf{w}\right\rangle$, since it is sure on $\partial\Omega_{\mathsf{w}}$ and the integrand vanishes where it is not. By transforming the surface integral back to the volume via the divergence theorem and employing Eq. (18), we rewrite the second line as

[TABLE]

Now, since $\mathsf{f}$ is sure over $\Omega$ and $\mathsf{N}^{\mathsf{T}}\mathsf{L}\left(\mathsf{B}\mathsf{w}-\mathsf{m}\right)$ is sure over $\partial\Omega_{\mathsf{t}}$ (it is the surface charge and traction), we can combine the ensemble average of Eqs.  (B.55) and (B.57) to obtain

[TABLE]

and hence

[TABLE]

Finally, the above manipulations allow us to write $\mathsf{w}\left(\mathbf{x}^{\prime}\right)$ as

[TABLE]

In essence, Eq.  (B.60) is the generalization of Eq.  (3.14) by Willis [57] to piezoelectric media, and as such, relies upon similar reasoning in its derivation. The effective operator $\tilde{\mathsf{L}}$ in our settings is obtained by substituting Eq. (B.60) into $\mathsf{h}\left(\mathbf{x}^{\prime}\right)=\mathsf{L}\left(\mathsf{B}\mathsf{w}-\mathsf{m}\right)\left(\mathbf{x}^{\prime}\right)$, namely,

[TABLE]

Finally, ensemble averaging Eq. (B.61) and comparing it with $\left\langle\mathsf{h}\right\rangle=\tilde{\mathsf{L}}\left(\left\langle\mathsf{B}\mathsf{w}\right\rangle-\mathsf{m}\right)$ delivers the following expression for $\tilde{\mathsf{L}}$

[TABLE]

Appendix C. Construction of the Green function

using Floquet solutions

Recall that in realization $y=0$, the governing equations are combined to obtain in each layer

[TABLE]

Since the homogeneous equation derived from Eq. (C.1) has periodic coefficients, its Green function is constructed using Bloch-Floquet solutions, namely [69],

[TABLE]

where

[TABLE]

Before we proceed to calculate the periodic parts $u_{\mathrm{p}}^{\pm}\left(x\right)$, it is worth showing that ensemble averaging of Bloch-Floquet functions reduces to volume averaging over their periodic part [57, 50, 18]. To this end, observe that the solution, say $\mathcal{U}$, to a differential equation with $l$-periodic coefficients and Bloch wavenumber $k_{B}$ governing realization $y$ is of the form

[TABLE]

We now recall that $\mathcal{U}_{y}\left(x\right)$ in realization $y$ is related to realization 0 via (Fig.  5)

[TABLE]

and therefore ensemble averaging over $\Omega_{\mathrm{p}}$ amount to

[TABLE]

thereby equals to volume averaging of the period part. It is straightforward to extend these arguments to the three-dimensional case.

We return to calculate $u_{\mathrm{p}}^{\pm}\left(x\right)$ by employing a standard transfer matrix approach [70, 71], as follows. The phase-wise solution of the Eq. (C.1) is

[TABLE]

where $\alpha^{\left(n\right)}$ and $\beta^{\left(n\right)}$ are integration constants,

[TABLE]

and to abbreviate notation, we momentarily suppress the subscript [math]. The integration constants $\alpha^{\left(n\right)}$ and $\beta^{\left(n\right)}$ are determined from the continuity and Bloch-Floquet conditions on $u(x)$ and the stress $\sigma(x)$, which are compactly written using the state vectors $\mathsf{s}^{\left(n\right)}\left(x\right)$ and transfer matrices $\mathsf{T}^{\left(n\right)}$

[TABLE]

Thus, the state vectors at the ends of each layer are related via

[TABLE]

where continuity conditions at the interfaces are

[TABLE]

which are combined to obtain

[TABLE]

The latter quantities are also related via the (Bloch) Floquet quasi-periodicity condition

[TABLE]

Eqs. (C.14)-(C.15) deliver together the eigenvalue problem

[TABLE]

The standard condition for non-trivial solutions provides the dispersion relation

[TABLE]

which, upon substitution back into the foregoing equations, provides the eigenmodes as functions of $k_{B}$, and specifically $u^{\pm}=u_{\mathrm{p}}^{\pm}\left(x\right)e^{\pm ik_{B}x}$.

In the calculations to follow is more convenient to use the Fourier expansion of of $u_{\mathrm{p}}^{+}$ and $u_{\mathrm{p}}^{-}$. Accordingly, in terms of the Fourier coefficients

[TABLE]

the ensemble average of $G$

[TABLE]

namely,

[TABLE]

confirming it depends on $x$ and $x^{\prime}$ solely via $x-x^{\prime}$. Its Fourier transform

[TABLE]

reads

[TABLE]

We proceed to calculate the rest of the terms in the Fourier transform of $\tilde{\mathsf{L}}.$ To this end, we define the Fourier coefficients

[TABLE]

and obtain

[TABLE]

where $\varsigma_{\mathrm{\left(1\right)}}$ and $\varsigma_{\mathrm{\left(2\right)}}$ denote two different properties from the set $\left\{C,\check{C},\dfrac{B}{A},\dfrac{B^{2}}{A}\right\}$, and

[TABLE]

Note that

[TABLE]

The process is exemplified using the calculation of $\left\langle C(x)\,G{{}_{,x}}\right\rangle(\xi)$. Firstly, note that in realization $y=0$ is

[TABLE]

The translated expression is thus

[TABLE]

It follows that

[TABLE]

and its Fourier transform is

[TABLE]

Appendix D. Explicit expressions for the effective

properties

Recall that in the present problem

[TABLE]

and hence

[TABLE]

In terms of the effective properties, we we also have that

[TABLE]

By substituting

[TABLE]

into Eq. (D.2) and comparing with Eq. (D.3) we obtain

[TABLE]

from $\left\langle\sigma\right\rangle$; from $\left\langle\phi_{,x}\right\rangle$ we find

[TABLE]

and finally from $\left\langle p\right\rangle$

[TABLE]

As we pointed out in the body of the paper

[TABLE]

Also note that $\dfrac{\tilde{B}^{\mathsf{T}}\tilde{W}}{\tilde{A}}\left(\xi\right)=\dfrac{\tilde{B}\tilde{W}^{\dagger}}{\tilde{A}}\left(-\xi\right)$. As mentioned, the drawback for using a single Green function in the absence of charge is one degree of freedom in calculating $\tilde{\mathsf{L}}$, which we eliminate by enforcing $\tilde{B}^{\mathsf{T}}=\tilde{B}$. The latter condition, together with Eqs. (D.5)-(D.7) deliver the effective properties $\tilde{C},\tilde{A},\tilde{\rho},\tilde{B}$ and $\tilde{W},\tilde{W}^{\dagger}$. These satisfy

[TABLE]

where relation between $\tilde{W}$ and $\tilde{W}^{\dagger}$ is similar to the relation between $\tilde{S}$ and $\tilde{S}^{\dagger}$.