{'fraction_non_alphanumeric': 0.031459270902866, 'fraction_numerical': 0.012553562584784995, 'mean_word_length': 4.8380892520427405, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiband homogenization of metamaterials in real-space: Higher-order nonlocal models and scattering at external surfaces

This work develops a dynamic homogenization approach for metamaterials. It finds an approximate macroscopic homogenized equation with constant coefficients posed in space and time; however, the resulting homogenized equation is higher order in space and time. The homogenized equation can be used to solve initial-boundary-value problems posed on arbitrary non-periodic macroscale geometries with macroscopic heterogeneity, such as bodies composed of several different metamaterials or with external boundaries. First, considering a single band, the higher-order space derivatives lead to additional continuity conditions at the boundary between a homogeneous material and a metamaterial. These provide predictions of wave scattering in 1-d and 2-d that match well with the exact fine-scale solution; compared to alternative approaches, they provide a single equation that is valid over a broad range of frequencies, are easy to apply, and are much faster to compute. Next, the setting of two bands with a bandgap is considered. The homogenized equation has also higher-order time derivatives. Notably, the homogenized model provides a single equation that is valid over both bands and the bandgap. The continuity conditions are applied to wave scattering at a boundary, and show good agreement with the exact fine-scale solution. Using that the order of the highest time derivative is proportional to the number of bands considered, a nonlocal-in-time structure is conjectured for the homogenized equation in the limit of infinite bands. This suggests that homogenizing over finer length and time scales -- with the temporal homogenization being carried out through the consideration of higher bands in the dispersion relation -- is a mechanism for the emergence of macroscopic spatial and temporal nonlocality, with the extent of temporal nonlocality being related to the number of bands considered.


Introduction
Metamaterials have the potential to display unusual properties, and have been the focus of much attention in mechanics and electromagnetism, e.g. [ GRAS19] and numerous others. A primary focus of much of this effort is to find the dispersion relations in the context of infinite periodic systems, and then tailor or optimize these relations to obtain desirable properties.
In this paper, we use the dispersion relation to develop homogenized models, in the context of linear periodic metamaterials. The homogenized models that we develop have two key features: (1) they are posed in real space and time, i.e., not in frequency space; and (2) they are posed in terms of a single equation that is valid over a large frequency range and across multiple bands. These features can enable the application of the models to problems posed on complex geometries and with complex time-dependent behavior.
Our approach is based on developing approximations to the dispersion relation and then inverting it to obtain the real-space model. This requires us to balance between two competing objectives: (1) an accurate representation of the dispersion relation, and (2) tractability of the representation to inversion to real space to obtain standard differential operators. For a single band, we find that a class of models that go back to Toupin [Tou62] provides a good balance. In brief, this class of models can be considered as approximating dispersion curves as rational functions; the inversion to real-space then provides standard but higher-order differential operators. Specifically, the homogenized equation includes strain gradients as well as spatial gradients of the inertial terms. When we extend this to higher bands, we find the appearance of higher-order time derivatives as well.
To summarize our overall procedure, we: (1) start with a given dispersion relation that is derived under assumptions of periodicity; (2) use a rational function approximation of the dispersion relation to develop a homogenized equation in space and time with constant coefficients; and (3) then assume that the homogenized equation is still valid even when the coefficients vary in space and the body is finite. Roughly, our procedure assumes implicitly a separation of scales between the periodic microstructure and the macroscopic heterogeneity.
To test and demonstrate the homogenized model, we turn to an important topic of current interest in the community: to predict the scattering of waves at boundaries between different metamaterials or boundaries between a metamaterial and a homogeneous material, e.g. [SW17,Wil20b,Wil20a,CGMP19]. We use problems extracted from this body of work as examples on which to test our approach, and discuss this in detail below. Our approach, in brief, is to identify the additional nonstandard continuity conditions that arise at a surface of discontinuity due to the higher-order derivatives in the homogenized equations. These higher-order continuity conditions then provide precisely the information required to uniquely solve for the wave scattering at the boundary. A schematic view of this approach is shown in Figure 1.

Prior Work on Scattering at Boundaries between Metamaterials and Homogeneous Materials.
A simplistic approach to model the scattering at a boundary might start from the usual assumption of infinite periodicity to obtain the homogenized properties of the materials on either side, and then use these properties in the nonperiodic setting with a boundary. However, as pointed out in [SW17], such an approach would consider only the propagating waves and miss the contribution from the evanescent waves that are present in the nonperiodic setting. They showed that if such homogenized models are applied to this problem, the evanescent modes cannot be accounted for, and the classical continuity conditions -continuity of displacement and traction -must be satisfied with propagating waves only; this leads to a violation of the energy flux conservation. [SW17] discuss an approach to deal with this problem. Their approach is to minimize the errors in displacement continuity and traction continuity, subject to the conservation of energy, by considering only the propagating modes. While simple and easy to use, it does not consider the microstructural information in formulating this criterion, i.e., it does not account for the details of the evanescent modes. Further, it can be difficult to apply this in a more general setting, for instance in a setting with complex geometry and time-dependence where a decomposition into propagating modes is difficult.
A different "bottom-up" approach to this problem, and analogous problems, has been proposed by other workers, e.g., [MMP18, CG20, CGMP19, PMM21, MM17, MM18, MPM19, GMOI19, MG18, MM16]. These approaches rigorously account for the lack of periodicity, e.g., by introducing boundary-corrector functions or using matched asymptotic expansions in 2-scale methods. Such methods typically require the coupled solution of the cell problem -solving for the fast variable in a periodic unit cell -along with the effective PDE describing the wave transmission, or the computation of the boundary-corrector functions. These requirements make the approaches challenging and computationally expensive, and consequently difficult to apply in a general setting. Further, these methods typically are restricted to a relatively limited frequency band, making them difficult to apply to general problems with complex time-dependent loading.
The Proposed Approach.The approach that we propose aims to sit between the approach of [SW17] and the rigorous homogenization approaches. Specifically, we aim for our homogenized models to be (1) applicable to a broad range of frequencies, (2) applicable to general geometries and time-dependent loadings, (3) computationally efficient and easy to compute with.
In the specific context of the insights provided by [SW17] on the importance of evanescent waves, our homogenized models are able to capture the evanescent modes and the bandgap. Specifically, in higher dimensions, there are typically 2 propagating scattered waves: a reflected wave and a transmitted wave. While the classical continuity conditions -i.e., that traction and displacement are continuous -enable us to approximately find these propagating waves, this approach neglects completely the effects of the evanescent modes. Further, [SW17] show that there exist important regimes wherein there are more than 2 propagating waves, and the 2 classical continuity conditions are insufficient to approximately find the propagating waves.
Our approach provides higher-order continuity conditions ( Fig. 1) at the interface that augment the standard continuity conditions. In general, these additional continuity conditions enable us to better account for the evanescent waves, thereby providing more accurate solutions. In the specific situations studied by [SW17], the additional conditions enable us to uniquely find the propagating scattered waves.
Emergence of Temporal Nonlocality with Multiple Bands. Existing dynamic homogenization techniques can provide a homogenized equation for a specific band of the dispersion curve by expanding about the frequency in that band, e.g. [GMOI19,CKP10,HMC16]. However, this typically does not provide a single equation that covers several bands. In this work, we propose a single homogenized equation that is applicable over a range of frequencies that covers multiple bands and the bandgaps. An important feature of this homogenized equation is that it has higher derivatives in time, suggesting the emergence of nonlocality in time due to homogenization as we increase the number of bands. Other researchers have developed the idea of nonlocality in time as an outcome of homogenization: foundational early work in this area includes the approach by Willis [Wil80,Wil81,Wil97], and the mathematical work of Tartar [Tar89,Tar91] and others following him [Ant93]. However, these and other prior works does not seem to have explored the path of accounting for higher-order dispersion bands as a mechanism for temporal nonlocality.
Organization. In Section 2, we formulate the exact fine-scale model and discuss the specific examples that we will use to compare the exact and homogenized models. In Sections 3 and 4, we describe the development of the homogenized single-band model in 1-d and 2-d respectively, and compare its predictions with those of the exact model. In Section 5, we describe the development of the homogenized model for two bands with a bandgap, and compare its predictions with those of the exact model. In Section 6, we apply our approach to the setting in which there can be more than 2 propagating waves, wherein the classical continuity conditions are insufficient to provide a unique solution. In Section 7, we discuss nonlocal equations as possible homogenized models that account for an infinite number of bands.

2.A. Notation
We will refer to the model that considers the detailed microstructure as the fine-scale model, and its solution as the fine-scale solution. Analogously, we will refer to the homogenized effective dynamical PDE as the homogenized model, and its solution as the homogenized solution.
The primary variable, i.e. the displacement in the mechanical setting, will be denoted u(x, t), where t denotes time and x is the spatial coordinate. In 1-d, x ≡ x; in 2-d, x ≡ (x 1 , x 2 ) in Cartesian coordinates. Our 2-d problems will be in the antiplane setting, so u is always a scalar.
For conciseness, we use the following representation for derivatives: and so on. We use the summation convention, e.g., ∂ ii ≡ ∂ 2 1 + ∂ 2 2 . We consider bodies wherein the material properties -either in the fine-scale model or in the homogenized model -can be discontinuous across boundaries. Consequently, various quantities can be discontinuous at these boundaries, and the jump in a quantity across a surface of discontinuity is denoted · . The surfaces of discontinuity are assumed to be fixed in space.
The Bloch wavevector in the fine-scale models is denoted K = (K 1 , K 2 ), and the wavevector in the homogenized setting is denoted k = (k 1 , k 2 ). The frequency is denoted ω, and will be normalized by the midgap frequency ω 0 .

2.B. Fine-scale Governing Equations
We will consider the linear wave equation in heterogeneous media. Denoting the domain of the body by Ω, we have: where the fine-scale material is assumed to be isotropic for antiplane shear, and µ(x) and E(x) denote moduli, and ρ is the density. Let S denote a stationary surface of discontinuity in Ω. The continuity conditions on S are given by: wheren(x) is the unit normal to S. The continuity conditions correspond to the classical conditions of displacement continuity and traction continuity respectively.

2.C. Problem Geometries and Boundary
Conditions a. One-dimensional Problems.In 1-d, we consider two problems, as shown in Figure 2. In both problems, a given traveling wave originates at x → −∞ in the homogeneous material and is incident at the boundary (a) Single boundary geometry in 1-d, with a semi-infinite laminate metamaterial adjoining a semi-infinite homogeneous material.
(b) Double boundary geometry in 1-d, with a laminate metamaterial of finite extent L (taken to be 10 unit cell widths) between two identical semi-infinite homogeneous materials. between the homogeneous material and the metamaterial. At x → ∞, a Bloch wave form is used in the semi-infinite metamaterial, and a traveling wave is used in the semi-infinite homogeneous material in the single-and double-boundary problems respectively. This formulation, and the solution, can be found in many works in the literature, e.g., Section 2 in [Wil16] and [CG20]. The primary quantities of interest are the magnitudes and energies of the reflected and transmitted waves at steady state; we aim to compare the homogenized model predictions of these quantities with the fine-scale model predictions for various frequencies. b. Two-dimensional Problem. In 2-d, we consider the geometry shown shown in Figure 3. We consider plane traveling waves originating at infinity in the semi-infinite homogeneous material, and oriented at an angle θ to the tangent of boundary of the metamaterial. The plane wave scatters on the boundary giving rise to transmitted and reflected waves. Because the boundary is infinite and planar, the transmitted and reflected waves are also plane waves and we apply a Bloch form at infinity in the semi-infinite metamaterial. This formulation, and the solution, follows [Wil16,SW17].
The primary quantities of interest are the magnitudes and energies of the transmitted and reflected waves at steady state, and we compare the homogenized model predictions and fine-scale model predictions for various frequencies and various angles of incidence θ.

2.D. Unit Cell Material Properties
Throughout this paper, we consider a periodic bilayer laminate metamaterial composed of two materials a and b; the geometry and notation is shown in Figure 4. Table 1 lists the numerical values of the material properties chosen for all the examples (both 1-d and 2-d) involving a single band. Table 2 lists the numerical values of the material properties chosen for the examples involving 2 bands in Section 5; the key difference is a bigger contrast to have a larger bandgap. The bandgap frequencies used for normalization for the properties from Table 1 and Table 2 are ω 0 = 3.448 685 rad s −1 and ω 0 = 8.646 75 rad s −1 respectively.
The homogeneous material has modulus and density denoted E h and ρ h respectively. We use a numerical value of ρ h = 1 throughout. For the 1-d single boundary problem, we use E h = 0.25 when we examine 1-band approximations and E h = 8 when we examine 2-band approximations; for the 1-d double boundary problem and the 2-d problem, we use E h = 2. These values are chosen to provide numerical results that can be easily compared when plotted.

2.E. Dispersion Relations for the Fine-scale Model
The dispersion relation for a periodic laminate in 1-d is classical and given in, e.g., [Wil16,Lek94]: Figure 5a plots this dispersion relation using the properties in Table 1.
The dispersion relation ω(K) for a periodic laminate in 2-d is given in, e.g., [Wil16,CG20]. It is plotted in Figure 5b as a function of the wavevector components K = (K 1 , K 2 ).

Single Band Homogenized Model in One Dimension
In this section, we examine the simple setting of a 1-d model with a single band. We begin by identifying an appropriate approximation to the dispersion relation; the approximation captures the primary feature   that waves above a certain frequency are not supported, and at the same time is amenable to be inverted to find the homogenized real-space-time dynamical equation (Section 3.A). We then identify the continuity conditions that correspond to the homogenized dynamical equation in a macroscopically-heterogeneous medium (Section 3.B). Finally, we apply these continuity conditions to study scattering at macroscopic boundaries in specific examples (Section 3.C).

3.A. Approximating the Dispersion Relation
We begin by comparing the dispersion relations of classical linear elasticity, strain gradient elasticity, and "microinertia" models in the context of homogeneous materials. Such models were introduced by Toupin

Dynamical equation Dispersion
Classical elasticity: The dispersion relations have similar behavior at k 2 → 0. The important distinction is the behavior at k 2 → ∞. In this limit, classical elasticity and strain gradient models have the undesirable feature that ω 2 → ±∞. In contrast, microinertia models have a finite limiting value of ω if we have equal orders of spatial derivatives on both sides of the momentum balance. We notice that strain gradient models lead to ω 2 (k) that are polynomials, and microinertia models lead to ω 2 (k) that are rational functions. Rational function approximations have the important advantage that they are more flexible, e.g. they can approximate curves with finite limits at k 2 → ∞, while also being able to easily invert to real space and time to give dynamical PDE posed in terms of standard differential operators.
Given this advantage of rational function approximations, our strategy begins by approximating a given exact dispersion relation such as (2.5) using a rational function approximation. Our approximation has the form: Some important features of this approximation are: 1. All coefficients in the rational function approximation are taken to be non-negative. This avoids singularities and the frequencies are real for all k. 2. The exponent of the highest-order terms in the numerator and denominator are balanced to ensure a finite limit at k → ∞. 3. The dispersion relation is an even function of k, and hence we will use only even powers of k in the rational function approximation. 4. While our rational function expression only approximately matches the exact given dispersion relation, we constrain the approximation to ensure that we exactly capture the long-wavelength and static behavior. Specifically, we require that our approximation satisfies where c is the long-wavelength wave velocity. This is ensured by taking D 0 as the mean density, and N 1 as the effective modulus obtained from static homogenization which is just the harmonic mean of the moduli. 5. We will denote the different rational functions by RF N N , where N refers to the exponent of the highest-order term in the numerator or denominator.

3.B. Homogenized Dynamical Equation and Boundary Continuity Conditions
We consider (3.4) with N = 6 for concreteness: Given an exact dispersion relation obtained from dynamic homogenization of a fine-scale metamaterial model, we select the coefficients N 0 , N 1 , N 2 , D 0 , D 1 , D 2 that best approximate the given exact relation.
We rewrite (3.6) as the dynamical equation in Fourier space: whereũ(k, ω) is the Fourier transform of u(x, t) in both space and time. We then use the inverse Fourier transform to obtain the (approximate) homogenized dynamical equation in real space and time: We emphasize that this Fourier inversion assumes implicitly that the coefficients are all constant in space.
We next construct the corresponding Lagrangian: and the action: where Ω is the body and [0, t] is the time interval of interest. The variation of S RF 66 with respect to u(x, t) recovers (3.7) as the Euler-Lagrange equation.
Here, we make the central assumption of our method. The dynamical equation (3.7) was homogeneous in space, i.e. the coefficients N 0 , N 1 , N 2 , D 0 , D 1 , D 2 , D 3 were constant. We now take (3.8) and (3.9) as the starting points of our homogenized model, and assume that they hold even if the coefficients are functions of x. Taking the variation of S RF 66 with respect to u, we obtain the dynamical equation: and the corresponding continuity conditions at discontinuities: The macroscopic homogenized model consists of (3.10) and (3.11)-(3.16). We can identify a conserved energy for the homogenized model by multiplying (3.7) by ∂ t u and integrating over Ω: The energy is clearly positive-definite and its time derivative can be readily computed to be 0 if there is no energy input from the boundaries of Ω.

3.C. Numerical Comparisons between the Homogenized and Fine-scale Models
We consider the metamaterial with the numerical values of the properties as chosen in Section 2.D (Table 1), and the dispersion relation in (2.5). The exact dispersion relation and several rational function approximations are shown in Figure 5a, and the values of the coefficients of the rational function approximations are given in Table 3. The rational function approximations are computed using a least-squares fit.
For both the single-and double-boundary problem, we consider both cases of lamina (a) and lamina (b) being adjacent to the homogeneous material. While the solution to the corresponding fine-scale models differs slightly, the homogenized model cannot -even in principle -capture these minor differences because we do not account for the details of the interface structure.   Table 3. Coefficients in the rational function approximations obtained using least squares fits to the exact dispersion relation (2.5). Figure 6 shows plots of the rational functions. N 1 = 1.190 476 19 is chosen to match the classical linear elastic limit as k → 0 for all of the approximations. For the homogeneous medium, N 1 = E h and D 0 = ρ h , and all other coefficients are zero.
a. Single-Boundary Problem.We consider the single-boundary problem described in Section 2.C ( Figure  2a). The numerical values of the homogeneous material properties are given in Section 2.D. We solve the homogenized model -(3.10) and (3.11)-(3.16) -by noticing that (3.10) has piecewise constant coefficients and hence the solution is simply a piecewise superposition of exponentials. The continuity conditions are then used to obtain the coefficients in the superposition solution. The solution method is easy and completely standard; the difference is that we need to solve for more coefficients than in the classical case, and we also have more continuity conditions to use. The comparison between the predictions of the exact fine-scale model and the homogenized model is shown in Figure 7. Specifically, we show the predictions of normalized energy transmission across the boundary as a function of frequency. The error in the homogenized model is remarkably small given the significant simplifying assumptions; the error is about 5% at about 80% of the bandgap frequency.   Figures 8, 9 and 10. Figure 8 shows the transmitted energy as a function of frequency. We notice the classical dips in the transmitted energy that correspond to the destructive interference between the transmitted waves and the multiple reflected waves. The homogenized model is able to capture the frequencies at which these dips occur and the overall qualitative structure of the curve well using the RF 66 approximation which significantly overestimates the frequency near the bandgap (Figure 6), with an error of about 20% at about 80% of the bandgap frequency.   Figure 2b. We highlight that these predictions are based on the RF 66 rational function approximation that has a significant overestimate of the bandgap frequency (Figure 6), yet it is able to provide good qualitative and fair quantitative predictions. Figure 9 shows the displacement field at a frequency about halfway to the bandgap, using the simplest RF 22 approximation and comparing with the exact fine-scale model. We find a very good agreement for the real and imaginary parts (i.e., the homogenized model captures the phase information correctly).
If we use a higher-order approximation, e.g. RF 66 , we are able to capture the bandgap much more accurately (Figure 6). Figure 10 shows an example of a wave propagating with a frequency ω/ω 0 = 1.025, i.e., it lies in the bandgap of the metamaterial. We notice that the evanescent modes at the boundary are captured by the homogenized model, though the decay is somewhat faster than that of the exact fine-scale model.

Single Band Homogenized Model in Two Dimensions
We next consider the 2-d setting, and follow the overall strategy of Section 3. For numerical comparisons, we consider the problem described in Figure 3. We begin by approximating the first band of the exact dispersion relation in Figure 5b using a rational function approximation of the type: where N (0) , D (1) are second-order tensors, N (1) , D (2) are fourth-order tensors, and so on. All of these tensors are symmetric and positive definite to avoid singularities in the approximation.
For the numerical calculations presented in this section, we keep only terms up to |K| 2 . A least-squares fit against the exact dispersion relation gives:

4.A. Homogenized Model and Continuity Conditions
Transforming (4.1) to real space and time gives the 2-d homogenized dynamical equation: The corresponding Lagrangian L 2d is: Dropping the assumption of homogeneity of the coefficients and using the principle of least action, we get the following continuity conditions on the boundary, wheren(x) is the unit normal to the boundary.

4.B. Numerical Comparisons between the Homogenized and Fine-scale Models
We use the 2-d problem described in Section 2.C and Figure 3, and with the numerical values for the material properties from Section 2.D and Table 1.
The method of solution of the homogenized equations is analogous to the 1-d case. We first notice that (4.3) is linear with piecewise constant coefficients, and hence the solutions are simply a piecewise superposition of exponential functions. The coefficients in the superposition are obtained by applying the continuity conditions (4.5)-(4.6). This gives the scattered displacement field and scattering coefficients for the homogenized model. Figure 12 compares the reflection coefficients predicted by the exact fine-scale model and the homogenized model as a function of frequency, for several frequencies and over the entire range of incident angles from grazing to normal incidence. There is very good quantitative agreement over all the frequencies and incident angles. Figure 13 plots the displacement field of the fine-scale model and the homogenized model for a specific frequency and θ = 30 • . Visually, the agreement is very good.
Our consideration of this problem is motivated by [SW17] that examined the same problem in a different way. They highlighted that when the metamaterial lamination is normal to the metamaterial boundary, there are a finite number of propagating waves and an infinite number of evanescent waves generated in the metamaterial. Standard dynamic homogenization can account only for the propagating waves because it uses the setting of an infinite periodic medium. If such homogenized models are applied to this problem, the evanescent modes cannot be accounted for, and the classical continuity conditions -continuity of displacement and traction -must be satisfied with propagating waves only; this leads to a violation of the energy flux conservation. In the case of multiple propagating waves in the 2-direction, the classical continuity conditions are not sufficient to determine the transmission and reflection coefficients.
The proposed approach of using higher-order homogenized models gives additional continuity conditions at the metamaterial boundary. Further, the higher-order homogenized model supports evanescent waves. These features address the issues raised in [SW17] and enable good quantitative predictions.

A Single Homogenized Equation to Describe Two Bands with a Bandgap
We now consider the question of modeling multiple bands, and focus on two bands with a bandgap to illustrate the idea in a simple setting. Prior work, e.g., [CKP10], has derived approximations for higher bands, but these approaches are restricted to a single band in the vicinity of an a priori selected frequency. In contrast to such approaches, we obtain a single homogenized equation that is seamlessly valid across the entire frequency range spanning both bands, without having to specify a priori the band or limit the frequency range of interest. Such a model is essential to enable the modeling of complex time-dependent transient loadings, such as shock loading, and waves that are composed of a broad range of frequencies.
To obtain a composite dispersion relation that can represent both bands together, we simply use the product of the individual dispersion relations: The resulting composite dispersion relation is easy to invert to real space and time to obtain the multiband homogenized equation: To simplify notation, we define: To find the Lagrangian L M B [u](t), we first multiply (5.4) by a test function v(x, t) and integrate over the spatial domain Ω and an arbitrary time interval [t 1 , t 2 ]. Using standard integration-by-parts operations, we can simplify to write: where we have ignored boundary terms in space and time, by assuming that the initial conditions are given and hence v and its time derivatives are zero at t = t 1 .
This can be used to write the action: where the Lagrangian is: Extremizing S M B without assuming that the coefficients are constant gives the homogenized equation: and the continuity conditions: In Appendix A, we use L M B and S M B to obtain the conserved energy (A5) and consequently an expression for the energy flux across a surface (A6). The energy flux expression is required to compute the transmitted and reflected energies. [Ost50] shows that the corresponding Hamiltonian (energy) has linear instabilities, referred to as Ostrogradsky's instability. However, our homogenized model is marginally stable in the sense of Lyapunov when the coefficients are chosen in the appropriate range.

5.A. Numerical Comparisons between the Homogenized and Fine-scale Models
We use the 1-d single-boundary problem described in Section 2.C and Figure 2a, with the numerical values for properties listed in Section 2.D (Table 2).
An important aspect of the homogenized model is that the metamaterial is described using higher-order time derivatives. To have a consistent description, we therefore require that the homogeneous material also be modeled using higher-order time derivatives of the same order. We achieve this by using considering the homogeneous material as a bilayer laminate metamaterial with identical layers, and constructing a two-band higher-order model (Fig. 14). Near Kh = 0, the linearity of the first band is captured well by the rational function approximation, however, the linearity of the second band is only captured approximately by the second-order rational function approximation. Even with this coarse approximation, we find good qualitative and quantitative agreement.
We solve the homogenized model by first using that (5.8) is linear with piecewise constant coefficients; hence, the solutions are the standard piecewise superposition of exponentials with unknown coefficients. The coefficients are obtained by applying the continuity conditions (5.9)-(5.12). A comparison of the transmitted energy obtained using the fine scale model and the homogenized model is shown in Figure 15. The band gap frequencies can be identified as the range of frequencies over which there is no transmission, and is captured very well by the homogenized model.

A Homogenized Model for Multiple Propagating Waves at High Frequencies
As we consider higher frequencies, there are multiple propagating waves in the metamaterial along with infinitely many evanescent waves. Figure 16 shows the K 1 h vs. K 2 h curves for different frequencies of the first two bands. We can see that there is a high frequency regime wherein there is a range of values of K 1 h for which we have two values of K 2 h; this implies that there are 2 propagating waves in the metamaterial at those frequencies. Figure 14. Dispersion relations for the metamaterial (Kh > 0) and the homogeneous material (Kh < 0). The periodicity of the artificial uniform laminate that represents the homogeneous material is chosen to be the same as the periodicity of the laminate metamaterial. The dispersion bands for the metamaterial and the artificially-periodic homogeneous material are shown by the black curves (exact) and red curves (rational function approximation).

6.A. Insufficiency of Existing Homogenization Models
For the problem described in Figure 3, where the interface is at right angles to layers in the metamaterial, the displacements fields for the incident, transmitted, and reflected waves are respectively: This solution uses the Bloch ansatz [SW17].ũ n (x 1 ),w m (x 1 ) are the Bloch modes in the metamaterial and the homogeneous medium respectively; K 0 is the wavevector of the incident wave; and the wavevector in the x 2 direction is denoted by K for the reflected field.
At low ω, K 2 is the only real-valued wavevector, and K (n) 2 is imaginary for all n > 1; this implies that there is only one propagating transmitted wave.
At higher ω (Fig. 16), K (n) 2 for n = 1, 2 are real-valued, and K (n) 2 is imaginary for n > 2; this implies that there are now 2 propagating transmitted waves in the metamaterial. In general, we could similarly have 1 or more reflected propagating waves in homogeneous medium, but for simplicity we consider only frequencies wherein there is a single reflected propagating wave.
To solve for the scattered propagating waves, we note that K 1 is conserved across the interface; i.e., the 2 propagating transmitted waves have the same real value of K 1 h, but distinct real values of K 2 . Even if Figure 15. Comparison of the transmitted energies predicted by the exact fine-scale model (black curve) and the homogenized model (red curve) for the approximation shown in Figure 14. The band gap is captured exactly as the rational functions are tuned to match the band gap exactly. we completely ignore the evanescent modes, we must determing 3 scattering coefficients: 2 transmission coefficients (T 1 , T 2 ) for the 2 transmitted propagating waves, and 1 reflection coefficient (R 0 ) for the reflected propagating wave. As discussed in Section 4.B, the classical conditions of displacement and traction continuity at the interface provide only 2 equations, but we have to find 3 coefficients.

6.B. Homogenized Model for Two Bands in Two-dimensions
We combine the ideas from Section 4 and Section 5 to obtain a multiband model for two bands in 2-d. The first and second bands of the exact dispersion relation are approximated using rational functions given bŷ ω a (K) andω b (K) respectively.
where the coefficients N , D, P , and Q are higher-order symmetric tensors. The 2-band dispersion relation can be formed by taking the product of the individual approximations, following the idea of (5.3).
Transforming the dispersion relation to real space and time provides the multiband homogenized dynamical equation in 2-d. Following prior sections, we construct a Lagrangian corresponding to the dynamical equation, and then use the principle of least action to derive the dynamical equation and the corresponding jump conditions: −(A (6) mn ∂ n u) n m = 0 (6.7) mnpq ∂ pq u n n = 0 (6.9) ∂ i u = 0 (6.10) Here,n(x) is the unit normal to the boundary. For conciseness, we have defined: A (1)

6.C. Quantitative Results
We begin by constructing the homogenized approximation of the metamaterial. To enable comparisons, we compare the rational function construction with a polynomial construction of the dispersion relation (Fig.  17). Using the simple rational functions of the form shown in (6.4) and (6.5), the coefficients are obtained from a least-squares fit to the exact dispersion relation: Consider the 2-d problem described in Section 2.C and Figure 3. Choosing a frequency ω/ω 0 = 1.218, at which we see two propagating transmitted waves in the metamaterial, we consider an incident wave with wavevector K 1 h = 2.9. Since K 1 h is conserved across the interface, the metamaterial has two propagating waves with K 1 h = 2.9 but different values for K 2 h. The homogeneous medium has just one propagating reflected wave. The fine-scale energy transmission is computed following Section 5 from [Wil16]; the transmitted energy (E trans /E inc ) from the fine-scale model is 0.7914, evaluated by using as many evanescent modes as needed for the error to be below 0.1%. The multiband homogenized model provides enough continuity conditions (6.7)-(6.10) to solve for the scattering coefficients. The transmitted energy computed from the homogenized model based on the polynomial approximation of the exact dispersion relation is 0.2844, while the model based on the rational functional approximation predicts a transmitted energy of 0.7001.
We draw 2 conclusions: first, the homogenized model based on rational function approximations does well quantitatively given the relative crudeness of the approximation; and, second, the comparison to the polynomial-based approximation suggests that a good approximation of the dispersion relation provides significant gains in accuracy. Based on these conclusions, we conjecture that further refinement of the dispersion approximation, through the use of higher-order models and the corresponding higher-order continuity conditions, will provide further quantitative accuracy.

Nonlocal Limit Models
Section 5 showed the emergence of 4th-order time derivatives when we considered 2 bands. In principle, it is straightforward to consider as many bands as one wishes by continuing the approach in (5.3). For instance, if we wish to consider N bands, we can write each band in terms of a rational function approximationω i (k), and construct the composite dispersion relation as: The resulting dispersion relation will include (even) powers of ω of up to order ω 2N , and the corresponding homogenized model in real time and space will include (even) time-derivatives of up to order 2N .
We recall that the fine-scale model will generically have an infinite number of bands. This raises the question of what limit homogenized model we might expect if we aim to represent all bands in our homogenized model. We could think of each higher-order time derivative as containing information from further out in time; the (heuristic) limit of an infinite order time-derivative might then correspond to a nonlocal operator in time.
Our heuristic reasoning parallels work in peridynamics [Sil00], wherein formal arguments based on Taylor expansions have been used to relate nonlocal integral operators in space to an infinite series of spatial derivatives of increasing order, following the idea of matching dispersion relations [WA05,Day17,BD18,SPGL09,SDP16]. From that perspective, one could think of peridynamics as a special case of the limit models discussed here, in that it is the limit of a single band model.
In [Day17], it was shown formally that the limit of microinertia models -corresponding to a single-band model in the perspective put forward in this paper -is consistent with a nonlocal integral equation of the form: (7.2) where C 1 (x, x ) and C 2 (x, x ) are spatial kernels. We note that the continuity conditions for nonlocal operators are much weaker than in models with higher-order derivatives.
Formally, if we consider the limit of time derivatives of higher order, we can arrive at a limit model that has nonlocality in both space and time is of the form: is a time-kernel, and T 0 is the fixed extent of non-locality in time. The structure of the term that is nonlocal in both space and time follows from (7.2), and the term that is nonlocal in time alone has a Taylor expansion that matches the structure of the terms in (7.1).
Remark 7.1 (Willis-type Nonlocal Homogenized Descriptions). Willis' approach [Wil80,Wil81] provides an alternative strategy to develop a nonlocal limit model; it is based on using Fourier space and averaging over space and time. In contrast, while we also work using Fourier space, we then homogenize over space and time through the mechanism of considering the dispersion bands.

Discussion
We have presented an approach to develop homogenized models for metamaterials that are posed in space and time, and consequently potentially applicable to general settings such as complex geometries and complex time-dependent loading and to a broad range of frequencies. The homogenized models have higher-order derivatives in space and time, leading to nonstandard continuity conditions at the boundaries between metamaterials.
While the predictions of the homogenized models compares well with exact fine-scale models, there are a number of important open questions. The most central question is a justification of the method in a more "bottom-up" fashion. While our approach provides models that are much simpler to use than the rigorous approaches, e.g. [CG20], it is important to understand the connection between our approach and these rigorous approaches; in particular, if it is possible to view our approach as a simplified version of the rigorous approaches under certain assumptions. Similarly, computational homogenization schemes developed for dynamic situations, e.g. [LR17, HO14, PKG13, SKG18] provide alternative, but potentially more expensive, approaches to the problems considered here. These may provide a route to connect the approach in this paper to finer scales using computational approaches; however, it is also important to demonstrate that these methods can effectively capture the evanescent waves at interfaces, which appears to be open at this time.
A second open question is how to perform practical numerical calculations with our model, given the nonclassical higher-order derivatives. In the context of finite element (FE) methods, two options appear feasible. First, second-order rational functions provided fairly good accuracy. They have spatial derivatives of 4th-order, making them amenable to Isogeometric FE methods which provides sufficient continuity for such derivatives [KB19]. Second, mixed FE methods are a possible approach for problems that require very smooth interpolations, though mixed methods can require care to formulate correctly. The implementation of numerical methods will enable application of the approach to numerous problems of current interest in the broad area of metamaterials.
An important avenue for future development is to go beyond problems with a scalar unknown field to consider vector problems that come up in elasticity and electromagnetism. While the overall principles do not change, one must account for the tensorial nature of problems with vector unknowns. In particular, the dispersion relation that must be inverted now requires tensorial quantities to be approximated as functions of the wavevector; a component-wise approximation is perhaps a feasible strategy to attack this problem. Similarly, we have restricted our work in this paper to 1 and 2 spatial dimensions and this must be extended to 3 dimensions; however, in our view, this is a relatively straightforward next step.
We highlight that this work, and the many others cited here, are all in the linear regime. Homogenization of wave scattering in the nonlinear regime is largely unexplored, though we mention a study of wave scattering at an interface between a phase-transforming solid and a linear elastic solid [AK92]. On the other hand, we also note that the linear dispersion behavior is important even in nonlinear problems, for instance to govern the effective dissipation [DB06] or stability [Abe14].
We note that some aspects of our overall approach -namely, the approximation of dispersion relations and subsequent inversion to real space and time -is broadly similar to strategies used in the construction of nonreflecting boundary conditions (NRBC) [BT80,GK90]. In the NRBC literature, such an approach has been used to construct wave equations that propagate waves in a tailored manner.
In the context of linear homogenization, we note other interesting works that obtain homogenized models that have features in common with our homogenized models. Specifically, the elimination of fine-scale degrees of freedom has been found to give rise to effective nonlocality in space by [PSP18,DW96], and also to give rise to very nonstandard dynamical behavior in [MW07,PZ12]. Further, [Wil80,Wil81,Tar91,Ant93,Tar89] find the emergence of nonlocality in time due to homogenization. Our work provides heuristic insight into an alternative mechanism: specifically, it suggests the emergence of temporal nonlocality as an outcome of spatial homogenization through the consideration of higher-order dispersion bands that corresponding to modes with finer spatial structure.
A consequence of our approach is also the perspective that peridynamics can be considered as the limit of a spatially averaged model that consider only a single dispersion band; the physical origin of the nonlocality in space is the coarse-graining of the fine structure. Analogously, a richer model that considers all dispersion bands will provide the full spatio-temporally nonlocal model that is equivalent to averaging in time and space, and can be useful for recent work in nonlocal modeling that includes wave propagation and dynamics [YYSD20,YYSD22].