## Abstract

Homogenization of the equations of motion for a three-dimensional periodic elastic system is considered. Expressions are obtained for the fully dynamic effective material parameters governing the spatially averaged fields by using the plane wave expansion method. The effective equations are of Willis form with coupling between momentum and stress and tensorial inertia. The formulation demonstrates that the Willis equations of elastodynamics are closed under homogenization. The effective material parameters are obtained for arbitrary frequency and wavenumber combinations, including but not restricted to Bloch wave branches for wave propagation in the periodic medium. Numerical examples for a one-dimensional system illustrate the frequency dependence of the parameters on Bloch wave branches and provide a comparison with an alternative dynamic effective medium theory, which also reduces to Willis form but with different effective moduli.

## 1. Introduction

The effective medium concept has proven to be very useful for modelling composite materials with periodic microstructure. Efficient methods exist for computing the effective elastic moduli under static and quasi-static conditions. A more challenging task is to define frequency-dependent dynamic effective constants which are capable of describing phononic band gaps of the wave spectrum (Liu *et al.* 2005) at the expense of non-classical features, such as negative (Avila *et al.* 2008) and tensorial (Milton & Willis 2007) density. Effective material models that can retain a dispersive signature of the initial material microarchitecture have significant potential application in the design of sonic metamaterials. The purpose of this paper is to provide an analytical formulation for the effective material parameters of a periodic elastic medium at finite frequency.

A general framework for describing dynamic effective medium theories has been developed by Willis, beginning with variational arguments for non-local behaviour in mechanical (Willis 1981*a*, 1984) and electromagnetic systems (Willis 1981*b*). The predicted constitutive relations modify both the elasticity (Willis 1983) and the inertia (Willis 1985), and include coupling between stress/strain, on the one hand, and momentum/velocity, on the other (see Willis (1997) for a review). The notion of anisotropic inertia has appeared in several mechanical contexts, such as stratified layers of fluids (Schoenberg & Sen 1983) and elastic composites with strongly contrasting constituents (Avila *et al.* 2005). Coupling between particle velocity and stress, although not usually considered in constitutive theories for continua, is to be expected for theories relating spatial means of these quantities because of the underlying inhomogeneity of the displacement and stress fields within the region being averaged. Interest in the Willis constitutive model has intensified with the observation (Milton & Willis 2007) that the non-classical material properties could be realized, in theory at least, by suitably designed microstructures with internal springs, masses and gyroscopes. Derivations of the Willis equations have been demonstrated for periodic systems in one (Willis 2009; Nemat-Nasser & Srivastava 2011; Shuvalov *et al.* 2011) and three dimensions (Srivastava & Nemat-Nasser 2012), and the structure of the equations has been rigorously proved for both electromagnetic and elastic materials (Willis 2011).

Despite the growing awareness that the Willis constitutive model is the proper dynamic effective medium formulation for periodic systems there is as yet no simple means of calculating the effective material parameters at finite frequency. Several different techniques have been proposed to address this deficiency, none of which is easily implemented for arbitrary three-dimensional problems. Thus, closed form expressions were obtained in Willis (2009) for one-dimensional stratification using Floquet modes, where the latter have to be determined numerically. A more general numerical technique has recently been developed and applied to one- (Nemat-Nasser & Srivastava 2011) and three-dimensional (Srivastava & Nemat-Nasser 2012) composites. This approach employs a background or reference material, and solves for the polarization field (difference in stress/momentum between the actual and that of a comparison body, Willis 1980) to arrive at expressions for the effective parameters in terms of system matrices. The numerical procedure, which uses a combination of plane wave expansions (PWEs) and discretization schemes, is complicated by the required selection of a comparison material, the choice of which has no bearing on the final outcome. A quite different technique (Shuvalov *et al.* 2011) suited to one-dimensional periodically stratified elastic solids uses the monodromy matrix (i.e. propagator for a single period) to explicitly define effective material parameters from its logarithm. All of these methods derive the density, stiffness and the coupling parameter, *S*, as functions of frequency and wavenumber. This means, perhaps surprisingly, that equations of Willis form but with different parameters can be obtained for the same medium, depending on the homogenization scheme employed, an aspect discussed in §6.

Our purpose here is to develop, in so far as possible, analytical expressions for the dynamic effective material parameters of a periodic elastic material. The homogenization scheme seeks the constitutive relations and the equations of dynamic equilibrium governing the spatial average of the *periodic* part of the Bloch wave solution. The initial heterogeneous system is here assumed to be defined by pointwise constitutive equations of the form proposed by Willis, which includes ‘classical’ elastodynamics as the case of isotropic density and zero coupling between stress and momentum. We show that the homogenized equations are also of Willis form with dynamic effective material parameters that are non-local in space and time, thus demonstrating that the Willis constitutive theory is self-consistent and closed. The homogenized effective parameters are conveniently expressed in terms of PWE of the original material parameters, and accordingly we call this the PWE method. This approach allows us to define effective constants for any frequency-wavenumber combination, including, but not restricted to, values of {*ω*,**k**} on the Bloch wave branches.

We begin in §2 with an overview of the Willis constitutive equations followed by a summary of the PWE homogenization results. Section 3 is devoted to the derivation of the PWE effective material parameters for a scalar wave model that is simpler than the fully elastic problem, but which exhibits the main features of the homogenization scheme. The elastodynamic effective parameters are obtained in §4 and some of their properties are examined in §5. Examples of the PWE effective moduli are provided in §6 based on numerical calculations for a one-dimensional system, and comparisons are made with the parameters predicted by the monodromy-matrix homogenization scheme (Shuvalov *et al.* 2011). Conclusions and final words are given in §7.

## 2. Problem set-up and summary of the solution

### (a) Willis elastodynamic equations for a heterogeneous medium

The field variables are the vectors of displacement, momentum and body force, **u**, **p** and **f**; the symmetric second-order tensors of stress, strain and strain source, ** σ**,

**,**

*ε***, with . They are related by the system of elastodynamic equations which we take in the Willis form (in order to demonstrate in the end that this model is closed under homogenization for periodic materials). Accordingly, the material parameters are the second-order tensor of mass density**

*γ***, fourth-order elastic stiffness tensor**

*ρ***and third-order coupling tensor**

*C***. The components of the above material tensors in an orthogonal basis of the coordinate space satisfy the symmetries 2.1The two former identities imply that the density and stiffness tensors are Hermitian:**

*S***=**

*ρ*

*ρ*^{+}and

**C**=

**C**

^{+}(* and

^{+}hereafter denote complex and Hermitian conjugation). Assuming both

**and**

*ρ***C**positive definite ensures positive energy (2.3), see below. Let

**S**be generally complex. Field variables are functions of position

**x**and time

*t*, while the material tensors in the initial heterogeneous medium are considered as functions of

**x**only (the homogenization leads to the PWE effective parameters which, after inverse Fourier transform, are functions of both

**x**and

*t*).

The equation of equilibrium and the generalized constitutive equations in Willis form (Milton & Willis 2007; Willis 2011) are, respectively,
2.2where the dot implies time derivative. Note that the terms of (2.2) with **S** taken in component form imply and , where the permutable indices are enclosed in parentheses (see §5*a* for further details). By (2.1) and (2.2) with ** γ**=0, the time-averaged energy density is
2.3It is obvious that equations (2.2) include ‘classical’ linearly anisotropic elasticity with

*ρ*

_{ij}=

*ρδ*

_{ij},

*C*

_{ijkl}=

*C*

_{klij}=

*C*

_{jikl}and

*S*

_{ijl}=0.

Equations (2.2) may be solved to find **u**, **p** and ** σ** for given forcing functions

**f**and

**. Accounting for the effect of two driving-force terms is an important ingredient of homogenization concepts (Willis 2011). In particular, dependence of the solutions on the forcing terms will turn out to be crucial in finding the PWE homogenized equations.**

*γ*### (b) The periodic medium

We consider the material with **T**-periodic parameters:
2.4where (*d*=1,2 or 3), , and the linear independent translation vectors define the irreducible unit cell (*t*_{j}∈[0,1]) of the periodic lattice. Let {**e**_{j}} with *j*=1,…,*d* be an orthonormal base in . Denote
2.5where **a**_{j}·**b**_{k}=*δ*_{jk} (centre dot (·) is the scalar product in ), superscript (^{T}) means transpose and is an element of the set *Γ* of reciprocal lattice vectors whose components {*g*_{j}} in the base {**e**_{j}} take all values . The Fourier transform maps a **T**-periodic function *h*(**x**) to a vector in the infinite-dimensional space *V* _{g} associated with **g**∈*Γ*, where the components of are Fourier coefficients of *h*(**x**):
2.6(all quantities in the Fourier domain are indicated by a hat hereafter). In particular, the average over the single cell is defined by
2.7Practical calculations in the Fourier domain imply truncation of vectors and matrices in *V* _{g} to a finite size, i.e. reducing *V* _{g} to finite dimension. This is tacitly understood in the following when using the concepts of matrix determinant and inverse as convenient shortcuts.

The governing equations (2.2) with periodic coefficients (2.4) will be solved for both the wave field variables and the force terms in Bloch form
2.8where *h*(**x**) implies *the periodic part* of the full dependence of *h* on **x** (this convention reduces subsequent notation), and the frequency *ω* and wavevector **k** are assumed real-valued. Note that one may think of **f** and ** γ** in the form (2.8) as phased forcing terms which ‘drive’ the solution.

### (c) Summary of homogenization results

For *h* expressed in Bloch wave form, (2.8), define the effective field variable
2.9Using the effective sources (**f**^{eff},*γ*^{eff})=(〈**f**〉,〈** γ**〉) e

^{i(k·x−ωt)}implies ignoring the influence from the forcing

**f**(

**x**),

**(**

*γ***x**) with zero average over a single cell, which is natural for the homogenization theory aimed at recovering the effective wave fields of the form (2.9). One of the principal results of the paper is that, given the periodic medium with the Willis form (2.2) of the governing equations (which incorporates ‘classical’ elastodynamics model), the equations describing relations between the effective wave fields

**u**

^{eff},

*σ*^{eff},

**p**

^{eff}and and the effective forcing

**f**

^{eff}and

*γ*^{eff}are also of Willis form: 2.10The effective parameters

**C**

^{eff},

*ρ*^{eff}and

**S**

^{eff}are non-local in both space and time (despite the fact that the exact parameters of the original periodic material are stationary). Correspondingly, they are functions of frequency

*ω*and wavenumber vector

**k**in the transform domain , where their explicit expressions are as follows: 2.11a 2.11b and 2.11cwith the summation on repeated indices hereafter implicitly understood. (Noting that the indices of tensors

**C**and

**S**are appropriately split between and , we assume for brevity that

*d*=3 and so all roman indices run from 1 to 3). The tensors , and are defined on where

*V*

_{g≠0}is the Fourier vector space restricted to

**g**∈

*Γ*\

**0.**Each ‘element’ , and with fixed indices

*i*…

*q*is, respectively, a vector and a matrix in

*V*

_{g≠0}. Their components are as follows: 2.11d 2.11e 2.11f and 2.11gwhere . Note that since the matrix is Hermitian (), the effective material tensors (2.11) established in this paper recover the symmetries (2.1) of the corresponding tensors of the initial periodic material (although not the sign definiteness, see §5

*a*for more details). The Willis formulation of elastodynamics equations is therefore closed under the PWE homogenization.

## 3. Effective parameters: derivation for a scalar wave model

Instead of proceeding with the elastodynamic equations (2.2) in full, we first provide a detailed derivation for the case of a scalar wave equation. The scalar case is more revealing as it avoids the multitudinous suffices of the equations of elasticity while retaining the essence of the homogenization. We therefore ignore the obvious fact that an uncoupled acoustic scalar wave in solids with **x**-dependent material parameters is restricted to with *d*=1 or 2, since the scalar problem is purely a methodological shortcut here. Once the scalar derivation has been completed, it is straightforward to generalize the results to the original elastodynamics system (2.2) in *d*=1, 2 or 3, see §4.

### (a) Exact solution of the scalar wave problem

For the scalar wave model, the displacement *u*, momentum *p*, body force *f* and density *ρ* are scalars, while the stress *σ*, strain ** ε**=

**∇**

*u*, strain source

**and coupling parameter**

*γ***S**are vectors and the ‘shear’ stiffness

*μ*is a

*d*×

*d*matrix in (disregarding the trivial case of

*d*=1). According to (2.1),

*ρ*is real and

*μ*Hermitian while

**may be complex. Let these material parameters be**

*S***T**-periodic in the sense (2.4) and take the wave variables and sources in Bloch form (2.8). Then the scalar-wave version of the governing equations (2.2) in component form, with

_{,j}≡∂/∂

*x*

_{j}, is 3.1Applying Fourier transform (2.6) reduces (3.1) to the following algebraic equations in

*V*

_{g}: 3.2It is convenient to pass to a matrix form in

*V*

_{g}while keeping explicit the coordinate indices in

**x**-space (this makes the eventual generalization to the fully elastodynamic problem clearer). Recall the notation for the vector in

*V*

_{g}with the components (see (2.6)), and introduce the matrices , , , (∈

*V*

_{g}⊗

*V*

_{g}when

*j*,

*l*are fixed) with the components 3.3Note that , etc., i.e. the matrices and are Hermitian, while is not once

*S*

_{l}(

**x**) is complex. Equations (3.2) can now be expressed in a matrix form on

*V*

_{g}as 3.4where and .

Reduction to a single equation for the displacement provides a definition of the impedance matrix in the Fourier domain:
3.5We note that the condition for existence of free waves in the absence of sources, defines the Bloch eigen-frequencies *ω*(**k**). Our interest is ultimately in homogenized equations of motion allowing for the Bloch waves. It is however simpler to first proceed under the assumption that is not singular (). The issue of how to handle a singular is resolved afterwards in §3*c*, by introduction of a slightly modified impedance.

Assuming is non-singular, the Green matrix function is defined as the inverse of . The exact scalar solution then follows from (3.5) as 3.6

### (b) Homogenized equations

Taking the average (2.7) of the wave solution (3.6) and the constitutive equations (3.4)_{2} yields
3.7Our objective is to rearrange equations (3.7) into a form relating the effective wave fields to the effective forcing terms with the latter viewed as independent variables. As a first step in this direction, we insert the effective sources (see equation (2.9)) and in equations (3.7) and rewrite the result
3.8where the scalars *a*, , the vectors ** α**, and the matrix are
3.9Note that

*a*is real and

**A**

^{+}=

**A**. Eliminating 〈

*f*〉 using (3.8)

_{1}allows us to recast equations (3.8)

_{2}and (3.8)

_{3}in a form reminiscent of the Willis constitutive relations, 3.10where

*X*

_{j}and

*Y*are at this stage arbitrary. Comparing (3.10) with the assumed structure of the effective constitutive relations applied to fields in the presence of a strain source, equation (2.10)

_{2}and bearing in mind (2.9) yields 3.11Equating the contributions to 〈

*σ*

_{j}〉 and 〈

*p*〉 in (3.10) and (3.11) from

**gives, in turn, 3.12Full equivalence requires**

*γ**Y*=

*ρ*

^{eff}and

**X**=

**S**

^{eff}, which implies two more identities found by equating the contributions from 〈

*u*〉 to 〈

*p*〉 and to 〈

*σ*

_{j}〉, respectively: 3.13The first equation provides an expression for

*ρ*

^{eff}and the second is an alternative to (3.12)

_{2}for

**S**

^{eff}(equivalence of (3.12)

_{2}and (3.13)

_{2}is noted below).

Combining equations (3.12) and (3.13)_{1} and reinstating the notations from (3.9) yields the effective parameters in the form
3.14a
3.14b
and
3.14cIn the derivation of (3.14), we have taken note that and used identities such as
3.15The equivalence of the expressions (3.12)_{2} and (3.13)_{2} for **S**^{eff} follows from using the explicit formulas (3.14) and identities of the above type.

### (c) Regularized form of the effective parameters

By definition (3.6), the Green function diverges on the Bloch branches *Ω*_{B} (3-parameter manifold in {*ω*,**k**}-space) defined by the Bloch dispersion equation
3.16where the null vector of is the Bloch eigenmode in the Fourier domain. At the same time, the effective dynamic parameters given by equations (3.14) generally remain finite on *Ω*_{B} (accidental exceptions such as degenerate and certain other points on *Ω*_{B} are disregarded for the moment and will be discussed in §3*d*). This can be shown by inspecting the limit of equations (3.14) as *ω*,**k** tend to a (single) point (*ω*,**k**)_{B} on a Bloch branch, from which it is seen that the right-hand side members of (3.14) that diverge on *Ω*_{B} in fact compensate and cancel each other. However, such limiting behaviour of equations (3.14) prevents their numerical implementation by evaluating each member independently before summing them up. The difficulty can be circumvented if we recast equations (3.14) in the analytically equivalent but explicitly different form that is rid of the terms diverging on *Ω*_{B}. Since the divergence in (3.14) is due to Green's function the idea is to extract the part of which is regular on *Ω*_{B}. Consider such a partitioning of in the general form where and as It is not unique. For instance, one standard way is to invoke the spectral decomposition
3.17where *ς* and are the corresponding eigenvalue and eigenvector of which satisfy and as and (the pseudoinverse of ) is regular on *Ω*_{B}. Substituting (3.17)_{2} for in (3.14) removes the diverging terms and yields an explicit form expressed through . Unfortunately, explicit calculation of is relatively laborious. In this regard, we advocate the use of another partitioning of into regular and diverging parts which requires merely singling out Fourier components with zero and non-zero **g**. This is an essential ingredient of our homogenization method that ensures both robust and straightforward numerical implementation.

With the above idea in mind, denote by *V* _{g≠0} the vector space that is associated with the set *Γ*\**0**, i.e. is orthogonal to the vector (see (2.5) and (2.7)). Let be the restriction of over this space *V* _{g≠0} and be the inverse of on *V* _{g≠0}:
3.18(note that as defined is *not* a restriction of on *V* _{g≠0}). Introduce the vectors , , ,
3.19Also denote and which is a slight generalization of the average (2.7) to include the projection of a matrix onto **g**,**g**′=**0**. Thus
3.20Using the notations (3.19) and (3.20), we note the identical decompositions
3.21where is adjugate of . Hence if as assumed in (3.18)_{2}, then
3.22a
3.22b
3.22c
and
3.22dwhere equation (3.22d) yields the sought-for partitioning into two parts, one regular and the other diverging on Bloch branches *Ω*_{B}. Note that generally (i.e. for ) the quantities and in (3.22d) are different from, respectively, the pseudoinverse and the eigenvector in (3.17).

Writing (3.22d) in the transparent form
3.23and using the notations (3.19) reduces equations (3.14) to the simpler form
3.24a
3.24b
and
3.24cIt is emphasized that the effective parameters explicitly defined in (3.24) are by construction the same as those in (3.14). At the same time, the terms in (3.14) divergent on *Ω*_{B} are eliminated in (3.24), in this sense (3.24) are uniformly convergent. The new expressions (3.24) are determined from the regular part of the Green function associated with **g**∈*Γ*\**0** which is straightforward to compute; moreover, their explicit form is simpler than that of (3.14).

### (d) Exceptional cases

The expressions for the effective parameters (3.24) defined by the matrix are obtained for the general case of non-singular . i.e. . If happens to be singular, then numerical implementation of (3.24) produces diverging terms (though equations (3.14) expressed via remain finite for so long as ). The properties of the wave solutions that occur specifically at values of *ω*,**k** rendering singular are of interest. Denote the set of these exceptional points (3-parameter manifold in {*ω*,**k**}-space) by *Ω*_{exc} so that
3.25Note that equations (3.22) are invalid on *Ω*_{exc}.

Consider first ‘generic’ points of *Ω*_{exc} which lie off the Bloch branches (3.16), i.e. (*ω*,**k**)_{exc} *Ω*_{B} (). Then by (3.20)_{2} and hence the averaged wave (3.8)_{1} is It is generally non-zero, which is not surprising as 〈*u*〉=0 would correspond to the zero effective response to the non-zero effective forcing. At the same time, the wave response to a pure force source (*γ*_{j}=0) at the points (*ω*,**k**)_{exc} *Ω*_{B} has a zero mean over the unit cell **T** in **x**-space. It is therefore somewhat satisfying and physically sensible that the effective parameters (3.24) indicate such exceptional points by becoming numerically divergent.

Now assume a generally non-empty set *Ω*_{exc}∩*Ω*_{B} (2-parameter manifold in {*ω*,**k**}-space) of the points of Bloch branches (3.16) at which happens to be singular. By (3.21)_{2}, simultaneous equalities and yield It is natural to assume that is non-zero. Then it is a null vector of and it is orthogonal to **w**. Hence, by (3.21)_{1}, is also a null vector of , i.e. it is a Bloch mode . Thus, if incidentally occurs at some point of a Bloch branch, then the Bloch mode at this point is , which is by construction orthogonal to **e** and so its original *u*_{B}(**x**) in **x**-space has zero mean 〈*u*_{B}〉=0 over the unit cell **T**. Let us call it a *zero-mean Bloch mode*. This is in contrast to ‘normal’ points of *Ω*_{B} where the Bloch modes are defined by (3.22c) as and thus have non-zero mean 〈*u*_{B}〉=1. Since zero-mean Bloch modes are zero (trivial) solutions of the effective equations for free waves, their explicit exclusion from the framework of the homogenization theory, as signalled by divergence of the effective material parameters (3.24), is physically consistent. Note that the divergence at the points of zero-mean Bloch modes is ‘genuine’ (analytical), i.e. it occurs regardless of which explicit expressions are used for the effective parameters (unless more restrictions are imposed apart from *Ω*_{exc}∩*Ω*_{B}, such as ).

In conclusion, we note another exceptional case for the effective parameters (3.24) is the possible existence of double points (self-intersections) on Bloch branches *Ω*_{B}. It can be shown, using an example the spectral decomposition and the pseudoinverse of , that a double point on *Ω*_{B} causes a genuine (in the sense mentioned above) divergence of the effective parameters, which conforms to the fact that a double point admits the construction of a zero-mean Bloch mode 〈*u*_{B}〉=0 from a pair of null vectors of .

### (e) Summary for a scalar wave problem

The equations governing the averaged effective fields *u*^{eff}, *σ*^{eff}, *p*^{eff} and *ε*^{eff}=**∇***u*^{eff} and their relation to the averaged applied body force *f*^{eff} and strain source *γ*^{eff} (see (2.9)) are of Willis form:
3.26where *ρ*^{eff}, and are non-local in space and time. They are defined in the transform domain by equations (3.24) with the following explicit form
3.27awhere is the inverse to the matrix with the components
3.27bNote that *μ*^{eff}=*μ*^{eff+}, *ρ*^{eff} is real and is complex-valued, see §5*a*.

The effective parameters defined by equations (3.27) are regular at any *ω*,**k** excluding possible double points (intersections) of Bloch branches, and the exceptional points (*ω*,**k**)_{exc}∈*Ω*_{exc} where (a numerical divergence appears in the vicinity of such points). If a point (*ω*,**k**)_{exc} incidentally occurs on a Bloch branch then the corresponding Bloch eigenmode *u*_{B}(**x**) has a zero mean 〈*u*_{B}〉=0 over the unit cell **T**. Zero-mean Bloch modes can also exist at the double point of a Bloch branch. Barring these exceptions, equations (3.27) are regular on Bloch branches which are defined by the dispersion equation (see (3.22a))
3.28The term which defines the dispersion relation is expressed as a more physically meaningful quantity in §5*b*.

## 4. Elastodynamic effective parameters

The procedure of adapting the scalar acoustic model to the vector elastodynamic model requires that the scalar and vector variables *u*, *p*, *f* and ** σ**,

**,**

*ε***of the acoustic case become, respectively, vectors and symmetric second-order tensors in . Accordingly, the scalar density**

*γ**ρ*, coupling vector

**S**and the stiffness matrix

**increase by two orders as tensor and become**

*μ***,**

*ρ***S**and

**C**defined in §2

*a*. The component form of the various quantities for the two problems are related as follows: 4.1Based on this equivalence the generalization of the scalar equations (3.24) proceeds by first replacing the single suffix

*j*in the former with

*ij*in the latter, then using definitions such as those in equation (3.19) to assign additional suffices, ultimately yielding the elastodynamic result (2.11).

Expressing the effective elastodynamic parameters (2.11) via the matrix ensures that they are generally regular on the Bloch branches. This is an important feature that deserves further scrutiny. By analogy with (3.16), denote the manifold of Bloch branches for vector waves by *Ω*_{B} (bold-lettered to distinguish it from *Ω*_{B} of the scalar case), so that
4.2where the matrix and the vector have components and in If and so is well-defined (see (2.11)), then
4.3This is similar to the scalar version (3.22a) and (3.22b) with the important distinction that now and are 3×3 matrices in and , . Their components in are defined by
4.4(the quantity **Z**^{eff} is expressed in more physical terms in §5*b*). Note also that the vector-case generalization of (3.22c) implies that the null vector of (=**Z**^{eff}) on *Ω*_{B} is equal to the average 〈**u**_{B}〉 of the Bloch eigenmode defined by (4.2)_{2}. Finally, consider the manifold *Ω*_{exc} of exceptional points where (cf. (3.25)) and the explicit expressions (2.11) diverge. For the points in *Ω*_{exc} *Ω*_{B}, i.e. lying off the Bloch branches, we have (cf. (4.3)_{1}) and hence the effective force **f**^{eff}(**x,***t*) (2.9) with the amplitude taken to be a null vector of produces at ** γ**=0 a zero effective response 〈

**u**〉=

**0.**This is again similar to the scalar case (except that the vector 〈

**f**〉 has to be specialized while 〈

*u*〉=0 ∀〈

*f*〉 at ,

**=0). If the intersection**

*γ*

*Ω*_{exc}∩

*Ω*_{B}occurs, i.e. occurs on a Bloch branch, then, contrary to the scalar case, the vector Bloch mode does not generally have a zero mean 〈

**u**

_{B}〉=

**0**unless the null vector of also satisfies (see equation (3.21)

_{1}).

## 5. Discussion of the effective model

### (a) Properties of the effective constants

As it was noted in §2*c*, the effective material tensors given by equations (2.11) retain the symmetries (2.1) of the corresponding tensors of the initial periodic material. At the same time, the effective density and elasticity tensors *ρ*^{eff} (at *ω*≠0) and **C**^{eff} are not sign definite. This is not at all surprising since (2.3) no longer describes a positive energy for the homogenized dispersive medium. Dispersive effective tensors defined by (2.11) correspond to the non-local properties of the homogenized medium in the space–time domain where **C**^{eff},*ρ*^{eff} and **S**^{eff} depend on ** ξ**=

**x**−

**x**

^{′}and

*τ*=

*t*−

*t*

^{′}. Note that the effective tensors are certainly not periodic in

**k**(even if taken on the periodic Bloch branches

*ω*

_{B}(

**k**)=

*ω*

_{B}(

**k**+

**g**)), and hence they are not periodic in

**as well as non-stationary. It is therefore clear that the initial periodic medium with non-zero coupling**

*ζ*,**S**(

**x**) in (2.2) cannot be seen as resulting from PWE homogenization of another periodic material, say with a much finer scale. At the same time, by proceeding from (2.2) and arriving at (2.10) with the same entries of the coupling terms

**S**and −

**S**

^{+}, the PWE homogenization confirms the closure of the Willis elastodynamics model. In this regard, note that our notation

**S**

^{+}defined as is equivalent to

**S**

^{†}of Srivastava & Nemat-Nasser (2012) (see also footnote 1).

Consider in more detail the conventional situation where the initial periodic medium represents a classically elastic solid with zero coupling **S**(**x**)=0 and real-valued density *ρ*(**x**) and stiffness **C**(**x**). It is evident from (2.11) with that *ρ*^{eff} and **C**^{eff} are even functions and **S**^{eff} an odd function of *ω* and hence of *τ*. Dependence on **k** satisfies *ω*_{B}(**k**)=*ω*_{B}(−**k**) and
5.1By (5.1), **C**^{eff},*ρ*^{eff} are real and **S**^{eff} is pure imaginary at **k=0** and any *ω*. It also follows from (5.1) that **C**^{eff}(*ξ*,*τ*) and *ρ*^{eff}(*ξ*,*τ*) are real and **S**^{eff}(*ξ*,*τ*) is pure imaginary in the space–time domain.^{1}

Let *ρ*(**x**) and **C**(**x**) be even functions of **x**. This condition on its own implies that *ρ*^{eff} and **C**^{eff} are even functions and **S**^{eff} is an odd function of **k** and hence of ** ξ**. Combination with (5.1) implies that a classically elastic solid with a periodically even profile is characterized by

*ρ*^{eff}and

**C**

^{eff}that are real and even (in each variable) functions of either set of variables (

*ω*,

**k**) and (

**,**

*ξ**τ*), while

**S**

^{eff}is a real-valued odd function of (

*ω*,

**k**) which vanishes at

**k**=

**0**and it is a pure imaginary odd function of (

**,**

*ξ**τ*).

The case of an even profile is one example of possible symmetry. The general case may be outlined as follows. Let the orthogonal transformation **R** (**R**^{−1}=**R**^{T}) be an element of the symmetry group of the initial material, such that *C*_{ijkl}(**x**)=*R*_{im}*R*_{jn}*R*_{kp}*R*_{lq}*C*_{mnpq}(**x**) ∀ **x** and also *ρ*(**x**)=*ρ*(**Rx**) and **C**(**x**)=**C**(**Rx**). Then *ω*(**k**)=*ω*(**Rk**) and
5.2Consequently, the effective tensors *ρ*^{eff}, **C**^{eff} and **S**^{eff} are invariant to **R** at **k**=**0**, and are invariant at **k**≠**0** to those orthogonal transformations which also satisfy **Rk**=**k**.

### (b) Effective impedance

Assume the homogenized medium with the effective elastodynamics equations in the Willis form (2.10). The problem of finding the effective response to effective sources (2.9) averaged over the unit cell of the original periodic medium is now seen as a forced-motion problem:
5.3In particular, uniform amplitudes 〈**f**〉 and 〈** γ**〉 can now be interpreted as Fourier harmonics of arbitrary sources in the homogenized medium. Equations (2.10) with reference to (5.3) can be recast as
5.4where

**Z**

^{eff}is the effective impedance matrix with the components in 5.5Equation (5.4) considered in the absence of sources implies 5.6Equation (5.6)

_{2}is the dispersion equation for free waves in the homogenized medium, and as such may be called the Christoffel–Willis equation owing to its analogy to the Christoffel equation in a uniform medium. Note that, strictly speaking,

**Z**

^{eff}has been introduced above as the impedance of the homogenized medium. We now show that it can be equally interpreted as the effective impedance for the initial periodic material. The consistency of such dual interpretation of

**Z**

^{eff}is an essential verification of the model.

For clarity, consider first the case of a scalar wave problem, where the impedance (5.5) is a scalar
5.7and (5.6) reduces to *Z*^{eff}(*ω*,**k**)=0. Substituting the effective material parameters (3.27) in (5.7) yields the identity
5.8Now recall that by (3.22a), and that, barring certain exceptions (see §3*d*), is equivalent to the dispersion equation (3.16) defining Bloch branches in the actual periodic material. Thus, we arrive at equation (3.28), where is now an equality, and we observe that *Z*^{eff}(*ω*,**k**)=0 is the same dispersion equation providing both the spectrum of eigenmodes for the homogenized medium and the Bloch spectrum (in the main) for the actual periodic material. It is instructive to note from (5.8) that the effective impedance *Z*^{eff} generally differs from defined by the statically averaged material constants (see (3.20)). This observation highlights the fact that the PWE dynamic homogenization of a periodic material must necessarily proceed from a forced-wave problem.

The same conclusion follows for the elastodynamic case. Substituting from (5.3) into (5.5) yields the same expression as obtained for in (4.4) and therefore verifies that , now an equality. By (4.3), is, apart from exceptional cases, equivalent to the Bloch dispersion equation. Comparison with (5.6) corroborates the desired consistency. Note in (4.4), which implies the similar observation as noted above for the scalar case.

The symmetry of the effective material tensors (see §5*a*) may lead in a standard fashion to quasi-diagonal or diagonal structure of **Z**^{eff}. To be specific, consider an example where **k**≠**0** is parallel to one of the orthogonal translations **a**_{j} of a three-dimensional simple cubic lattice of spherical inclusions of locally isotropic material in a locally isotropic matrix material. Then, using elementary considerations of symmetry, the tensors *ρ*^{eff}, **C**^{eff} and **S**^{eff} must satisfy the tetragonal symmetry (4 *mm*) and hence **Z**^{eff} in the basis {**a**_{j}} is diagonal with two equal components. Their vanishing indicates the polarization of the averaged Bloch eigenmodes defined by (5.6)_{1} and also shows that the Bloch spectrum *ω*_{B}(**k**) defined by equation (5.6)_{2} contains some branches that are doubly degenerate for any non-zero **k**∥**a**_{j}.

A remark is in order concerning the effective properties for **k**=**0**. Normally, equations (5.5) and (5.6)_{2} imply that (*ρ*^{eff}=0 in the scalar-wave case) at the points **k**=**0**, *ω*(**0**)≠0 of the Bloch branches. It is however noted that **Z**^{eff}(*ω*,**0**)=−*ω*^{2}*ρ*^{eff} may not hold if *ω*(**0**) belongs to exceptional points of divergence of the effective tensors (see §§3*c* and 4). For example, assume a model case of periodic **C**(**x**) and constant *ρ*≡*ρ*_{0}. It is clear that **Z**^{eff}(*ω*,**0**)≠−*ω*^{2}*ρ*_{0}**I,** since the equality would mean that the Bloch spectrum has no open stopbands at **k**=**0** which is certainly not the case for an arbitrary periodic **C**(**x**). In fact, spatially averaging the pointwise equilibrium equation
5.9for **k**=**0** implies that any Bloch mode with *ω*≠0 must have zero mean, 〈**u**_{B}〉=**0**. In conclusion, if *ρ* is constant while **C**(**x**) is periodic, **k**=**0** definitely belongs to the exceptional subspace where **C**^{eff}(*ω*,**0**) and/or **S**^{eff}(*ω*,**0**) must diverge. A numerical illustration of this is given in §6.

### (c) Low frequency limit of the density

The form of the homogenized effective moduli of (2.11) while relatively compact belies their intricate dependence on the fixed material properties of the periodic system as well as on frequency and wavevector. Here, we touch on one small corner of the parameter space, the quasi-static limit, with particular attention to the effective inertia. Assume a classical initial medium (with **S**=0 and positive *ρ*(**x**)). Taking it is clear from (2.11) that to leading order in *ω* the effective moduli reduce to their static form, the effective density is the average, and **S**^{eff}=*O*(*ω*). The first correction to the effective density is positive-definite tensorial,
5.10This illustrates explicitly the departure of the homogenized density from that of classical elastodynamics, as expected on the basis of prior results for Willis equations parameters for periodic systems (Shuvalov *et al.* 2011; Srivastava & Nemat-Nasser 2012).

## 6. Numerical illustrations and discussion of the homogenization scheme

Numerical examples are presented for the one-dimensional scalar wave problem. All results reported are for the periodic medium based on the three-layer single cell defined by table 1. The PWE effective parameters *μ*^{eff}, *ρ*^{eff} and *S*^{eff} are taken on the Bloch branches *ω*_{B}(*k*) which reduces them to functions of a single variable, *k*. The first six branches *ω*_{B}(*k*) are shown in figure 1, and the corresponding values of effective parameters are displayed in figure 2.

Setting apart the obvious case of the origin of the first branch which starts with *S*^{eff}=0 and *ρ*^{eff}=〈*ρ*〉 at *k*=0 (see §5*c*), several other features are noteworthy in figure 2. In accordance with (5.1), Re *S*^{eff} vanishes at *k*=0 on all branches, while Im *S*^{eff} starts at *k*=0 with non-zero values that alternate in sign beginning with the positive value on the second branch. The effective density *ρ*^{eff} vanishes at *k*=0 on all branches above the first one, as expected from the form of the dispersion relation *Z*^{eff}(*ω*,*k*) = 0, which becomes, using (5.7),
6.1As *k* grows, *ρ*^{eff}(*k*) exhibits both positive and negative values, following the alternating trend shown by Im *S*^{eff}. The magnitude of *ρ*^{eff} remains within the range of the constituent material densities, with the exception of the fifth branch on which it is negative with for a range of *k* (here is the largest density in the periodic medium ). The calculated PWE *μ*^{eff} is always non-negative; however, its magnitude is exceedingly large near *k*=0 for higher order modes, e.g. the fifth and sixth. This can be partly understood from the dispersion relation (6.1) which permits *μ*^{eff} to become large in magnitude as . The enormous values on some branches, e.g. *μ*^{eff}>10 000 in figure 2*e*, may indicate proximity to an exceptional point related to a small value of (see §3*d*). The magnitude of *S*^{eff} is intermediate between that of the density and the stiffness. It is worth noting that (2.3) evaluated for the one-dimensional scalar Willis model displays negative values on the third and fifth branches. The proper form of the energy density for the dispersive effective medium is a delicate topic (Ruppin 2002), best left for a separate discussion.

Figure 3 illustrates the emergence of singular effective parameters at *k*=0, *ω*(0)≠0 when the density is uniform, as predicted in §5*b*. Although the curves shown are only for the second branch, the same dependence is observed in all higher branches.

An alternative route to homogenization of one-dimensional periodically stratified elastic body was proposed in Shuvalov *et al.* (2011). The basic idea is to compare the logarithm of the monodromy matrix (MM) **M**(*y*+*T*,*y*), which is the propagator over one period *T* along the stratification direction *Y* , with the matrix of coefficients of the system of elastodynamics equations for a uniform medium. In this way, one can identify the dynamic effective material parameters at any frequency wavevector on a Bloch branch as long as the logarithm is well defined. Unlike the PWE, the MM method does not require introducing force terms and is not based on averaging. It is worth emphasizing that the PWE and MM approaches are not two ‘technically’ dissimilar derivations of the same result but imply different models of a homogenized medium. In this light, the more significant is the result that the MM method also necessarily leads to the effective elasticity of Willis form (classical elasticity will not suffice) (Shuvalov *et al.* 2011). Note that the MM homogenization involves the choice of the reference point *y*=*y*_{0} at which the initial data for **M**(*y*_{0}+*T*,*y*_{0}) is prescribed. The MM homogenized media associated with different *y*_{0} have the same eigenfrequency spectrum *ω*(*k*) which by definition coincides with the Bloch branches of the actual periodic medium. At the same time, the values of MM effective constants depend on the choice of *y*_{0}. Once any *y*_{0} is fixed, the exact displacement and traction are recovered at periodically distant points *y*_{0}+*nT*. As a result, the MM homogenization can provide an exact solution of boundary-value or reflection problems for a finite structure containing several periods or for a periodic half-space in contact with other regions (Shuvalov *et al.* 2011).

As an example, figure 4 shows the effective parameters of the same material (table 1) which are calculated for one of the Bloch branches (the fifth) by the MM method. The calculation uses the formulas (4.2) and (4.7)_{1−3} of Shuvalov *et al.* (2011) with a standard definition of the 2×2 propagator **M**(*y*_{0}+*T*,*y*_{0}) of scalar waves through a trilayer. The two sets of data on figure 4*a*,*b* correspond to two different choices of the reference point *y*_{0} within a period. Both sets reveal the same general features of MM homogenization (Shuvalov *et al.* 2010, 2011) such as vanishing of *ρ*^{eff} and *μ*^{eff} at *k*=*π* and a pure imaginary value of *S*^{eff}. The latter stands in stark contrast to the generally complex-valued *S*^{eff} in the PWE homogenization model. At the same time, despite the quite disparate natures of the PWE and MM methods, it is noteworthy that they do reproduce the enormous values of *μ*^{eff} occurring at small *k* on the fifth Bloch branch. This may be seen as evidence that such behaviour of *μ*^{eff} is not a numerical artefact.

## 7. Conclusion

The principal result of the paper is a set of explicit formulas for calculating the homogenized material parameters using PWE expansions of the underlying elastic moduli and density, equation (2.11). The homogenized governing equations retain important physical properties of the original periodic system. Thus, they provide (i) the averaged response to averaged source(s), and (ii) the eigen-spectrum *ω*(**k**) for Bloch waves. The effective parameters of equation (2.11) are convergent for virtually all frequency and wavevector combinations, including {*ω*,**k**} pairs on the Bloch wave branches. The exceptional cases yield singular values of the effective parameters and correspond to zero-mean Bloch waves for the scalar wave system.

Some comments are in order on the PWE homogenization method and its outcome. A non-zero strain source is essential for the derivation of the effective parameters, although the strain source may be put to zero after the homogenization is complete. Inclusion of both body force and strain source is consistent with the findings of Willis (2011) that homogenization procedures based on Green's functions can lead to under-determined systems unless a full set of sources is included from the outset. The PWE homogenization procedure derives equations governing the average of the periodic part of the Bloch wave solution. The PWE scheme is not the only approach for defining an effective medium for periodic system. Another approach which applies for media with one-dimensional periodicity is the MM homogenization of Shuvalov *et al.* (2011). The effective governing equations are again of Willis form with explicit expressions for the effective parameters that are distinct from the PWE parameters, as the numerical example in §6 demonstrates. The fact that different sets of effective Willis equations can be obtained for the same medium points to a non-uniqueness in the definition of dynamic homogenization. The MM scheme was originally proposed as a method for layered structures that faithfully retains the exact phase relationship over multiple periods, although at only one point within the period (Shuvalov *et al.* 2011). The PWE scheme provides predictions for averaged quantities and does not yield pointwise exact solutions.

We have shown that the homogenization of a classically elastic medium (**S**=0) yields a Willis effective medium with either complex-valued **S**^{eff}(*ω*,**k**) (satisfying (5.1)) or with pure imaginary (§6). The results indicate the closure property for the Willis model with strictly complex **S** under homogenization.

The present formulation, in starting with a heterogeneous set of Willis elastodynamic equations, and resulting after homogenization in a set of homogeneous constitutive equations of the same form, confirms the completeness of the Willis model. This type of closure complements the fact that the Willis constitutive equations are closed under spatial transformation (Milton *et al.* 2006). Finally, we note that the present results only partially open the door on our understanding of the effective dynamic properties of phononic crystals. Much remains to be investigated regarding analytical aspects, numerical implementation and applications of these properties.

## Acknowledgements

A.N.N. acknowledges support from CNRS and from ONR.

## Footnotes

↵1 This observation makes consistent the use of −

**S**^{+}in (2.2) in place of**S**^{†}=**S**^{T}in Milton & Willis (2007, eqns (7.3) and (7.5)). Equation (5.1)_{3}implies that the term −**S**^{eff}^{+}in the homogenized equations (2.10) can be interpreted as**S**^{eff}^{+}(*ω*,−**k**) which is the formal adjoint of**S**considered as a differential operator, i.e.**S**^{†}, see Willis (2009, 2011). The variety of notation in the literature for the coupling**S**^{eff}does not make things simpler. In this regard, notation in Srivastava & Nemat-Nasser (2012) agrees with (2.10) and (5.1); notation in Willis (1997, eqn (3.26)) is similar to (2.10) although the conjugate of*M*_{ijk}(equivalent to*S*_{ijk}here) does not appear.

- Received November 28, 2011.
- Accepted January 13, 2012.

- This journal is © 2012 The Royal Society