## Abstract

All metamaterial applications are based upon the idea that extreme material properties can be achieved through appropriate dynamic homogenization of composites. This homogenization is almost always done for infinite domains and the results are then applied to finite samples. This process ignores the evanescent waves which appear at the boundaries of such finite samples. In this paper, we first clarify the emergence and purpose of these evanescent waves in a model problem consisting of an interface between a layered composite and a homogeneous medium. We show that these evanescent waves form boundary layers on either side of the interface beyond which the composite can be represented by appropriate infinite domain homogenized relations. We show that if one ignores the boundary layers, then the displacement and stress fields are discontinuous across the interface. Therefore, the scattering coefficients at such an interface cannot be determined through the conventional continuity conditions involving only propagating modes. Here, we propose an approximate variational approach for sidestepping these boundary layers. The aim is to determine the scattering coefficients without the knowledge of evanescent modes. Through various numerical examples we show that our technique gives very good estimates of the actual scattering coefficients beyond the long wavelength limit.

## 1. Introduction

Metamaterials are artificially designed composite materials which can exhibit properties that are not found in naturally occurring materials. These properties can be electromagnetic [1,2], acoustic [3–6] or elastodynamic [7]. In the context of electromagnetism, these properties refer to magnetic permeability and electrical permittivity. For acoustic metamaterials, they refer to bulk modulus and density and for elastodynamic metamaterials they refer to moduli (bulk, shear and anisotropic) and density. Irrespective of the different properties which metamaterials research in different fields target, their final aim is the same. Metamaterials research seeks to design composite materials for the fine-tuned, predominantly frequency-dependent control of the trajectory and dissipation characteristics of the applicable waves.

Metamaterial properties are generally achieved through an appropriate dynamic homogenization technique which relates a microstructure to its frequency-dependent homogenized properties. Currently, there are two main ways of doing this. The first is based upon asymptotic methods [8–16] and the second is based upon field averaging methods [17–27]. In addition to these, there are also scattering measurements based analytical and experimental techniques. While simple, in principle, these run the risk of resulting in properties which violate basic thermodynamic laws [2,28]. Within the established routes of dynamic homogenization, the process is the following: periodic boundary conditions are assumed over a unit cell which gives rise to wave solutions of the Bloch form. The fields resulting from these Bloch waves are then homogenized which gives rise to frequency-dependent effective properties for the composite metamaterial.

Dynamic homogenization serves as a route for realizing the challenging properties required by the application areas of metamaterials research (transformation acoustics [29], elastodynamics [30,31], etc.) The assumption is that the regions which require a certain set of metamaterial properties can be realized through a composite whose dynamically homogenized properties are the same as the desired properties. A deeper assumption here is that the homogenized properties which were initially calculated for infinite domains can now be applied to non-infinite domains. This assumption is not always correct. The free space homogenized properties may or may not apply to non-infinite cases and this has been explicitly shown to be the case by various researchers [32–34]. The reason for the failure of this assumption is subtle. As dynamically homogenized properties are calculated for free space waves, they do not allow for any evanescent modes. However, the composites which are supposed to realize these properties in a finite setting support evanescent modes. When such composites are interfaced with other regions then these evanescent modes are integral in satisfying displacement and stress continuity across the interface. Without these evanescent modes, however, the stress and displacement fields are discontinuous. Therefore, if a complete correspondence is to be maintained between the composite and the dynamically homogenized region which it is supposed to represent (in the sense that scattering from the two should be equivalent) then the relationships between the displacement and stress fields at the interface between the homogenized region and its surrounding region are not ones of simple continuity. In other words, the interface conditions are indeterminate and the true scattering coefficients cannot be determined through simple displacement and stress continuity relations. This is clearly a fundamental issue which results from ignoring evanescent modes in the process of dynamic homogenization.

This issue has been recognized in the electromagnetics community for several decades. There have, therefore, been efforts to account for the effect of these evanescent waves. The predominant technique for doing so is through the inclusion of the so-called Drude transition layers [2,35,36]. A transition layer is an artificial layer which is placed between two materials (a metamaterial and a homogeneous material for instance) and which has suitably chosen material properties. The choice is geared towards producing the scattering coefficients which would have resulted had the evanescent modes not been ignored in the original problem. The concept of introducing a transition layer with simple physical properties to simulate the effect of the actual boundary layers created by evanescent waves is an appealing one but defining how to obtain a suitable set of parameters is not straightforward except in the ‘quasi-static’ or ‘homogenization’ range of frequencies, for which it was first introduced. Furthermore, it is not clear how these transition layers can be made to handle cases where there are multiple propagating modes. More recently, efforts have been made to account for this issue through the use of asymptotic methods [37]. However, these techniques are applicable only in the long wavelength limit and their extension to frequencies far beyond the limit would be open to question.

In this paper, we propose a different approach for taking into account the effect of the evanescent waves. The approach is illuminated through its application to a model problem which has been considered in recent papers [38–40]. It is based upon the observation that at an interface, the energy flux balance only contains contributions from the propagating modes. We, therefore, propose to determine the scattering coefficients of the propagating modes by insisting that they exactly satisfy the energy flux balance while minimizing the displacement discrepancy across the interface. Thus, the energy flux is exactly satisfied but the continuity conditions are only approximately satisfied. We show that the process works very well in estimating the ‘exact’ scattering coefficients for a wide range of cases including those which are far beyond the long wavelength limit. The process is also general enough to potentially apply to more complex two- and three-dimensional cases.

## 2. Bloch waves in the laminate

Following [38], we define our laminate as a periodically layered structure in the *x*_{1}-direction with the layer interfaces parallel to the *x*_{2}–*x*_{3} plane. In the direction of periodicity, the laminated composite is characterized by a unit cell *Ω* of length *h* (0≤*x*_{1}≤*h*). For our purposes, the unit cell is composed of two material layers with shear moduli *μ*_{1},*μ*_{2}, densities *ρ*_{1},*ρ*_{2}, and thicknesses *h*_{1},*h*_{2}, respectively. If anti-plane shear waves are propagating in the laminate, then the only non-zero component of displacement is taken to be *u*_{3} which has the functional form *u*_{3}(*x*_{1},*x*_{2},*t*). Within the *i*^{th} layer (*i*=1,2), it satisfies the following equation of motion:
*σ*_{13}(*x*_{1},*x*_{2},*t*),*σ*_{23}(*x*_{1},*x*_{2},*t*). The shear stress component *σ*_{13} and displacement *u*_{3} are continuous at the material interfaces. Across an interface between layers *i* and *i*+1 at *x*_{1}=*x*^{i}:
*Ω*; specifically, for the displacement field, we have
*k*_{2} must be continuous across the layers to satisfy Snell’s Law. The other non-zero stress component *σ*_{23} has a similar Bloch-periodic form but *x*_{2},*ω* dependence suppressed):
^{−iK1h}. Quantities in the above equation depend upon assumed values of *ω*,*k*_{2}. The solutions to the eigenvalue problem above furnish the wavenumber *K*_{1} and the modeshape for which (2.3) satisfies the governing equation. The wavenumber solutions themselves come from the following equation:
*K*_{1} is a solution then so are ±(*K*_{1}±2*nπ*/*h*) for all integer *n*. Consider two different solutions of the current problem (details in [38]):
*μ*:

## 3. Normal mode decomposition

We now consider an interface between a homogeneous medium with shear modulus *μ*_{0} and the layered composite. The interface itself can be placed at any angle with the layers but presently we assume that it is along *x*_{2}=0 (figure 1). The layered medium is in the region *x*_{2}>0 with the layers being parallel to the *x*_{2}-axis. A plane harmonic wave is incident at the interface from the homogeneous medium. This wave sets up an infinite number of transmitted and an infinite number of reflected waves. A finite number of these are propagating waves and the rest are evanescent waves. The incident, transmitted and reflected fields are written down as
*U*_{n}(*x*_{1})=e^{−i2nπx1/h} and *x*_{2}<0. The wavenumber components *M* propagating modes and infinitely many evanescent modes in the *x*_{2}-direction. The problem of determining the scattering parameters can be solved to any required degree of precision by considering a sufficiently large number of terms in the normal mode expansions. To facilitate calculations, we can restrict the reflected modes to a range of −*N*≤*n*≤*N* and transmitted modes to a range of 0≤m≤2*N* such that 2*N*+1>*M*. This allows us to consider all propagating transmitted modes in the expansion. With this, the displacement (*u*_{3}) and stress (*σ*_{23}) continuity are given by (exponential terms suppressed):
*N*+1) equations in as many variables through the application of the orthogonality condition (2.8). Specifically, we have
*S* is a column vector of size 2(2*N*+1) with elements *M*_{i}] are square matrices of sizes (2*N*+1)×(2*N*+1) with the following non-zero elements:
*I* is a column vector of size 2(2*N*+1) with elements
*i*th mode.

## 4. Evanescent field as a boundary layer

To illustrate the role of evanescent waves, we consider the general example that was treated in [38]. The homogeneous medium is taken to be aluminium (*μ*_{0}=26 GPa, *ρ*_{0}=2700 kg m^{−3}) and the laminated composite is composed of two materials of thicknesses *h*_{1}=0.003 m (epoxy: *ρ*_{1}=1180 kg m^{−3}, *μ*_{1}=3 GPa) and *h*_{2}=0.0013 m (steel: *ρ*_{2}=8000 kg m^{−3}, *μ*_{2}=80 GPa). We take the frequency of excitation to be 200 kHz in which case the first propagating band is fully developed (figure 2). At this frequency, there is one fully propagating mode in the laminate and infinitely many modes which are propagating in the *x*_{1}-direction and non-propagating in the *x*_{2}-direction. For the interface configuration under consideration, it is the set of non-propagating modes which forms the evanescent boundary layer adjacent to the interface and on the side of the laminate. There is another evanescent boundary layer formed on the side of the homogeneous material which is determined by those reflected modes for which *κ*^{(n)} is negative and imaginary. Concentrating for the present on the laminate side of the interface, evanescent modes in the *x*_{2}-direction have the functional dependence *k*_{1}*h*=1.5, the first evanescent mode has *x*_{2}=1.4 mm from the interface (*t*_{e}, defined (for the purpose of this discussion) as the distance from the interface at which the evanescent wave with the smallest magnitude of its wavenumber reaches 10% of its amplitude at the interface. For the present case when *t*_{e}=0.7437*h*. At locations which are more than *t*_{e} away from the interface, the scattered field can be taken to consist only of the propagating solutions. On the homogeneous material side, the boundary layer thickness is similarly determined by that imaginary *κ*^{(n)} which has the smallest magnitude. The boundary layer thickness changes with both the angle of incidence and the frequency. It can be arbitrarily large at those frequencies and incidence angle combinations where the first evanescent mode has a vanishingly small imaginary part of the wavenumber. In metamaterial applications, one generally seeks to replace finite composite samples with their free space homogenized constitutive properties. Inherent in this process is the assumption that the evanescent modes which are invariably generated at an interface can be neglected away from it. It is clear from the above that even at low frequencies there may exist cases where no such replacement is possible because the influence of the evanescent modes may persist throughout the sample due to a large value of *t*_{e}.

### (a) Continuity conditions and energy conservation

While the evanescent modes are required for the satisfaction of the boundary conditions (3.4), they do not enter the energy conservation equation (3.12). Conservation of energy must, therefore, emerge from the satisfaction of boundary conditions. To show that this is indeed true, we split the energy conservation equation into the *m*_{t} transmitted and *n*_{r} reflected propagating components (since non-propagating components do not contribute):
*U*_{n},*U*_{n′}〉=*δ*_{nn′},

### (b) Boundary layers and the role of evanescent waves in satisfying continuity conditions

Within the boundary layers, the relevant stress and displacement components vary continuously with *x*_{2}, for any fixed *x*_{1}, between the two values which correspond to the free space propagating waves on both the reflected and the transmitted sides. To show this, we use the modified form of the interface conditions (4.2) and term the left and right sides of the displacement equation *u*_{t} and *u*_{r}, respectively. Similarly, *σ*_{t} and *σ*_{r} refer to the analogous *σ*_{23} stress components. Figure 3 shows the variation of the displacement and stress fields (absolute values) as a function of the distance from the interface. Positive values of distances are into the laminate and all calculations are done over a line *x*_{1}=*const*. which bisects any one of the steel laminae, material 2. We have used an angle of incidence of 30° and for this case the boundary layer thicknesses in the laminate and the homogeneous medium are 0.0084 m and 0.0053 m, respectively. To clarify the effect of the boundary layer, we have also plotted the absolute values of only the scattered components in the homogeneous medium (dashed curves.) Since this case only has one propagating transmitted mode and one propagating reflected mode, the absolute values of the scattered field beyond the boundary layer are constants. This is evident from the dashed curves on the homogeneous side and the solid curves on the laminate side wherein the absolute values of the scattered field stabilize to constants beyond the boundary layer. The sinusoidal variations of the absolute values of *u*_{r},*σ*_{r} are due to the combination of two waves with opposing wavenumbers in the *x*_{2}-directions. At the interface the purpose of the evanescent waves is to match *u*_{r} and *σ*_{r} with *u*_{t} and *σ*_{t}, respectively. We define [*u*]=*u*_{r}−*u*_{t} and [*σ*]=*σ*_{r}−*σ*_{t} as measures of how well the continuity conditions are satisfied. These are functions of *x*_{1} and are also dependent upon the magnitudes of the displacement and stress terms. To understand how well continuity conditions are being satisfied we average and normalize these measures as *a*,*b* plots *c*,*d* shows *a*,*b* shows why it is far easier to satisfy displacement continuity than stress continuity. For 200 kHz, 30° incidence, and 30 evanescent modes, we plot *u*_{r},*u*_{t},*σ*_{r},*σ*_{t} as a function of the location along the interface. While *u*_{t} is continuous along the unit cell, *σ*_{t} is discontinuous at the layer transitions. This is not surprising as while the displacement is required to be continuous across the layers, there is no such continuity requirement on *σ*_{23} across the layers. *σ*_{r} on the other hand will necessarily be continuous as it is composed of trigonomentric functions. With allowance for 30 evanescent modes, *u*_{t} and *u*_{r} are virtually indistinguishable but the match between *σ*_{t} and *σ*_{r} is visibly imperfect. To get a good match of *σ*_{r} with *σ*_{t}, it is clear that evanescent modes with high wavenumbers are required which can adequately approximate the discontinuous jump. On the other hand, only relatively few evanescent waves are required to match *u*_{r} with *u*_{t} which are both necessarily continuous functions. It is further notable that the coefficients for the propagating modes do not change significantly as the number of evanescent modes in the expansion is increased (figure 5*c*). This behaviour is seen for frequency and incidence angle combinations mentioned in figure 4*c*,*d*. From the above arguments and especially for higher angles of incidence, it is clear that the inclusion of a greater number of evanescent modes in the expansion will serve primarily to improve the stress continuity.

## 5. A variational approach for sidestepping the boundary layers

It is clear from the above treatment that beyond the boundary layers the scattered field on both sides of the interface can be safely taken to comprise only the propagating modes. Therefore, in these zones the material behaviour can be described by dynamically homogenized free space effective properties. The issue is the complicating presence of the boundary layers. Given two homogeneous materials and an interface between them the behaviour of a wave which impinges upon it is determined by the stress and displacement continuity relations across the interface. However, given a homogenized metamaterial, a homogeneous medium (or another homogenized metamaterial) and an interface between them, how should one determine the behaviour of an impinging wave? It is clear from figure 3 that stress and displacement are not continuous across the interface if one neglects the boundary layers on either side. The above question can be rephrased into the following: given only the propagating modes of the homogeneous and homogenized media and an interface between them, how can we determine the scattered field resulting from an impinging wave?

This question is of practical importance in the area of metamaterials research, the primary concern of which is to achieve extreme material properties through appropriate dynamic homogenization techniques. These homogenized properties are almost always calculated for free space propagating waves and then applied to finite or semi-infinite samples. The boundary effects, therefore, are generally neglected. The free space homogenized properties may or may not apply to non-infinite cases and this has been explicitly shown to be the case by various researchers [32–34]. In the area of electromagnetism, the boundary effect is generally taken into consideration through the inclusion of Drude transition layers [2,35,36]. The problem of assigning to such a layer a suitable set of parameters is not straightforward except in the ‘quasi-static’, or ‘homogenization’ range of frequencies, for which it was first introduced. Simovski [2] has reported progress in identifying layers suitable for higher frequencies, though his work considered only a dipole lattice approximation for electromagnetics. His proposed layers were aimed for use in what he termed the ‘metamaterial’ range of frequencies, at which the wavelength in the matrix material was significantly greater than the Bragg wavelength but the possibility of resonance of the dipoles could occur. He considered only normal incidence and suggested but did not prove that exactly the same parameters might be applicable also for oblique incidence. It is not clear how the Drude layer concept could be developed to accommodate frequencies at which there are more than one propagating transmitted and/or reflected waves.

We propose to determine the scattering coefficients (allowing only for propagating waves) through an indirect route. First we note that the energy conservation equation (4.1) only consists of contributions from propagating modes. Second we note from figures 4 and 5 that displacement continuity is satisfied more easily than stress continuity and that higher evanescent modes play a significant role in satisfying stress continuity. With these observations, we propose to solve the following minimization problem in search of the appropriate scattering coefficients:
*ϕ* is the energy constraint:
*m*_{t} propagating transmitted modes and *n*_{r} propagating reflected modes, then this system is expressed in a matrix form:
**M**,**N** are square matrices of size *m*_{t}+*n*_{r}, **I** is a column vector of length *m*_{t}+*n*_{r}, and λ is the Lagrange multiplier. We have
**N** is diagonal with components *m*_{t}+*n*_{r}+1 equations in *m*_{t}+*n*_{r}+1 unknown variables. Being nonlinear this system is solved through established gradient descent algorithms. In the following subsections, we compare the scattering coefficients which we calculate from the above minimization process with those which are calculated from equation (3.6). The former considers only propagating modes whereas the latter considers both propagating and evanescent modes.

Before proceeding, it is relevant to refer to figure 6, which is similar to fig. 2 of [38] but with slightly different ranges of frequency. Figure 6*a* displays equifrequency contours in the *K*_{1}–*k*_{2} plane in the lower range of frequencies. Except for *f*=270 kHz there is only one propagating mode, to which we assign the label *m*=0. These waves undergo positive refraction (both components of group velocity are positive). At *f*=270 kHz, however, there are two transmitted modes for *K*_{1}*h* greater than about 2.2. For incidence from an aluminium half-space, this corresponds to an angle of incidence of approximately 69.4°. This additional mode is negatively refracted and is assigned the label *m*=1. It appears first (at *K*_{1}*h*=*π*, *k*_{2}*h*=0) at *f*≈261332 Hz. Figure 6*b* displays a range of higher frequencies in which there are two propagating modes, *m*=0, positively refracted and *m*=1, negatively refracted. The range of *K*_{1}*h* values over which the mode *m*=1 exists increases as the frequency increases. The figure does not show it but the frequency above which there are two propagating modes for all *K*_{1}*h*∈[0,*π*) is 296 630 Hz.

### (a) Examples with positive refraction

As the first example, we compare the scattering coefficients calculated from the two approaches at 100, 200 and 260 kHz. At these frequencies and for all angles of incidence less than 90°, there exist one transmitted propagating mode (*T*_{0} mode) and one reflected propagating mode (*R*_{0} mode). Figure 7 shows the above-mentioned comparison for the absolute values of the scattering coefficients as functions of the angles of incidence. In this and later figures, the values calculated by allowance for 30 modes are plotted as solid lines, while those calculated from the optimization scheme are shown as diamonds. Note that there is no reason why the phase information (real and complex parts) of the scattering coefficients should match for the two approaches. This is due to the fact that the imposed constraint is on energy which depends only upon the magnitudes of the scattering coefficients. The definition of the phase for any wave in the laminate is arbitrary, in any case: if *θ*. The good agreement is not surprising in the case of the lowest frequency because 100 kHz is not far beyond what may be regarded as the ‘homogenization’ range, in which the evanescent modes contribute little, the transmission and reflection coefficients are real and the energy balance equation provides a legitimate substitute for the equation giving continuity of traction. It can be seen from the figure that optimized results are also close to the ‘exact’ results at the two higher frequencies, especially for incidence directions away from normal. They are, in addition, good enough to be useful, even close to normal incidence: the reflection coefficient is small and the reflected energy is proportional to the square of its magnitude. This is demonstrated in figure 7*d*, which shows plots of the approximate and ‘exact’ energy fluxes (the individual terms in (4.1)), for *f*=200 kHz. Although the relative error may be large, the absolute error is small. Thus, the scattered energies for both propagating modes are well estimated without considering the boundary layers by the variational scheme, at all angles of incidence.

### (b) Examples with negative refraction

Figure 8 shows the calculated comparisons for frequency *f*=285 kHz, at which there is a switchover from one to two transmitted modes, at an angle of incidence around 30°. The lower two of the plots shows the comparisons in terms of energy flux. Good performance of the optimization scheme is again demonstrated, even close to the angle of incidence at which switchover occurs. It should perhaps be noted that the wave *m*=1 exists but is evanescent at smaller angles of incidence, with a rate of decay that approaches zero as the switchover angle is approached, corresponding to the boundary layer becoming arbitrarily thick. Our variational approximation is thus severely tested around this angle of incidence. Figure 9 confirms good performance for *f*=300 kHz, at which frequency there are two propagating transmitted modes for all angles of incidence.

This section is concluded with some results for a different system, which was discussed by Srivastava [40]. The laminated material is the same but the aluminium half-space is replaced by one with shear modulus *μ*_{0}=0.4818 GPa and density 3000 kg m^{−3}. Incidence from this half-space generates a single-transmitted mode, which is refracted negatively for angles of incidence greater than 30°, and there are two propagating reflected modes. The relevant plots are shown in figure 10.

### (c) Interface parallel to the layers

Consider now the special case when the interface is parallel to the layers. The location of the interface is taken to be at *x*_{1}=*d* where 0≤*d*≤*h*. The relevant fields in this configuration are
*x*_{1} wavenumber component *K*_{1} and a corresponding modeshape *x*_{1}=*d* gives

#### (i) Energy flux continuity

The energy flux component in the laminate relevant to the present configuration is given by the *x*_{1}-component of the Poynting vector:
*x*_{1}. Given that *x*_{1} in the laminate, the balance of energy for the transmission/reflection problem is easily deduced to be
*K*_{1} is real.

Our variational approximation was designed so as to avoid dealing with evanescent waves that would appear in an exact formulation. In the present case, however, there are no such waves and the exact solution (5.9) is easily found. The minimum value of the (squared) magnitude of the displacement jump is zero and the energy balance is attained when *u*][*u*]*+[*σ*][*σ*]* subject to the energy flux constraint. It reproduced the coefficients calculated directly in figure 11. The behaviour at other frequencies is similar. In this figure, the steel lamina is positioned centrally. It is least noteworthy and perhaps surprising that the reflection and transmission coefficients are symmetric about the centreline. Thus, for instance, a composite whose first lamina is epoxy of thickness *h*_{1}/4 transmits exactly the same flux of energy as one whose first lamina is epoxy of thickness 3*h*_{1}/4. The mathematical reason is that the mode shape *T* itself is not. Physically, we simply observe that our calculation is one just for the steady state which remains once the transients have propagated away.

## 6. Concluding remarks

In this paper, we clarify the emergence and purpose of the evanescent waves which appear at the interface between a metamaterial and a homogeneous region in a model problem. We show that these evanescent waves form boundary layers on either side of the interface and that outside of these boundary layers the composite can be represented by appropriate infinite domain homogenized relations. We show that if one ignores the boundary layers, then the displacement and stress fields are not continuous across the interface. Therefore, the scattering coefficients at such an interface cannot be determined through the conventional continuity conditions involving only propagating modes. We propose an approximate variational approach for sidestepping these boundary layers. The aim is to determine the scattering coefficients without the knowledge of the evanescent modes. Through various numerical examples, we show that our technique gives very good estimates of the actual scattering coefficients, not only for the long wavelength region but far beyond it as well. The scattered energy is well estimated for all modes and at all angles of incidences—even in cases where multiple transmitted or reflected modes were present. The technique works well even in the case where negative refraction is occurring.

The example with the interface parallel to the layers provides a hint that minimizing a (suitably normalized) combination of [*u*][*u*]*+[*σ*][*σ*]* subject to the energy constraint could be useful for other problems: for instance, when the composite comprises a matrix with a two- or three-dimensional array of inclusions and the interface between the composite and the homogeneous medium is wholly between the matrix and the homogeneous medium, so that the traction components associated with the propagating modes would be continuous. Study of such possibilities would require more extensive computations than those reported here.

A further question that could be raised is this. Suppose the exact position of the interface relative to the microstructure is not known. What then would be the ‘best estimates’ for the transmission and reflection coefficients? The only certainty is the energy balance, and the variational scheme highlighted in this work would be expected at least to yield a unique approximation, so long as the exact formulation would contain an infinite number of evanescent waves. Its virtue, which is also its weakness, is that it need take no account of the properties of the material at the interface. The problem of the laminate with its boundary parallel to the layers provides an extreme example. The transmission and reflection coefficients show significant sensitivity to the exact thickness and composition of the first layer. For the frequency shown, adopting an approximation such as taking the average values would be prone to error but would remain indicative of the actual response. The error reduces to zero as the low-frequency homogenization limit is approached. The value of *K*_{1}*h* for the example shown was 2.504, demonstrating significant departure from that limit.

## Authors' contributions

Both authors contributed equally to the research and to the writing.

## Competing interests

We declare we have no competing interests.

## Funding

A.S. acknowledges the support of the NSF CAREER grant no. 1554033 to the Illinois Institute of Technology.

- Received October 12, 2016.
- Accepted March 24, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.