## Abstract

Materials that require coupling between the stress–strain and momentum–velocity constitutive relations were first proposed by Willis (Willis 1981 *Wave Motion* **3**, 1–11. (doi:10.1016/0165-2125(81)90008-1)) and are now known as elastic materials of the Willis type, or simply Willis materials. As coupling between these two constitutive equations is a generalization of standard elastodynamic theory, restrictions on the physically admissible material properties for Willis materials should be similarly generalized. This paper derives restrictions imposed on the material properties of Willis materials when they are assumed to be reciprocal, passive and causal. Considerations of causality and low-order dispersion suggest an alternative formulation of the standard Willis equations. The alternative formulation provides improved insight into the subwavelength physical behaviour leading to Willis material properties and is amenable to time-domain analyses. Finally, the results initially obtained for a generally elastic material are specialized to the acoustic limit.

## 1. Introduction

Metamaterials are engineered composite materials that derive their interesting properties from hidden degrees of freedom [1–3]. One of the most interesting aspects of metamaterials is that their properties often surpass those of naturally occurring and man-made composite materials. For example, materials have been made with negative dynamic mass density [4], negative dynamic modulus [5] and negative phase speed [6]. Models of materials displaying generalized elastodynamic response which couples stress and momentum density to both strain and particle velocity have also been explored [7–10]. The latter class of materials is known as Willis materials, and their description requires the definition of an additional material parameter in the momentum and dynamic stress relationships [11]. In attempting to design materials with extreme properties, such as Willis materials, it is of critical importance to know the constraints placed on the material properties by fundamental physical laws. These constraints provide the scientist with an understanding of what material properties are physically meaningful and provides a deeper understanding of the material. It also gives both theoreticians and experimentalists a means to check that predicted or measured effective material properties do not violate underlying physical principles.

Restrictions on the range of material properties are derived from assumptions about the material itself. For example, standard elastic materials are assumed to be passive and causal. As passive elastic and acoustic metamaterials are composed of passive materials, they must also be passive and causal. Similarly, if all of the constituent materials are reciprocal, reciprocity must also be enforced [12]. The restrictions imposed on general Willis-type material properties due to passivity have been derived by Srivastava & Nemat-Nasser [9] and Srivastava [13]. The restrictions imposed by causality have been derived for general Willis materials [13], but the implications of these restrictions were not fully explored. The consequences of restrictions placed on Willis coupling parameters due to passivity, reciprocity and causality are the primary contribution of this work.

The purpose of this paper is to derive general restrictions on the material properties of a linear Willis material, which are described in §2. Reciprocity and passivity in Willis materials are discussed in §3 and then combined. Reciprocity is shown to impose a symmetry on the material properties corresponding to the major symmetries of the stiffness tensor. Passivity in a reciprocal material is shown to restrict the sign of the imaginary part of the density and stiffness tensors and the real part of the Willis coupling parameters. Section 4 investigates the requirements of causality in a local medium using arbitrary history functions. Considerations of low-order dispersion in the light of requirements imposed by causality suggest an alternative representation for Willis coupling which is more appropriate for time-domain analyses. Sections 5a,b present requirements on dispersion relations in two separate limiting cases, and §5c specializes in the restrictions on the material properties to the acoustic case for application to heterogeneous fluids. Section 6 then summarizes the work and discusses its relevance to the fields of acoustic and elastic metamaterials.

## 2. Constitutive equations and minor symmetries

A linear Willis material is a composite material which, when homogenized, requires the stress–strain and the momentum–velocity constitutive equations to be coupled. This non-intuitive coupling arises from both non-local effects and from material asymmetries [7–10]. Using standard index notation and assuming a time-harmonic motion (e^{−iωt} time convention), the generalized elastodynamic constitutive equations take the form [7,8,10]
*σ*_{ij} is the Cauchy stress (hereafter called the stress), *C*_{ijkl} is the stiffness tensor, *u*_{i,j} is the displacement gradient, *μ*_{i} is the momentum density (hereafter called the momentum) and *ρ*_{ij} is the mass density tensor. The third-order tensors *S*_{ijk} and *f*(*ω*,** k**). Notice that equations (2.1) are written in the more general form using the displacement gradient rather than the strain [8].

Willis coupling encompasses all linear phenomena which generate stress from a net translation and linear momentum from local strain [11,14]. For example, in reciprocal media this coupling has been demonstrated to arise from asymmetry of a representative volume element (RVE) [10,14,15] and through non-local phenomena, such as lattice effects [16,17]. Willis coupling has also been demonstrated to arise in non-reciprocal media, although not described as such [18]. Willis coupling is mathematically analogous to bianisotropic materials in optics, and so chiral materials may also be considered a subset of general Willis materials [19]. For the sake of simplicity, the following derivation will assume that the Willis material is irrotational (so the displacement gradient *u*_{i,j} may be replaced by the symmetric strain tensor *ε*_{ij}) and only accounts for local interactions, unless explicitly stated otherwise.

For this set of assumptions, various symmetries of the material properties follow from standard elasticity. Assuming no external body forces, the stress tensor is symmetric [20], so the stiffness and the stress–velocity coupling tensor satisfy the symmetry conditions

## 3. Reciprocity and passivity

Reciprocity is qualitatively described as being able to exchange a source and a receiver without changing the measured signal, and is a common characteristic of linear elastic and acoustic media [12]. A passive medium cannot generate energy. This is another way of stating that a passive material cannot amplify mechanical disturbances such as elastic waves as they propagate through it. Restrictions on the Willis material properties imposed by assuming that a medium reciprocal, passive or both passive and reciprocal will be derived in this section. It is noted that the constraints due to reciprocity and passivity have been derived elsewhere using an alternative approach [9].

### (a) Reciprocity

As mentioned above, Willis material properties are a result of the homogenization of inhomogeneous elastic composites. As linear elastic materials are reciprocal, it follows that the emergent Willis material properties should also exhibit reciprocity [12]. The restrictions placed on the material properties of a linear Willis material, by assuming that they are reciprocal materials, are derived in this section. The present analysis differs from many analyses associated with reciprocity, in that we assume that the material is reciprocal to derive material symmetries, whereas previous work assumed material symmetries to derive reciprocity [12].

The dynamic equation for a time-harmonic linear elastic system may be written as the real part of
** f** e

^{−iωt}is a time-harmonic body force density,

*σ*

_{ij,j}represents the divergence of the stress field and the over-dot denotes a time derivative. Note that for time-harmonic systems

*u*

_{i}. The common e

^{−iωt}may be neglected by noting that the stress, body force and momentum density can be complex-valued and represent the time-averaged fields.

Consider two elastodynamic source distributions in the same material, denoted A and B. These two source distributions are independent of each other, and arbitrary except that they satisfy the dynamic equations

Using the relations *F* may be expanded as
*M* may also be expanded as

As the choice of force source distributions was arbitrary, the strain and velocity fields must be considered independent fields. Therefore, the first and third terms in equation (3.9) must be individually equal to zero. Consequentially, the second term must also be zero, yielding

### (b) Passivity

A passive mechanical medium is one that does not supply any mechanical energy. It can only respond to externally provided stimulus. This is expressed in terms of energy with the statement that net energy flux *into* a region of a passive material must be equal to or greater than zero. As with reciprocal media, any heterogeneous material constructed of passive materials must also be passive. This is true for standard elastic materials and generalized elastodynamic media like Willis materials. The restrictions imposed by passivity on Willis materials have also been derived by other authors (see [9,13]), but they are provided in this section for completeness. The present derivation will follow the approach of Banerjee [22].

The net energy flux per period into a finite-sized region *Ω* with boundary *Γ* is the real part of the complex power [23], while the imaginary part is known as the reactive power, which does not contribute to the energy transmitted to *Ω*. The complex power in the regions *Ω*, denoted by *P*, is given by
*n*_{i} is the unit outward normal, and the equation has been simplified using the divergence theorem and the dynamic equation given in equation (3.1). Using the Willis relations, equation (3.11) may be expanded as

Through the careful choice of boundary conditions, the velocity and strain in the region *Ω* may be imposed independently, assuming that *Ω* is small enough relative to a wavelength that no significant wave motion may occur within the domain. Therefore, one may set the particle velocity *v*_{i} equal to zero to yield an inequality involving the stiffness, or one may set the strain equal to zero to yield an inequality involving the mass density, written respectively as

### (c) Reciprocal and passive media

If the medium is both reciprocal and passive, the major symmetries of the stiffness and density tensors can be exploited to simplify equations (3.13) to yield
*Y* ′≡ℜ{*Y* } and *Y* ′′≡ℑ{*Y* }. In words, equations (3.15) mean that for a passive material and using the e^{−iωt} time convention, the imaginary part of the stiffness tensor must be negative definite and the imaginary part of the density tensor must be positive definite. Equation (3.14) may also be simplified to give

The result that the coupling tensors must be purely imaginary in the lossless case is interesting, given the fact that complex coupling tensors have emerged from a derivation that only considers lossless media [10,15]. This is counter to the stiffness and density tensors which are purely real unless losses are present. For those properties, the imaginary component is non-zero because the material response lags the excitation in lossy media. On the other hand, the imaginary part of the coupling tensors have been associated with asymmetry at the inclusion scale [10,15], while the real part has been observed to be due to non-local coupling phenomena from ‘weak spatial dispersion effects associated with the finite phase velocity along the array’ [16]. This explains the apparent discrepancy between the present derivation, which does not account for such non-local phenomena, and the complex-valued predictions of the above-mentioned derivation which assumes a periodic lattice.

The bounds imposed by reciprocity may be derived for coupling tensors accounting also for non-local phenomena. As shown by Alù [16] and Sieck *et al*. [17], local coupling phenomena are accounted for by coupling tensors which are *even* in wavenumber which requires, for example, that *S*_{ijk}(*ω*,** k**)=

*S*

_{ijk}(

*ω*,−

**). On the other hand, first-order non-local coupling phenomena are accounted for by coupling tensors which are**

*k**odd*in wavenumber; for example,

*S*

_{ijk}(

*ω*,

**)=−**

*k**S*

_{ijk}(

*ω*,−

**). In general, the coupling tensors may be a combination of even and odd components (labelled by superscripts), such that**

*k*## 4. Causality

In general, the strain and velocity fields at a given time and location in a material depend upon the fields at other times (material memory) and locations (non-locality). By assuming that nothing in the future may affect the present state of the system (that is, the system is causal) one restricts the admissible functions that describe the material response [24,25]. Welters *et al.* showed that causality is guaranteed for linear and passive materials for electromagnetic wave propagation [26], which is similar to elastodynamic wave propagation in many respects. While causality is well known and often used in electromagnetism to provide restrictions on the spectral shape of material properties for special cases, the principle of causality has primarily been used in the fields of acoustics and elastodynamics for measurement purposes [27–29]. Therefore, we will start our discussion of the effects of causality in a linear, local Willis material by deriving the Kramers–Krönig relations.

For the present analysis, it is assumed that only local effects are significant, and that all of the results will hold point-wise throughout the material. Despite the fact that lattice effects are known to be important in Willis materials [16], an analysis on the effects of non-local causality is beyond the scope of this work. For this case, a constitutive equation for a linear material relating a fields *s*=0, and are purely real. Further, the history functions are assumed to go to zero as ^{iωt} (the complex conjugate of our time convention e^{−iωt}), and time averaging yields
*A*_{n}, as
*A*_{n}(*ω*)=*A*_{n}′(*ω*)+*iA*_{n}′′(*ω*), where *A*′_{n}(*ω*) and *A*′′_{n}(*ω*) are the real and imaginary parts of *A*_{n}(*ω*), respectively, defined as
*ω* and at *ω* plane. Thus, in the limit *A*′′_{n}(*ω*)=0. Then, as *ω* plane, Cauchy’s residue theorem allows one to express the material property, *A*_{n}, as
*A*′′_{n}(*ω*) to *ω*. Finally, as *A*′_{n}(*ω*) is even and *A*′′_{n}(*ω*) is odd with respect to *ω*, the integrals in equations (4.6) may be reduced to only positive values of *λ* to give

### (a) An informative limiting case

Equations (4.8) provide explicit relationships between the real and imaginary parts of the dispersive properties of a causal Willis material. Unfortunately, the generality of these expressions does not lend itself to easy implementation in determining general theoretical restrictions imposed on Willis material properties by assuming causality. An important question, then, is: How can the Kramers–Krönig relations be used to learn more about Willis materials? Some answers to that question can be obtained by first acknowledging assumptions underlying the implementation of equations (4.8) and then investigating limiting cases. In order to use these relationships, one must first acknowledge the assumptions implicit in this formulation. The primary assumption relevant to this study is that the history functions that are valid and associated material properties describe a material that is well represented as a continuum. That is to say, one must assume that all small-scale physics underlying the dispersive nature of the medium are accurately represented by the history functions at every material point. For example, one does not need to explicitly consider molecular dynamics associated with relaxation in a viscoelastic medium to model the causal stress–strain response of the medium using a stiffness history function. Under the continuum hypothesis one must therefore assume that the upper limit of the integrals in equations (4.8) is only infinite in the sense that it is sufficiently high to not restrict the accuracy of the relationships between real and imaginary parts of the material property.

Beyond these basic assumptions, it is also important to note that it is rarely, if ever, possible to know the real or imaginary parts of any material property at all frequencies. Numerous authors have addressed this difficulty by limiting their analysis to certain frequency ranges where simplifying assumptions may be made such that knowledge of the material response outside of that range is unnecessary to provide accurate relations between the real and imaginary parts of a material property. Common examples include investigation of regions of low loss near poles [24] and limiting analysis to regions of weak dispersion [27]. These specific cases and their implications on dispersive Willis material properties will be discussed in more detail in §5a,b.

One specific limiting case of interest is that of a passive, reciprocal and lossless Willis material in the long wavelength limit. This case is of specific interest to the field of acoustic and elastic metamaterials, which focuses on effective properties of a medium that result from deeply subwavelength microstructure. In the long wavelength limit, one can evaluate the Kramer–Krönig relations if the upper limit of equations (4.8) is set to a frequency that is simultaneously much greater than the frequency of interest and well below the lowest microstructural resonance frequency. This interpretation of the integral limits is permissible because of the continuum approximation [31]. In this case, the material properties are local in time and space, and all losses are negligible. While the upper frequency limit described here is difficult, if not impossible, to define from a mathematical standpoint, the value of this approximation to aid in physical reasoning based on a causality argument will become clear in the following paragraphs. For this case, equation (3.17) provides the passivity requirement that the imaginary parts of the stiffness and mass density tensors are zero and the real parts of the coupling tensors are also zero at all frequencies considered in the Kramer–Krönig relations. The requirements of causality summarized in equations (4.8) then simplify to
*et al*. [10] used Bloch wave analyses to determine the effective coupling tensors for a one-dimensional elastic lattice, Sieck *et al.* provided estimates for the analogous one-dimensional acoustic system [17], and Milton calculated coupling tensors for a two-dimensional periodic lattice [32]. In all of these works, complex valued coupling tensors have been predicted in the low-frequency limit when all constituent materials of the microstructure were assumed to be lossless. One thus concludes that while equations (2.1) provide an accurate and insightful representation the behaviour of materials with asymmetric microstructure (i.e. they properly describe the motion of the system), they are not a universally causal representation. While a causal representation may not be necessary for frequency domain analyses, effective medium theory that does not provide causal effective properties may pose problems for time-domain analyses and to interpretation of experimental data [27–29]. It is therefore important to investigate whether a causal representation can be found.

One common aspect of the studies mentioned above is that the lowest order dispersion of the imaginary parts of the coupling tensors is proportional to −*iω*. Willis also demonstrated this proportionality analytically using a Green’s function approach [33,34]. As time derivatives in the frequency domain may be substituted by multiplication by −*iω*, this form suggests a time derivative operation on the field variable may be hidden within the effective Willis coupling coefficient. This conclusion is supported by a simple model of a one-dimensional multi-degree of freedom ‘particle’ with inherently asymmetric microstructure, such as the one shown in appendix A. Analysis of that particle predicts that stress is related to the *acceleration* rather than the velocity, and the momentum density is related to the *strain rate* rather than the strain. Furthermore, Nassar *et al.* also suggested that the low-order dispersion predicted by their model of Willis coupling was actually a time derivative on the field quantities being wrapped into the material property in the frequency domain [14]. Following these suggestions, the reciprocal Willis equations (equations (2.1) with *Ψ*_{ijk}=−*S*_{ijk}/*iω*, one obtains a new set of constitutive equations
*Ψ*_{ijk} is a purely real material property in the low frequency limit that represents Willis coupling. This time domain formulation is particularly interesting because it clearly indicates that the generalized relationships of elastodynamics for stress and momentum degenerate to the common representations if strain and velocity are constant. In other words, the time domain formulation of equations (4.11) implies that Willis coupling is the result of microstructure that does not influence the static stress–strain response of the heterogeneous medium or rigid body dynamics, but that becomes observable under time-varying excitation. This is particularly interesting because it indicates that Willis coupling does not violate Galilean invariance of the stress field in a medium with uniform and constant velocity. As mentioned above, a similar set of constitutive equations were proposed by Nassar *et al*. after analysing a one-dimensional periodic system of masses and springs [14]. Their proposal was based on the requirement that the imposed stress field be invariant under Galilean transformations and then writing the definition of the momentum in similar manner. Lastly, it is important to note that despite the fact that the time domain formulation was deduced from causality using a highly idealized dispersionless limiting case, *Ψ*_{ijk} is not restricted to be dispersionless in general. Indeed, previous work on the topic clearly indicates that this material property may be strongly dispersive when resonant subwavelength material inhomogeneities are present in a medium [9,17].

As equations (4.11) are simply an alternative representation of the Willis equations, the restrictions on *S*_{ijk} due to passivity may be easily converted to restrictions on *Ψ*_{ijk} (the restrictions imposed by reciprocity are implicitly already included). In particular, for a reciprocal and passive system we require
*Ψ*_{ijk} is the imaginary part, analogous to the lossy parts of the stiffness and density tensors. Furthermore, using an analysis similar to that given above, the Kramer–Krönig relations for *Ψ*_{ijk} may be written as

## 5. Special cases

In this section, additional requirements of material properties for three special cases of particular interest are considered. First, the case of media supporting wave propagation with very low loss is considered. Then, the case of media without any resonances is discussed. Finally, the requirements on the fully general elastodynamic material properties are reduced to the acoustic limit, where neither shear stress nor strain are supported.

### (a) Low loss regions

Propagation with low losses in wave energy is an important subject, and much has been written on the topic [24]. If the net power flux through a region is set equal to zero, then

While the presence of a band structure precludes the use of the lossless restrictions, Landau & Lifshitz [24] demonstrated that interesting requirements may be found in the lossless regions. If *A*′′(*λ*)=0 for *λ* in some passband, then the principal integral on the right-hand side of equation (4.7a) may be reduced to a normal integral for *ω* in the passband. Then, taking a derivative with respect to *ω*, one may write
*A*′′(*λ*). But for the stiffness and density, passivity requires the imaginary part of the material properties to keep the same sign for all frequencies. Therefore, the sign of ∂*A*′/∂*ω* is the same as the sign of *A*′′(*ω*) in the passbands. It should be noted that because passivity only places a restriction on the *magnitude* of the coupling tensors, a similar statement for the coupling tensors cannot be made.

### (b) Relaxing media

In the regions of loss the Kramers–Krönig relations may still yield information about the slope of *A*′(*ω*). Assuming non-negligible *A*′′(*ω*), O’Donnell *et al*. [27] were able to expand the principal value integral in equation (4.7b) about the pole in an infinite power series equivalent to
*ω* or for *A*′′(*ω*) approximately proportional to *ω* (a property of relaxing media), the higher-order terms become negligible, and we find the sign of ∂*A*′/∂*ω* has the opposite sign of *A*′′(*ω*). This holds for all of the material properties, including the coupling terms.

The cases of Landau & Lifshitz [24] and of O’Donnell *et al*. [27] lead to opposite conclusions, so the situations in which they are applicable should be clearly delineated. First, the two approaches start from different underlying assumptions, and thus are described using different relationships. As Landau & Lifshitz assume very small *A*′′(*ω*), they start with the Kramers–Krönig relation for *A*′(*ω*). On the other hand, O’Donnell *et al.* make no assumption on *A*′′(*ω*), and use the Kramers–Krönig relation for *A*′′(*ω*). Furthermore, O’Donnell *et al.* tacitly assume the dominant contribution to the principal value integral comes from the pole, which implies there are no resonances of *A*′(*λ*) for *λ* close to *ω*. An example where this assumption breaks down is a low-loss material displaying band structure due to sub-wavelength resonances. In a passband far from any band gap, the dominant contribution to the integral in equation (4.7b) is the pole, and the limit of O’Donnell *et al.* holds. However, close to the edge of the passband, the relatively large values of *A*′′(*ω*) in the band gap dominate, and the limit of Landau & Lifshitz becomes appropriate. Figure 1 provides an illustration of these two approximations and their approximate domains of applicability.

### (c) Acoustic limit

Because of the recently increased interest in acoustic metamaterials, the causal form of the Willis equations, equations (4.11), is specialized here to the fluid case where shear is not supported. Instead of using the Cauchy stress and strain tensors, the acoustic pressure *P* and the volumetric strain *ε*_{V}, which are both scalars, will be used.

The volumetric strain is defined as the trace of the strain tensor, and the pressure is defined as negative one-third of the trace of the Cauchy stress
*σ*_{ij}=−*Pδ*_{ij} and *ε*_{ij}=*ε*_{V}*δ*_{ij}/3, respectively. Then equations (4.11) become
*κ*≡*C*_{iijj}/9 and *ψ*_{i}≡*Ψ*_{jji}/3 so that we may write
*κ*′′≤0, the requirements on *ρ*_{ij} do not change, and the passivity restrictions on the coupling parameter are given by

## 6. Conclusion

Restrictions on the range of physically meaningful Willis material properties have been derived based on the principles of reciprocity, passivity and causality. These restrictions are summarized in table 1. Reciprocity leads to minor symmetries of the stiffness and coupling tensors. Passivity leads to restrictions on the imaginary parts of the stiffness and mass density tensors and the real parts of the coupling tensors. Causality results in the Kramers–Krönig equations and relate the dispersion relations of the real and imaginary parts of the material properties. The standard representation of the linear Willis constitutive equations is shown to be acausal, and potential alternative representations have been presented.

## Authors' contributions

C.F.S. derived the reciprocity and passivity results for acoustic media, and M.B.M. derived the same and the restrictions due to causality for elastic media. A.A. and M.R.H. provided valuable insight which significantly guided the derivations given.

## Conflict of interests

We declare we have no competing interests.

## Funding

This work was supported by ONR through MURI Grant No. N00014-13-1-0631.

## Acknowledgements

The work for this paper was performed independently by the authors.

## Appendix A. One dimensional asymmetric fluid

Willis coupling (even) can arise from inherent asymmetry of the microstructure of a material [10,15]. To demonstrate this, consider an element represented by two unequal masses and a spring between them (figure 2). Let the first mass be *m*_{1}=*ρ*_{1}*S*Δ*x*/2 at *x* and the second mass be *m*_{2}=*ρ*_{2}*S*Δ*x*/2 at *x*+Δ*x*, and let the spring have a spring constant *k*. In these expressions, *S* is the cross-sectional area and Δ*x* is the width of the fluid element.

There are two equations of motion for each of the masses within the fluid element are given by
*ξ*(*x*) is the displacement at the field point, *P*(*x*) is the pressure field, and the over-dots denote a time derivative. These equations may be rearranged to yield the equivalent equations
*y*=*x*+Δ*x*/2 and neglect terms *O*[(Δ*x*)^{2}]
*ξ*′(*y*) as the volume strain *ε*(*y*) and *v*(*y*), and noting that −*P*′(*y*) is the time derivative of the momentum density *y* is not shown)

From this form, we readily identify the approximate effective density *ρ*=(*ρ*_{1}+*ρ*_{2})/2, the effective bulk modulus *κ*=−*k*Δ*x*, and the (reciprocal, passive and causal) coupling parameter *ψ*=Δ*x*(*ρ*_{2}−*ρ*_{1})/4. Note that *k*=*O*(1/Δ*x*), so the coupling terms are higher-order expansions of the effective constitutive equations

- Received July 29, 2016.
- Accepted September 26, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.