## Abstract

A method to derive homogeneous effective constitutive equations for periodically layered elastic media is proposed. The crucial and novel idea underlying the procedure is that the coefficients of the dynamic effective medium can be associated with the matrix logarithm of the propagator over a unit period. The effective homogeneous equations are shown to have the structure of a Willis material, characterized by anisotropic inertia and coupling between momentum and strain, in addition to effective elastic constants. Expressions are presented for the Willis material parameters which are formally valid at any frequency and horizontal wavenumber as long as the matrix logarithm is well defined. The general theory is exemplified for scalar SH motion. Low frequency, long wavelength expansions of the effective material parameters are also developed using a Magnus series, and explicit estimates are derived for the rate of convergence.

## 1. Introduction

Elastic waves in periodically layered or continuous (functionally graded) elastic media of general anisotropy have been studied extensively using different methods. Among them is the sextic formalism of Stroh, which incorporates the elastodynamics equations into a first-order ordinary differential equation for the displacement-traction state vector with system matrix **Q** composed of material parameters (Ting 1996). The wave-field propagator matrix along the stratification direction *y*, **M**(*y*, 0), is given by the Peano series of multiple integrals of products of **Q**(*y*) (Pease 1965). This is essentially a power series in distance-to-wavelength ratio, which is therefore particularly well suited in tackling the problem of approximating a periodically stratified medium by an effective homogeneous medium. The Stroh formalism clarifies the meaning of zero-order homogenization, or static averaging, of a periodic medium by revealing the zero-order effective material parameters (Behrens 1968; Grimsditch & Nizzoli 1986) as nothing more than the matrix **Q**(*y*) integrated over the period *T*, which is the leading term of the logarithm of the Peano series for **M**(*T*, 0) (Norris 1992). Static averaging implies a non-dispersive effective medium. Generalization to a higher order effective homogeneous medium, which must be dispersive, is less obvious. Its derivation is commonly based on the long-wave dispersion of the fundamental Bloch or Floquet branches. Their onset in arbitrary anisotropic periodically stratified media was analysed by Norris (1992); based on this, the scalar-wave equation for a transversely isotropic dispersive effective medium was modelled by Norris (1992, 1993) and Norris & Santosa (1992) and in the subsequent literature (e.g. Andrianov *et al.* (2008) and its bibliography). A semi-analytical approach for general anisotropy (Wang & Rokhlin 2002; see also Potel *et al.* 1995) fits the long-wave Floquet dispersion to statically averaged effective constants, such that the effective medium is seen as a ‘continuum of non-dispersive media’ that are different for different frequency and propagation direction.

In this paper, a new method is proposed for finding the dynamic effective constitutive equations at finite frequencies and wavelength. This is achieved by explicit construction of effective spatially constant material coefficients that exactly reproduce the monodromy matrix, **M**(*T*, 0). The effective constitutive theory is exact in the sense that it gives the correct displacement-traction field at the unit-cell interfaces over arbitrarily long distance of propagation. Two key steps distinguish the method advocated. First is the idea, based on the Floquet theorem, of defining the effective medium such that the sextic system of elastodynamics equations in this medium has the matrix of coefficients, **Q**_{eff} equal to i**K**, where **K** is the Floquet wavenumber matrix with an exact definition . For the low-frequency long-wave range, a matrix logarithm admits an expansion called the Magnus series (Blanes *et al.* 2009). Restricting it to the leading-order term leads to the statically averaged effective model (see above). Taking the next-order term(s) of the Magnus expansion for **K** reveals that, unless the variation of material properties over a period is symmetric, the above-defined dispersive effective medium cannot be fitted to the standard form of elastodynamic equations, in which the frequency dispersion and non-locality would be fully accounted by the dependence of effective density and elasticity on *ω* and *k*_{x}. This motivates the second significant step in the present method, which is identifying a class of constitutive models that does fit i**K** with the system matrix, **Q**_{eff} of a homogeneous medium. We demonstrate by construction that the model described by the Willis constitutive relations with a dynamic stress-impulse coupling tensor (Milton & Willis 2007; Willis 2009) provides such a class of materials. Expansions of the Willis material coefficients based on the Magnus series for **K** are obtained in the low-frequency long-wave range where they are analytical in *ω*,*k*_{x}. Explicit estimates of the dependence on *ω*,*k*_{x} are found that enable closed-form asymptotics of the Willis coefficients with a desired accuracy. At the same time, the definition and the fitting of the matrix, **Q**_{eff}≡i**K** to the effective Willis material is formally not restricted to the low-frequency long-wave range.

A significant outcome of the proposed approach is that fairly explicit expressions are obtained for the effective material parameters. This is particularly evident for the example of SH (shear horizontal) waves discussed in detail in §4. In this regard, the approach is distinct from that of Willis (2009) which leads to expressions for the parameters in the (spatial) transform domain. Note also that the results here apply to a single realization of the layered medium, no ensemble averaging is invoked. This point is discussed further in §4.

The paper proceeds as follows. The problem is formulated in §2 where the sextic formalism for periodically stratified media is outlined, the Floquet wavenumber matrix **K** introduced and its Magnus expansion examined (see also appendix A, electronic supplementary material). The main results of the paper are derived in §3. It starts by observing that the matrix i**K** viewed as a sextic-system matrix **Q**_{eff} for a homogeneous effective medium cannot be associated with a medium from the class of anisotropic elastic materials but it does fit the Willis model. Using the *ansatz* that the sought effective medium is described by the spatially homogeneous equations for a Willis material, the corresponding system matrix **Q**_{eff} is constructed and equated with i**K**. Under certain assumptions, a prescription for unique definition of the Willis effective medium is put forward. The remainder of §3 discusses the general properties of the Willis material parameters. The example of SH wave motion in a periodic structure is considered in §4. It illustrates the method for defining the effective coefficients of the Willis model beyond the Magnus series expansion (which in its turn is detailed for SH waves in appendix B, electronic supplementary material). The explicit expressions obtained are used to solve a reflection–transmission problem at the interface of the effective medium. Conclusions are presented in §5.

## 2. Background

### (a) Stroh formalism and the wavenumber matrix *K*

*K*

We consider a Cartesian elastic medium with density *ρ*=*ρ*(*y*) and stiffness tensor *c*_{ijkl}=*c*_{ijkl}(*y*). Basic notations used include the superscripts ^{T}, ^{+} and * for transposition, Hermitian and complex conjugation, respectively, and **T** for the matrix with zero diagonal and identity off-diagonal blocks.

Taking the Fourier transforms of the equilibrium and stress–strain equations
2.1in all variables except *y* leads to an ordinary differential problem for the quasi-plane modes with the phase factor *e*^{i(kxx−ωt)}, where *ω* is the frequency and *k*_{x} the wavenumber in an arbitrarily chosen direction *X* orthogonal to *Y* (rotating *X* causes all 21 elastic constants *c*_{ijkl} to appear). Denote the unit vectors parallel to *X* and *Y* by **m** and **n**, so that *x*=**m**·**r**, *y*=**n**·**r** and let **A**(*y*) and **F**(*y*) be the amplitudes of displacement **u** and traction **n**** σ**, respectively. Then equations (2.1) combine into the Stroh system
2.2for the state vector

**incorporating**

*η***A**and

**F**(Ting 1996). Taking it in the form

**=(**

*η***A**,i

**F**)

^{T}defines the 6×6 system matrix as 2.3via the 3×3 blocks

**N**

_{J}of the Stroh matrix 2.4which is composed of the matrices with elements (

*nn*)

_{jk}=

*n*

_{i}

*c*

_{ijkl}

*n*

_{l}, (

*nm*)

_{jk}=

*n*

_{i}

*c*

_{ijkl}

*m*

_{l}=(

*mn*)

_{kj}and (

*mm*)

_{jk}=

*m*

_{i}

*c*

_{ijkl}

*m*

_{l}(note that

**N**

_{2}is negative-definite). The usual indicial symmetry of

*c*

_{ijkl}used in equation (2.4) leads to a Hamiltonian structure

**N**=

**TN**

^{T}

**T**and

**Q**=

**TQ**

^{T}

**T**of

**N**and

**Q**. Since

*ω*,

*k*

_{x}and

*ρ*,

*c*

_{ijkl}are real and hence

**Q**(

*y*) is imaginary, it follows that

**Q**=−

**TQ**

^{+}

**T**. The latter identity on its own suffices to ensure energy conservation. Alternative definitions of

**and hence of**

*η***Q**may be chosen. In general, the eigenvalues of

**Q**(

*y*) are first-degree homogeneous functions of

*ω*,

*k*

_{x}, which implies the absence of dispersion.

Given the initial condition at some *y*_{0}(≡0), the solution to equation (2.2) is ** η**(

*y*)=

**M**(

*y*, 0)

**(0), where**

*η***M**(

*y*, 0) is the 6×6 matricant evaluated by the Peano series (Pease 1965) 2.5Suppose now that

*ρ*,

*c*

_{ijkl}and hence

**Q**depend on

*y*periodically with a period

*T*. Denote 2.6where 〈 … 〉 is the static average over a unit cell. It is understood hereafter that the wave-path distance

*y*includes a large enough number

*n*of periods, which is when the present development is of interest. By equation (2.5), the matricant

**M**(

*T*, 0) over [0,

*T*], which is called the monodromy matrix, is 2.7

The wavenumber matrix **K** is introduced by denoting the monodromy matrix as
2.8In the following, unless otherwise specified, **K** is understood as defined in the first Brillouin zone, which implies the zeroth Riemann sheet of with a cut . Using equation (2.8), the matricant **M**(*y*, 0) can be written as
2.9where with . Equation (2.9)_{2} represents the Floquet theorem. Denote the eigenvalues of **M**(*T*, 0) and **K** by e^{iKαT} and *K*_{α} (*α*=1,…,6), respectively. In the general case where **M**(*T*, 0) and **K** have six linear independent eigenvectors **w**_{α}, the Floquet theorem implies that the wave field ** η**(

*y*)=

**M**(

*y*, 0)

**(0) with the initial data expanded as takes the form 2.10**

*η*The identity **Q**=−**TQ**^{+}**T** yields **M**^{−1}(*y*, 0)=**TM**^{+}(*y*, 0)**T**, which in turn leads to and
2.11for **K** defined in the first Brillouin zone. If the unit-cell profile is symmetric, i.e. the variation of material properties within the period *T* is symmetric about the middle point so that is even about then the above identities are complemented by **M**(*T*, 0)=**TM**^{T}(*T*, 0)**T** and hence **K**=**TK**^{T}**T**; so, with reference to equation (2.11), **K** is real. Thus
2.12

### (b) Expansion of *K* in the Magnus series

*K*

The logarithm of the monodromy matrix **M**(*T*, 0) can be expanded as a Magnus series (Blanes *et al.* 2009; see also Wang & Rokhlin 2004):
2.13where [**Q**(*x*),**Q**(*y*)]=**Q**(*x*)**Q**(*y*)−**Q**(*y*)**Q**(*x*) is a commutator of matrices depending on successive integration variables. Each Magnus series term **K**^{(m)} is a (*m*+1)-tuple integral of permutations of *m* nested commutators involving products of (*m*+1) matrices . A commutator-based form may be anticipated by noting that all **K**^{(m)} for *m*>0 must vanish in the trivial case of a homogeneous material with a constant **Q**≡**Q**_{0} and hence with . For practical calculations it is convenient to use the recursive formulas provided in Blanes *et al.* (2009). In obvious contrast with the Peano expansion, the Magnus series converges in a limited range: the sufficient condition for its convergence is 〈∥**Q**∥_{2}〉<*π*/*T*, where ∥·∥_{2} is the matrix norm (Moan & Niesen 2008). This condition implies that the eigenvalues *K*_{α}(*ω*^{2},*k*_{x}) of **K**, defined as continuous functions such that *K*_{α}(0, 0)=0, satisfy the inequality ∥Re *K*_{α}∥<*π*/*T*.

The Magnus series for **K** is a low-frequency long-wave expansion. Actually *ω* and *k*_{x} are two independent parameters for **Q** and hence for **M** and **K**. It is however essential that the dependence of 3×3 blocks of **Q** on *ω* and *k*_{x} is homogeneous (see equation (2.3)). Therefore, the blocks of each *m*th term of the Peano and Magnus series are homogeneous polynomials of *ω*, *k*_{x} of degree one greater than the same block of the (*m*−1)th term. This is what enables introducing a single long-wave parameter *ε*≡*kT* with a suitably defined wavenumber *k* (see equation (A3), electronic supplementary material). Thus, the Magnus series (2.13) is basically an expansion in powers of *ε*. Taking small enough *ε* enables its approximation by a finite number of terms. At the same time, it should be borne in mind that the Magnus series as an expansion of logarithm may converge relatively slowly. Explicit estimates expressed in terms of *ω*, *k*_{x} and 〈**N**〉 which ensure a desired accuracy of truncating the Magnus series up to a given order are formulated in appendix A, electronic supplementary material.

The structure of polynomial dependence of the Magnus series terms **K**^{(m)} on *k*_{x} and *ω*^{2} is
2.14where the real matrices **a**^{(m)} in **K**^{(m)} are (*m*+1)-tuple integrals of appropriate commutators; for instance, in **K**^{(1)} are
2.15The series terms **K**^{(m)} of odd and even order *m* are imaginary and real, respectively, and each term **K**^{(m)} on its own satisfies equation (2.11); therefore
2.16as taken into account in equation (2.14). According to equation (2.12),
2.17Note the pure dynamic imaginary terms proportional to ±i*ω*^{m+1}, which appear in the diagonal blocks of the series terms **K**^{(m)} of odd order *m* unless these are zero for a symmetric by equation (2.17).

### (c) Dynamic homogenization

According to the Floquet theorem (2.9), the wave field variation over a large distance *y* is characterized mainly by the function (which is an exact wave field at *y*=*nT*). Formally, with is a solution to equation (2.2) with the actual matrix of coefficients **Q**(*y*) replaced by the constant matrix i**K**. This motivates the concept of an effective homogeneous medium, whose material model admits the wave equation in the form (2.2) with a constant system matrix **Q**_{eff}≡i**K**. Confining the Magnus series (2.13) to the zero-order term defines the statically averaged **Q**^{(0)}_{eff}=〈**Q**〉, which fits equation (2.3) with **N**_{eff}=〈**N**〉 and hence yields the non-dispersive effective density and stiffness in the well-known form *ρ*^{(0)}=〈*ρ*〉 and
2.18see Norris (1992). It is evident that the statically averaged completely ignores the dynamic effects and is inadequate to describe waves at finite frequency over long propagation distance. Dynamic properties are realized by taking **Q**_{eff}=i**K** beyond the zero-order term 〈**Q**〉 (§3). Note that, in contrast to 〈**Q**〉, a dispersive **Q**_{eff}=i**K** generally depends on where the reference point *y*=0 of the period interval [0,*T*] is chosen.

## 3. A dispersive effective medium with **Q**_{eff}=i**K**

### (a) The constitutive equations

Our purpose is to take into account the full nature of the wavenumber matrix in **Q**_{eff}=i**K**. With this in mind, compare the structure of the dispersive effective matrix as given by the Magnus expansion, with that of **Q**(*y*) given by equation (2.3). They differ in two ways. First, **Q**_{eff} (≠〈**Q**〉), is no longer imaginary, and hence the identity, , which leads to equation (2.11) (and ensures energy conservation), is no longer compatible with a Hamiltonian structure for **Q**_{eff}. This is a well-known feature of dispersive models (e.g. Portigal & Burstein (1968)). The second, more significant, dissimilarity is that, in contrast to equation (2.3), **Q**_{eff} has pure dynamic terms on the diagonal, already at the first-order i**K**^{(1)} (see equation (2.14)). This means that assuming dispersive density and elastic constants does not suffice for the constitutive relations of the dispersive effective medium to be written in the standard form of equations (2.1). Recalling that the upper rows of the sextic system (2.2) imply the traction-strain law, it is seen that the latter must be complemented by a purely dynamic term, which implies different constitutive relations than those of the inhomogeneous medium itself.

On this basis, following the studies of Milton & Willis (2007) and Willis (2009), the equations of equilibrium and the constitutive relations (2.1) are replaced by the more general form proposed by Willis
3.1The vector **p** generalizes the normal notion of momentum density, and the elements of the Willis coupling tensor satisfy *S*_{ijk}=*S*_{jik} by assumption, ensuring the symmetry of the stress tensor. A principal objective is to show that setting **Q**_{eff}=i**K** leads inevitably to dispersive effective matrix density *ρ*^{(eff)} and stiffness and, on top of that, to the Willis form of the effective constitutive relations with stress-impulse coupling.

### (b) The effective Willis medium

Denote by **S**_{n} and **S**_{m} the matrices with components
3.2The same derivation that led from equation (2.1) to the sextic system (2.2) with the coefficients (2.3) now leads from equations (3.1) to (2.2) with the system matrix
3.3with . The identity is assumed in order to ensure that the effective medium is, like the inhomogeneous periodic medium, non-dissipative (energy conserving). This implies hermiticity constraints on the complex-valued material parameters: **c**^{(eff)}=**c**^{(eff)+}, *ρ*^{(eff)}=*ρ*^{(eff)+}, , , where **c**^{(eff)} is the 6×6 stiffness matrix in Voigt’s notation. The coupling tensor *S*_{ijk} is therefore purely imaginary. Note that the blocks of the effective Stroh matrix (≠**TN**^{T}_{eff}**T**) consist of sub-matrices (*nn*)^{(eff)},(*nm*)^{(eff)},(*mm*)^{(eff)} built from **c**^{(eff)} according to the definition (2.4) but with (*mn*)^{(eff)}=(*nm*)^{(eff)+} so that .

Equating the matrix **Q**_{eff} introduced in equation (3.3) to the matrix i**K**(*ω*^{2},*k*_{x}) with the block structure (2.11) yields the blockwise equalities
3.4Identification of the effective parameters based on these identities is ambiguous given that they may depend upon both *ω* and *k*_{x}. A unique and arguably the simplest solution is obtained by first assuming a purely dynamic *S*_{ijk}=*S*_{ijk}(*ω*). Then (3.4)_{1} and (3.4)_{2} yield
3.5where the blocks **K**_{J}(*ω*^{2},*k*_{x}) of **K** are taken at *k*_{x}=0 as indicated by the notation **K**_{J}(0)≡**K**_{J}(*ω*^{2}, 0) used here and subsequently. Note that the validity of equation (3.5) requires the block **K**_{2} to be invertible, which is assumed in the following. For **K** defined by the Magnus series (2.13), the conditions that ensure existence of in a certain low-frequency long-wave range and the estimates that enable truncation of its expansion are established in appendix A, electronic supplementary material. In view of equation (3.5), the remaining identity (3.4)_{3} becomes
3.6Pursuing the logical extension of the assumed dependence of the inertial coupling tensor **S** on frequency alone, we assume that the inertia tensor *ρ*^{(eff)} is also purely dynamic. This leads to a unique solution of equation (3.6) since *ρ*^{(eff)}=*ρ*^{(eff)}(*ω*) is found by setting *k*_{x}=0,
3.7The limit may be achieved in terms of derivatives of matrices **K**_{J} at *k*_{x}=0, whose existence is guaranteed for instance within the range of convergence of the Magnus expansion. Accordingly, the solutions of (3.6) are
3.8where .

In summary, equations (3.5) and (3.8) provide unique material properties for the Willis model with
3.9and *S*_{ijk}=*S*_{jik}. The lack of dependence of the inertial parameters on *k*_{x} means that non-local effects are confined to the elastic moduli **c**^{(eff)}. The result reduces to non-dispersive statically averaged moduli of equation (2.18) with *ρ*^{(eff)}=〈*ρ*〉**I** and **S**=**0** when i**K** is restricted to the zero-order term 〈**Q**〉 of equation (2.13).

Regarding computation of the Willis parameters from equations (3.5) and (3.8), it is assumed that the wavenumber matrix **K**(*ω*^{2},*k*_{x}) defined as is known either in the form of long-wave low-frequency series, or from the direct definition of matrix logarithm in some neighbourhood of a given point *ω*, *k*_{x} (see the example in §4). In particular, one may first evaluate the matricant **M**(*T*, 0) numerically and **K** then follows from the matrix logarithm. The matrix **K**^{′}(0), required for the solution of equations (3.8)_{2}, involves evaluating the derivative of in *k*_{x}. It may be expressed using either series or integral representation of as
3.10which leads to
3.11where **A**=**M**(*T*, 0)−**I** at *k*_{x}=0, and the derivative of the matricant itself is (Pease 1965)
3.12The sufficient conditions for the range of validity of the above series and integral definitions of **K**^{′}(0) are specified in appendix A, electronic supplementary material.

### (c) Discussion

#### (i) The Willis equation and its inertial quantities

Consider the above results (3.5) and (3.8) in more detail. Anisotropic density and the coupling coefficients that relate particle momentum and stress are unknown in ‘standard’ models of solids. Here they appear as inevitable ingredients of a model that replaces periodic spatial inhomogeneity with a spatially homogeneous but dispersive and non-local theory. The departure from normal elasticity is evident from the equations of motion for the displacement that follows from equation (3.1), 3.13This in turn leads to an energy conservation equation of the form , where the real-valued energy density and flux vector are 3.14

In order to gain some insight into these new dynamic terms, consider a layered transversely isotropic medium with the principal axis along **n** identified as the 2-direction. The effective density is of the form and the only non-zero elements of the coupling tensor (up to symmetries *S*_{ijk}=*S*_{jik}) are *S*_{112}=*S*_{332}, *S*_{211}=*S*_{233} and *S*_{222}. Only one combination of the three independent coupling elements has impact on the equations of motion,
3.15and it has no influence on pure SH wave motion (polarized in the plane orthogonal to **n**). Long-wave expansions of and *S*_{ijk} are presented in equation (3.16). Further detailed discussion for SH waves is provided in §4.

More generally, the absence of generating functions for means that some elements of the Willis coupling tensor *S*_{ijk} should be set to zero in order to complete its definition. The relevant elements are necessary to determine stress and momentum but do not enter into the equation of motion (3.13) and the sextic system (2.2) with (3.3) because the purely imaginary property of the coupling tensor means that *m*_{i}(*S*_{ijl}−*S*_{ilj}) are the elements of . Consider the two-dimensional situation with indices taking only two values so that, on account of the symmetry *S*_{ijk}=*S*_{jik}, there are at most six independent elements. Four of these may be found from equation (3.5)_{1}, and one more follows from equation (3.8)_{2} using the symmetry property. The single element **m**·**S**_{m}**m** is undefined and may be set equal to zero. In the three-dimensional situation all but four combinations of the 18 independent elements of *S*_{ijk} are obtainable. Taking an orthonormal triad {**m**_{1},**n**,**m**_{2}}, the following elements of the coupling tensor are not defined by the effective medium equations and are therefore set to zero: **m**_{α}·**S**_{mα}**m**_{β}+**m**_{β}·**S**_{mα}**m**_{α}, *α*,*β*∈{1,2}. To be explicit, let **n** lie in the 2-direction, then the **S**_{n} equation (3.5)_{1} defines the nine elements *S*_{2jk}; these combined with the equations (3.8)_{2} yield *S*_{112}, *S*_{132}, *S*_{32}, and the equations with the above prescriptions give *S*_{113}=−*S*_{131}, *S*_{313}=−*S*_{331}, *S*_{111}=0, *S*_{333}=0.

#### (ii) Expansion of the Willis parameters

Explicit insight into the structure of the Willis parameters can be gained from their expansion obtained via the Magnus series for the wavenumber matrix **K**. In view of equations (2.14) and (2.16), **K**_{1}(0) is imaginary and expands as with odd *m* and while **K**_{2}(0) is real and expands as with even *m* and . This confirms that **S**_{n} is imaginary and vanishes at *ω*=0. It is easy to check that the right-hand sides of equations (3.5)_{3} and (3.6) are zero at *k*_{x}=0 and at *ω*,*k*_{x}=0, respectively. Based on the forms of **K**_{J}(*ω*^{2},*k*_{x}) as generated by the Magnus expansion, and evident from equation (2.14) for the leading order contributions, it may be demonstrated that equation (3.8)_{2} is consistent with imaginary **S**_{m} that is zero at *ω*=0. It is noteworthy that keeping the density *ρ*^{(eff)} as 〈*ρ*〉 or as any other scalar is generally not possible as this would contradict the pure dynamic term on the right-hand side of equation (3.6). Finally, it is emphasized that, by virtue of equation (2.17), the stress-impulse tensor defined as a pure dynamic quantity *S*_{ijk}=*S*_{ijk}(*ω*) vanishes in the case of a unit cell with any symmetric heterogeneity profile (regarding the ‘inaccessible’ part see §3*c*(i)).

Equations (3.5) and (3.8) with polynomials **K**_{J}(*ω*^{2},*k*_{x}) given by the Magnus series (2.13) imply that the elastic moduli are rational functions of *ω*^{2},*k*_{x} while the density and coupling terms *ρ*^{(eff)}, *ωS*_{ijk} are functions of *ω*^{2}, defined by the series
3.16with real *ρ*^{(m)} and imaginary proportional to *ω*^{m}, real or imaginary depending on whether *m* is odd or even, respectively. These series are similar to equation (2.13) in that they are majorized by the power series in long-wave parameter *ε*. The Magnus series with *M* terms enables finding *M* terms of the series (3.16). It is apparent from equations (3.5)_{1}, (3.8)_{1} and (2.3) that **S**_{n} and *ρ*^{(eff)} depend only upon **N**_{2} and *ρ*, thus
3.17with
3.18in which are suppressed (as kept tacit hereafter) and dependence of co-factors on the successive integration variables is understood. The remaining part of the coupling tensor is only obtainable through the combination of **S**_{m}+**S**^{+}_{m}, and it depends upon **N**_{1}, **N**_{2} and *ρ*, with
3.19

For example, the expansions of the inertial quantities for the transversely isotropic layered medium discussed in §3*c*(ii) are *ρ*^{(eff)}=〈*ρ*〉**I**+*ρ*^{(2)}+…
3.20and , have, respectively, the same form as , with *c*_{22} replaced by *c*_{66} in equation (3.20)_{1,2}.

#### (iii) Effective medium defined from the Floquet dispersion

Modelling a dispersive effective medium may be based on a more relaxed approach that abandons fitting the matrix i**K** to the coefficients of sextic system of wave equations and deals instead with the asymptotic secular equation for the eigenvalues i*K*_{α} or *e*^{iKαT} of i**K** or **M**(*T*, 0), which is a dispersion equation for the onset of fundamental Floquet branches *K*_{α}(*ω*,*k*_{x}) or *ω*_{α}(*k*_{x},*K*) analysed in Norris (1992, 1993) and Norris & Santosa (1992). This gives the same secular equation as that for the i**K** matrices, and hence preserves the long-wave Floquet dispersion but not the displacement-traction vector **w**_{α} at the period edges (see equation (2.10)). By not fitting all of the physical properties, this type of approach to homogenization modelling introduces extra degrees of freedom. In particular, a ‘modified’ effective medium may be defined that is asymptotically similar to i**K** but the matrix has no pure dynamic terms in the diagonal blocks, and hence matches the Stroh-like forms (2.3) and (2.4) (though now with (3.9)_{1}), i.e. satisfies the standard form of the governing equations (2.1) with dispersive effective coefficients.

For instance, in the one-dimensional case *k*_{x}=0, the matrix
3.21has asymptotically (to the order of this matrix itself) the same secular equation as the matrix
3.22The latter ‘skips’ (by construction) the Willis coupling tensor and leads to the same definition of the matrix of second-order elastic coefficients (*nn*)^{(2)} as in equation (3.5), while the second-order density matrix following from equation (3.22) is generally different from *ρ*^{(2)} in equation (3.17) owing to non-commutativity of 〈**N**_{2}〉 and . See also the SH example in appendix B, electronic supplementary material.

## 4. Effective medium coefficients for SH waves

### (a) The wavenumber matrix

Consider SH waves in an isotropic medium with periodic density *ρ*(*y*) and shear modulus *μ*(*y*). The SH state vector ** η**(

*y*)=(

*A*,i

*F*)

^{T}, where

*A*and

*F*are the amplitudes of

*u*=

*u*

_{3}and

*σ*

_{23}(the indices correspond to

**u**∥

*X*

_{3},

**n**∥

*X*

_{2},

**m**∥

*X*

_{1}), satisfies equation (2.2) with the system matrix 4.1The 2×2 case leads to some simplifications not available for higher algebraic dimensions. In particular, the two eigenvalues of the monodromy matrix

**M**(

*T*, 0), which are the inverse of one other (since due to the isotropy), are defined by the single quantity

*tr*

**M**(

*T*, 0). The implications are explored in Shuvalov

*et al.*(2010) and only the necessary equations are cited here. The main result is that the wavenumber matrix, and hence the effective system matrix

**Q**

_{eff}(

*ω*)=i

**K**has semi-explicit form, and 4.2where , and ±

*K*(no subscript) are the eigenvalues of

**K**.

### (b) Willis equations and effective coefficients

Following the general formalism of §3, the effective material is assumed to have constitutive equations described by the Willis model, which in this case has only a single momentum component *p*_{3} and the usual stress components for SH waves in elasticity. Noting that *S*_{53}=0, on account of the transversely isotropic axis **n**, we have
4.3These constitutive relations imply, using equation (3.1)_{1}, that the governing equation for the SH displacement is of the form
4.4where ^{′} means d/d*y*. The coupling term *S*_{43} is absent from the equation of motion, as expected from the Willis equations (3.1) for a scalar problem. At the same time, equation (4.3) leads to the state-vector system matrix in the form
4.5where has been used.

Setting **Q**_{eff} of the Willis model equal to that of equation (4.2) gives the material parameters
4.6with *K*_{J}(0)=*K*_{J}(*ω*, 0). These may be expressed directly in terms of the elements of the monodromy matrix, using the form (4.2) along with ,
4.7where *M*_{J}=*M*_{J}(*T*, 0) are functions of *ω* and *k*_{x}, and (0) means evaluated at *k*_{x}=0. Note that the expressions for *ρ*^{(eff)} and also follow from the equation of motion (4.4) and its solution *u*(*y*)=*u*(0)*e*^{iKy}, using *k*_{x}=0 for *ρ*^{(eff)}(*ω*).

Explicit expressions for the low-frequency long-wave expansion of the material parameters may be found in the same manner as in §3*c*(ii) for the general case. The starting point is the Magnus expansion **Q**_{eff}=〈**Q**〉+i**K**^{(1)}+i**K**^{(2)} for the SH wavenumber matrix. Details of the analysis and a summary of the results are presented in appendix B, electronic supplementary material.

### (c) Examples and discussion

#### (i) A bilayered unit cell

The general formulation is illustrated by the case of a two-component piecewise constant unit cell. Specifically, consider a periodic structure of homogeneous isotropic layers *j*=1,2, each with constant density *ρ*_{j}, shear modulus *μ*_{j} and thickness *d*_{j}. The monodromy matrix has the well-known form
4.8where is the phase shift over a layer and *γ*_{j}=*μ*_{j}*ψ*_{j}/*d*_{j} (see Shuvalov *et al.* (2010)). Figures 1 and 2 show the computed parameters for the case of layers of equal thickness, *d*_{1}=*d*_{2}=1/2, with *ρ*_{1}=1, *c*_{1}=1; *ρ*_{2}=2, *c*_{2}=2, where *c*_{j} is the shear wave speed (*c*^{2}=*μ*/*ρ*). Figure 1 shows the effective parameters for propagation normal to the layers (*k*_{x}=0). The vanishing of both and *ρ*^{(eff)} at the band edge at *ω*=*ω*_{1}≈2.6 is expected on the basis of the fact that **Q**_{eff} is singular at the band edge and scales as (*ω*−*ω*_{1})^{−1/2} near it (Shuvalov *et al.* 2010). Referring to the 12-element in equation (4.5), this implies first that and then, from the 21-element and the finite value of , that *ρ*^{(eff)}∝(*ω*−*ω*_{1})^{1/2}. The square root decay of both and *ρ*^{(eff)} is apparent in figure 1.

The wavenumber is finite, *k*_{x}=1, in figure 2. This has the effect of increasing the frequency of the band edge, and introducing a range of frequency from *ω*=0 up to the cut-on at *ω*≈1.7 in which the effective wave is non-propagating. Note that *ρ*^{(eff)} and *S*_{43} are unchanged from figure 1 while the elastic modulus is different, and tends to zero at the new band edge as expected. The non-zero *k*_{x} leads to non-zero , and the parameter becomes complex-valued at the *k*_{x}=0 band edge. Only the real parts of the quantities are shown in both figures. No attempt is made here to discuss their imaginary components, which requires careful analysis of the branch cuts and is a topic for separate study.

#### (ii) Reflection and transmission of a half-space of effective material

As an example of the type of boundary problem that can be solved using the effective medium equations, consider reflection–transmission of SH waves at a bonded interface *y*=0 between the half-space of the periodically stratified medium (*y*>0) and a uniform half-space (*y*<0) of isotropic material with *ρ*_{0}, *μ*_{0} and . An SH plane wave is incident from the uniform half-space with propagation direction at angle *θ* from the interface normal. The total solution is taken as
4.9The reflection and transmission coefficients *R* and *T* follow from the continuity conditions for particle velocity and traction at the interface. They may be expressed in the standard form using SH impedances defined as . The impedance in the uniform half-space is . The impedance *Z*_{+} for the effective medium follows from equation (4.3) as
4.10This is identical to the impedance of the periodically stratified half-space because they both imply a ratio of components of the outgoing eigenvector **w**, which is common to **M**(*T*, 0) and **K**. In these terms, the continuity conditions for displacement and traction yield the exact result
4.11Figure 3 shows |*R*(*ω*)| and |*T*(*ω*)| calculated for normal incidence from a uniform half-space with *ρ*_{1}=1, *c*_{1}=1 on a periodic structure of two layers with *ρ*_{1}=1, *c*_{1}=1, *ρ*_{2}=2, *c*_{2}=2, which was used in figures 1 and 2. As expected, |*R*|≤1 with total reflection in the stopband.

The explicit dependence of the reflection coefficient on the effective medium parameters , and *S*_{43} means that, in principle, measurement of *R* via experiment can provide useful knowledge for their determination.

#### (iii) Uniform normal impedance

It is instructive to consider the particular case of *k*_{x}=0 with independent of *y*, i.e. *z*≡*z*_{0}. The 2×2 matrix **Q**(*y*) is then a scalar multiple of a constant matrix, and *M*_{1}(0)=*M*_{4}(0), (where (0) stands for *k*_{x}=0). As a result, by equation (4.6), the effective parameters at any *ω* retain their values obtained from static averaging: *S*_{43}=0, and *ρ*^{(eff)}(*ω*)=〈*ρ*〉. This simplification is a consequence of the fact that constant *z* implies constant eigenvectors of the SH matricant **M**(*y*, 0) and hence no reflection of SH waves normally propagating through a periodic structure, which is in accordance with the physical meaning of the impedance *z*. Consistency is also observed in that the effective impedance *z*^{(eff)} defined through the above effective parameters is equal to *z*_{0},
4.12The only effect of the inhomogeneity is to speed up or retard the advancing waves according to the effective speed *c*^{(eff)}=*z*^{(eff)}/*ρ*^{(eff)}, which in the present case follows from equation (4.12) as *c*^{(eff)}=*z*_{0}/*ρ*^{(eff)}=〈*c*^{−1}〉^{−1}.

#### (iv) Discussion

In the case of purely unidimensional motion, *k*_{x}=0, the system (4.3) involves only the first three parameters of equation (4.7): , *ρ*^{(eff)} and *S*_{43}. Willis (2009) derived expressions for the same quantities for a laminated medium. His coefficients (Willis 2009, eqn. (3.30)) relate weighted means of strain and velocity (〈*we*〉, to ensemble the means of stress and momentum density (〈*σ*〉, 〈*p*〉), where *w* is a general weighting function first introduced in Milton & Willis (2007), such that the ensemble means correspond to *w*=1. The Willis parameters derived here, e.g. equation (4.3), concern strain and velocity at the single point *y*=0 in the unit period, and therefore correspond to the specific weight function *w*(*x*)=2*Lδ*(*x*) in the notation of Willis (2009). It is important to note, however, that the stress and momentum density used here are not ensemble averages but are quantities associated with the same point in the unit period. This identification, for instance, means that the solution of the reflection–transmission problem of §4*c*(ii) is in fact the deterministic solution. In summary, while the governing equations are the same in both cases, the Willis parameters developed here do not bear a one-to-one correspondence with those in Willis (2009).

## 5. Conclusion

A fully dynamic homogenization scheme has been developed for periodically layered anisotropic elastic solids. In the process, the dispersive and non-local Willis model has been shown to provide an optimal constitutive setting for the effective medium. The crucial point of the present method is the insistence that the matrix of coefficients **Q**_{eff} of the sextic system of elastodynamics equations, for whatever homogeneous effective medium is considered, must exactly match the Floquet wavenumber matrix **K** of the periodic system. This is not a low-frequency long-wave approach, so long as **K** is defined at the given frequency *ω* and horizontal wavenumber *k*_{x}. The wavenumber matrix **K**=**K**(*ω*,*k*_{x}) is an analytical function of *ω* and *k*_{x}, which may be explicitly defined via the Magnus series expansion that is guaranteed to converge below the first Floquet stopband at the edge of the Brillouin zone. The choice of constitutive model for the effective medium is critical. We have demonstrated that the standard anisotropic elasticity theory does not suffice as it cannot provide a **Q**_{eff} to properly account for dynamic terms appearing in the wavenumber matrix, **K**(*ω*,*k*_{x}). On the other hand, the Willis model for the effective medium, which includes coupling effects, can allow us to associate elements of the effective system matrix **Q**_{eff} with elements in **K**. The main results are contained in equations (3.5) and (3.8), which infer the material parameters of the effective Willis medium from **K**. Invoking the Magnus series, explicit expressions for the low-frequency long-wave expansion of these effective Willis parameters have been found, and the accuracy for their truncated asymptotics has been estimated.

The example of SH plane wave reflection and transmission considered in §4*c*(ii) indicates the type of application possible using the Willis effective medium. The point is not so much to provide new solutions for layered media, although it is simpler to formulate and solve such problems using equations for a homogeneous model. The potential power of the dynamic effective medium model is that it is possible to relate the effective properties of the Willis material to the measurable dynamic quantities. Thus, the reflection and transmission problem illustrates how the reflection coefficient *R* depends on a certain combination of the Willis parameters. Measurements of *R*=*R*(*ω*,*k*_{x}) provide a means to characterize periodic layered systems as equivalent homogeneous but dispersive materials. Other problems that may be considered are, for instance, surface wave propagation in a periodically layered half-space, waveguides comprised of periodic layers and point forces.

## Acknowledgements

This work has been supported by the grant ANR-08-BLAN-0101-01 from the ANR (Agence Nationale de la Recherche) and by the project SAMM (Self-Assembled MetaMaterials) from the cluster AMA (Advanced Materials in Aquitaine). A.N.N. is grateful to the Laboratoire de Mécanique Physique (LMP) of the Université Bordeaux 1 for the hospitality.

- Received July 22, 2010.
- Accepted December 15, 2010.

- This Journal is © 2011 The Royal Society