On band gaps of nonlocal acoustic lattice metamaterials: a robust strain gradient model

We have proposed an “exact” strain gradient (SG) continuum model to properly predict the dispersive characteristics of diatomic lattice metamaterials with local and nonlocal interactions. The key enhancement is proposing a wavelength-dependent Taylor expansion to obtain a satisfactory accuracy when the wavelength gets close to the lattice spacing. Such a wavelength-dependent Taylor expansion is applied to the displacement field of the diatomic lattice, resulting in a novel SG model. For various kinds of diatomic lattices, the dispersion diagrams given by the proposed SG model always agree well with those given by the discrete model throughout the first Brillouin zone, manifesting the robustness of the present model. Based on this SG model, we have conducted the following discussions. (I) Both mass and stiffness ratios affect the band gap structures of diatomic lattice metamaterials, which is very helpful for the design of metamaterials. (II) The increase in the SG order can enhance the model performance if the modified Taylor expansion is adopted. Without doing so, the higher-order continuum model can suffer from a stronger instability issue and does not necessarily have a better accuracy. The proposed SG continuum model with the eighth-order truncation is found to be enough to capture the dispersion behaviors all over the first Brillouin zone. (III) The effects of the nonlocal interactions are analyzed. The nonlocal interactions reduce the workable range of the well-known long-wave approximation, causing more local extrema in the dispersive diagrams. The present model can serve as a satisfactory continuum theory when the wavelength gets close to the lattice spacing, i.e., when the long-wave approximation is no longer valid. For the convenience of band gap designs, we have also provided the design space from which one can easily obtain the proper mass and stiffness ratios corresponding to a requested band gap width.


Introduction
During the last few decades, researchers have been exploiting various potential applications of metamaterials [1][2] . Metamaterials refer to some composite materials that have artificially designed structures and exhibit extraordinary physical properties that the source materials do not have. For example, acoustic metamaterials can have negative effective mass or negative effective modulus, leading to forbidden band gaps, negative refractive index, stealth, etc. [3][4] It also has a wide range of applications in vibration reduction, noise isolation, acoustic control, waveguide, and so on. Thus, the developments of acoustic metamaterials have greatly expanded the application range of acoustic engineering. Intensive research efforts have been devoted to the analyses of the propagation and dispersion characteristics of elastic waves in acoustic metamaterials [5][6][7][8] with negative effective properties [9][10][11] .
The so-called band gap of metamaterials refers to a certain frequency range in which the wave propagation is forbidden. For metamaterials with local resonator, band gaps generate near the resonant frequency of the internal structure, which is hardly found in materials with simple monatomic lattice structures [12] . Plenty of metamaterials with such a characteristic have been designed for vibration isolation and sound insulation. Many mechanics-based models have been proposed to analyze the band gap structures of structured materials. Huang and Sun [13] studied the band gap structure of locally resonant acoustic metamaterials, and illustrated the effects of mass and spring parameters on band gaps. The band gaps of acoustic metamaterials with multiresonators were investigated by Zhou et al. [14] . Tan et al. [15] proposed a kind of dual-resonator microstructure design to obtain the optimal effective negative mass in acoustic metamaterials. Chen et al. [16] studied how to control band gaps in an active elastic metamaterial with negative capacitance piezoelectric shunting. Liu and Reina [17] investigated the effect of hierarchical microstructures on the band gap structure of periodic lattice systems with local resonators. An et al. [18] discussed the band gaps and vibration properties of two-dimensional (2D) disordered lattice acoustic metamaterials. The modified metamaterial system with local resonators coupled was studied by Hu et al. [19] and Zhao et al. [20] . With the developments in three-dimensional (3D) printing technology, acoustic metamaterials with various microstructures can be readily made.
As a benchmark problem, the dispersion of one-dimensional (1D) periodic mass-spring system has been investigated not only as the discrete lattice model, but also as the continuum model with microstructures considered. Mindlin [21] proposed the strain gradient (SG) elastic theories to account for microstructural effects, which has influenced the solid mechanics community for several decades. To date, the SG model has been successfully used to capture size effects in mechanical properties, when the specimen size gets close to the material's characteristic length [22][23][24] . Parallelly, it has also been applied to dynamic problems when the wave-length and the material's characteristic length are on the same order of magnitude. Polyzos and Fotiadis [25] checked two simple 1D models to explicitly confirm the effectiveness of Mindlin type SG continua. For 1D periodic lattice, the material's characteristic length mentioned above is explicitly related to the atomic spacing. As a follow-up discussion, De Domenico et al. [12] adopted a new truncating strategy for Taylor expansions of the displacement field, and found that the continuum model works well at the long-wave extreme and suffers from an instable issue when the wave length gets close to the atomic spacing. A proper continuum model is helpful in efficiently capturing dispersive characteristics instead of spending too substantial computational costs.
Continuum models have also been built for metamaterials with mass-in-mass resonators. Huang et al. [11] proposed the multi-displacement (MD) continuum model for 1D mass-in-mass metamaterials. This model can accurately predict the dispersion relation in the long-wave domain. Zhou et al. [26] extended the MD continuum modeling to 2D metamaterial plate. Metrikine and Askes [27] built a nonlocal continuum model for the same purpose. As a more complex option, the nonlocal gradient continuum models were proposed [28][29] , which are unconditionally stable but fail to provide proper dispersive characteristics in the whole Brillouin zone. Dispersions of the 1D periodic lattice to be studied here have been a benchmark problem and attracted extensive research interest. Even though it can be easily solved by the discrete model, a robust continuum model capable of giving correct results throughout the whole first Brillouin zone is still lacking [12,30] .
It is also meaningful to investigate the effects of nonlocal springs installed in lattices. Challamel et al. [31] studied the dynamic behaviors of generalized axial lattices with direct and indirect interactions by considering extra higher-order boundary conditions. Zhang et al. [32] discussed the vibration of the multiply connected bar-chain system with direct and indirect interactions and emphasized the effects of long-range interactions on mechanical vibration. Ghavanloo and Fazelzadeh [33] investigated the effect of long-range interactions on the wave propagation in 1D acoustic metamaterials.
In this study, we attempt to construct an exact SG continuum model to predict dispersive characteristics of diatomic lattice metamaterials with local and nonlocal interactions throughout the first Brillouin zone. This paper is structured as follows. In Section 2, we employ a discrete model to provide dispersive diagrams of mass-in-mass metamaterials. In Section 3, by modifying the Taylor expansion of displacement, we propose a wavelength-dependent SG continuum model to capture dispersion characteristics of 1D diatomic lattices. By combining the MD model presented by Huang et al. [11] , such a model is extended to mass-in-mass lattices. In Section 4, the proposed model is validated by comparing with other continua and discrete results. In Sections 5 and 6, the effects of stiffness ratio, mass ratio, and SG orders are respectively investigated. In Section 7, the influence of nonlocal interactions on 1D diatomic lattices is studied. The conclusions are given in Section 8.

Discrete model of 1D diatomic lattice metamaterials
We consider an infinitely long diatomic lattice chain of 1D acoustic metamaterials, as shown in Fig. 1(a), where multiple particles and springs are placed periodically. The second atom composed of mass M 1 and mass M 2 as well as the first atom composed of mass m 1 and mass m 2 is arranged adjacently in space. We consider one-and two-neighbor interactions between different atoms, which are represented by spring stiffnesses K 1 , K 1 , and K 3 , K 3 , respectively. The stiffness coefficients of springs that connect the inner and outer mass particles of the first atom and the second atom are respectively K 2 and K 2 . u (j) 1 and u (j) 2 represent the displacements of the outer and inner mass particles in the atom located at j, respectively. In addition, we assume that the springs here are all massless. The distance between each two neighboring identical atoms is a.
For particles with different masses in the nth cell which is enclosed by the red dash rectangle in Fig. 1(a), the corresponding equations of motion (EOMs) are written as [30] One-neighbor interactions (stiffness K 1 ) One-neighbor interactions (stiffness K ' ) Two-neighbor interactions (stiffness ).
The solutions to the above equations are where A 1 , A 2 , B 1 , and B 2 are the amplitudes, k is the wave number, and ω is the angular frequency. Substituting Eq. (5) into Eqs. (1)-(4) yields To ensure that A 1 , A 2 , B 1 , and B 2 have non-zero solutions, the coefficient determinant where the dimensionless parameters are defined as follows: (10) is the dispersion equation of the corresponding diatomic lattice chain.
Particularly, letting K 3 = K 3 = 0 in Eq. (10) leads to the dispersion relation for the diatomic lattice chain only with local interactions. From Fig. 2(a), we find that there are four wave branches in such a diatomic lattice simply with one-neighbor interactions: the fourth, third, optical, and acoustic branches, suggesting a substantial difference from the single-atom counterpart. Furthermore, three band gaps appearing in the dispersive diagram allow a more flexible design of acoustic diatomic lattices.  For each wave branch, it is interesting to discuss the relative motions of particles in a periodic cell. When K 3 = K 3 = 0, Eqs. (6)-(9) are significantly simplified and serve as those for the lattice with only one-neighbor interactions. For such a case, the amplitude ratios of particle vibrations can be determined as where the two expressions for B1 A1 , i.e., Eqs. (11) and (12), can be proven to be identical. Thus, we only employ Eq. (11) for the following discussion.
We take the case with α = 4, β = 2, α 0 = 2, and β 0 = 0.5 to demonstrate relative motions of particles. Based on Eqs. (11), (13), and (14), we plot the vibration amplitude ratios in Fig. 3. Figure 3(a) gives the result for the acoustic branch in Fig. 2(a). At the extreme of ka 2π → 0, all amplitude ratios get close to 1, indicating that for an infinite wavelength, all particles in a unit cell move together like a rigid unit. Throughout the first Brillouin zone, A2 A1 0 suggests that particles M 2 and M 1 always keep the same vibration direction. While particles m 1 and m 2 can move in the direction opposite to M 1 when the reduced wave number falls in some range around 0.5. In Figs. 3(b)-3(d), the most surprising feature may be that the amplitude ratios do not get close to 1 even when the wave number tends to zero. Furthermore, ratios can possess considerable negative values over most of the first Brillouin zone. Keeping in mind that a negative vibration amplitude ratio means opposite motions of two related particles, the continuity of displacement field in each periodic cell is naturally violated. This fact explains why the MD method employed in Section 3 is necessary for diatomic lattices.   The above formulation reduces to that of the diatomic lattice without internal particle by letting M 2 = m 2 = 0, K 2 = K 2 = 0, K 1 = K 1 , and K 3 = K 3 = 0. Now, Eqs. (6)- (9) give for acoustic mode, for optic mode.
When setting the mass ratio M1 m1 = 2, we obtain the dispersion and amplitude ratio diagrams in Fig. 4, which totally agree with those presented by Born and Huang [34] . Additionally, when ka 2π = 0.5, both the denominator and numerator in Eq. (15) tend to zero; but according to the L'Hospital's rule, we confirm that ka 2π = 0.5 is not a singular point of B1 A1 . Optic mode (Born and Huang [34] ) Optic mode (discrete model) Acoustic mode (Born and Huang [34] ) Acoustic mode (discrete model) Optic mode (Born and Huang [34] ) Optic mode (discrete model) Acoustic mode (Born and Huang [34] ) Acoustic mode (discrete model) 3 Robust SG continuum model with a wavelength-dependent Taylor expansion We develop the exact SG continuum model by taking the diatomic chain only with oneneighbor interactions shown in Fig. 1(a) as an example. A detailed discussion on the SG continuum model with nonlocal interactions is provided in Section 7. The chosen unit cell contains two whole atoms, two one-neighbor springs, and two internal springs, as enclosed by the dash line rectangle in Fig. 1(b). In comparison with existing studies [30,33,35] , we focus on a modified Taylor expansion of displacement to obtain the wavelength-dependent SG continuum model.
For the four-mass particles of the unit cell in Fig. 1 where the superscript "(i)" means "internal mass", i.e., m 2 and M 2 . Then, the displacements of the (2n ± 2)th and (2n − 1)th mass particles can be expressed as To enhance the accuracy of SG continuum models, we update the Taylor expansion of displacements corresponding to the nearest and second nearest mass points. Therefore, Eq. (17) can be rewritten as where the value of d I is related to the wavelength and not necessarily equal to atomic spacing a. u 1 , are respectively the first-, second-, third-, · · · , nth-order derivative of the u 1 field, and the same notation rule applies to u 2 . For 1D problems, u 1 = ∂u1 ∂x can be taken as the normal strain, and thus u 1 and u 1 correspond to the first-and second-order derivatives of strain, or equivalently SGs. When the wavelength λ is much larger than the microstructure spacing a, we adopt d I = a, as most continuum theories have been done under the so-called long-wave approximation. Nevertheless, our present aim is to study the first Brillouin zone, and the wavelength can get close to the microstructural size. Keeping in mind the periodicity of wave motion, the displacement at the position x ± λ is equal to that at x, i.e., From Fig. 5, if we attempt to obtain the displacement at the position x ± a by the Taylor expansion with respect to x, to guarantee the accuracy, we need to judge which one is the nearest to x ± a between x and x ± λ, and then set d I as the distance from the nearest to x ± a. Due to the periodicity of motion, displacements at positions of x + lλ (l = 0, ±1, ±2, ±3, · · · ) are equal to that at x. The Taylor expansion should be carried out according to the nearest location among all x + lλ.
Based on the above considerations of motion periodicity, we get (3) For the selected unit cell 1 in Fig. 1(b), the potential energy density and kinetic energy density are, respectively, where Aa is regarded as the volume of the unit cell containing four complete mass particles, which does not depend on the choice of springs. The same is true in the following discussion. Inserting Eqs. (16)-(18) into Eqs. (21) and (22) and truncating the achieved expressions after By applying the Hamilton variation principle, i.e., we can get the following EOMs: M 2ü m 2ü Keeping in mind that u 1 (x, t), u 2 (x, t) are continuous displacements, Eqs. (26)-(29) provide their evolutions in terms of x and t. In this sense, it can be taken as a continuum model for the studied diatomic lattice, just as done by Mindlin [21] , de Domenico et al. [12] , and Polyzos and Fotiadis [25] .
The solutions to the above equations are Substituting Eq. (30) into Eqs. (26)-(29) yields the dispersive equations of the new SG con-tinuum model Adopting the dimensionless parameters mentioned and considering the existence of nontrivial solutions in Eqs. (31)- (34) yield where Vibration amplitude ratios can be given in the form of where η + 1 =

Validation of the present model
The SG model developed above is validated and used to study the effects of various factors in diatomic lattices. In brief, all results given by the present continuum model agree well with those given by the discrete model.
Compared with models in the literature, the present model is expected to achieve a higher predicting accuracy, solely due to the adjustment of the Taylor expansion of displacement. For the demonstrating purpose, we set α = 4, β = 2, α 0 = 2, and β 0 = 0.5. Figure 6 shows the results provided by the MD method and the theory presented by Huang et al. [11] as well as the present model. The sixth-order and tenth-order truncations of the SG continuum model given by Zhou et al. [30] are called ZWLT6 and ZWLT10, respectively. The following trends can be found.
(I) Both the discrete model and continua can produce four dispersive branches; when the reduced wave number gets close to zero, or equivalently at the long-wave limit, all continua agree with the discrete model, which is deemed to be exact in this study. Throughout the first Brillouin zone, these models predict three band gaps, while different continua possess different accuracies.
(II) The MD model can capture the dispersion characteristics of the discrete model only in the partial wavenumber range. When ka 2π < 0.25, the three dispersive branches from bottom to top coincide reasonably well with the discrete results. However, with the increase in the wavenumber, there is an obvious deviation between them. Particularly, the fourth branch given by MD suggests an extraordinary instability having not been observed.
(III) ZWLT6 and ZWLT10 given by Zhou et al. [30] mainly manifest the role of SG orders. On one hand, with a higher SG order, the ZWLT SG model has a broader workable range out of the first Brillouin zone. On the other hand, the model with a higher order suffers from a more severe instability. That is, once the wavenumber exceeds the above workable range, the dispersion curves deviate abruptly away from the discrete counterparts. An unacceptable fact caused is the failure in predicting band gaps. For example, both ZWLT6 and ZWLT10 predict a band gap width between the optic and third branches which is far narrower than the discrete counterpart and gets close to zero. Polyzos and Fotiadis [25] and De Domenico et al. [12] also realized this instability issue of SG continua. Finally, at certain wavenumber values, the acoustic branches of ZWLT6 and ZWLT10 intersect the coordinate axis, which means that they have predicted an imaginary frequency since then. The same phenomenon also appears when predicting the dispersion relationship of 1D monoatomic chains [12,27] .
(IV) Considering the periodicity of the displacement function, we modify the Taylor expansion, note the relationship between d I and wavelength λ, and propose the present model. The figure demonstrates that the present model can successfully capture the dispersion relation of the original discrete model in the entire first Brillouin zone ka 2π ∈ (0, 1) , and the limitation of long wave approximation is eliminated. This illustrates that our method is effective. It is worth noting that we have not emphasized the order too much. We only keep the result to the eighth power of d I , and we will further explain the effect of order in the following section.  Even though both the discrete model and the above SG model can give correct dispersion diagrams, developing continuum models is still necessary. After all, not all lattices can be considered to be infinite, and strictly obey the periodic boundary conditions. When we face a finite lattice composed of an extremely large number of masses with boundary conditions particularly prescribed, the discrete model quickly becomes too time-consuming, and the continuum model becomes a better option.

Effects of mass ratio
: 1D diatomic chain with local interactions as an example As already shown above, there are three band gaps in the dispersion diagram of the 1D diatomic lattice. We study how the band gap behavior is modulated by changing the mass ratio α = M2 M1 = m2 m1 and the spring constant ratio β = K2 K1 = , as shown in Figs. 7 and 8, respectively. Furthermore, variations of both the band gap widths and the lower band bounds are presented in Fig. 9. For all examples, the SG order is set as eight, the reason for which will be explained in the next section. Figure 7 reveals the effect of mass ratio on dispersion. With the increase in the mass ratio, the four branches tend to be positioned at a lower frequency. Relatively, the influence of mass ratio on the fourth branch is the weakest. And the curve of acoustic branch becomes smoother.
When the values of M2 M1 and m2 m1 increase from 1 to 4, the band gap range changes greatly; while for those increase from 8 to 12, the band gap range changes slightly. From Figs. 7(a)-7(d), the band gap between the optical branch and the acoustic branch becomes narrower, while that between other branches changes slightly.  When designing band gap metamaterials, the frequency level and gap width are two key parameters. As a preparation, we define the lower bound and width of each band gap in Table 1. For example, for the band gap between the acoustic and optical branches, the lower bound is defined as the maximum acoustic frequency ω max ac normalized by ω 0 ; while its band width is equal to the difference between the optic minimum ω min op and the acoustic maximum ω max ac normalized by ω 0 if there does exist a gap. For this purpose, we have provided Fig. 9, that is, changes of width and starting frequencies of band gaps with mass and stiffness ratios are presented. The following trends can be concluded.
(i) Figures 9(a), 9(b), and 9(c) describe the effects of mass and stiffness ratios on the lower bounds of frequency. For a fixed α, three lower bounds increase with the increase in β. For a fixed β, it decreases with the increase in α. Consequently, when the mass ratio is smaller and the stiffness ratio is larger, the larger lower bound frequency can be obtained. Among the three, the third lower bound changes most abruptly in the studied parameter space, as shown in Fig. 9(c). Based on these diagrams, requested lower bounds can be achieved by selecting optimal α and β.
(ii) Figures 9(d), 9(e), and 9(f) depict the variation of every band gap width as the function of α and β. The first and second band widths increase and the third one decreases as β increases. For certain ranges of α and β, three band widths tend to become zero, suggesting the disappearance of the gaps. However, under the same conditions, the change to the third band gap is opposite, which means that the width of the third band gap tends to be a larger value. Relatively, the second band width has the widest adjustable range and has a maximum more than 4.  Thus, frequency levels and widths of band gaps can be programmed by adjusting mass and stiffness ratios. For example, a lattice with a lower M2 M1 = m2 m1 and a higher K2 K1 = K 2 K 1 possesses higher frequency band gaps. It serves as a guide to construct lattices to prevent the propagation of waves with specific frequencies.
6 Influence of orders of SG continua: 1D diatomic lattices with local interactions as an example Even though it is widely recognized that higher-order SG continua possess a higher accuracy, it is necessary to confirm whether unexpected issues may arise due to inclusions of higher SG terms, and determine up to which order of SG we can sufficiently capture satisfactory dispersion properties. To this end, we pick ZWLT series of continua to make a comparison. Because the only difference between ZWLT and the present model lies in whether to adopt the wavelengthdependent Taylor expansion, such a comparison is suitable to validate the proposed continuum. Four continua are examined, namely, fifth-, sixth-, eighth-, and tenth-order ones. Performances of continua are analyzed as follows.
(I) ZWLT continua show a strong dependence on SG orders. The increase in the SG order has two effects. On one hand, the workable range of continuum is extended by a higher order truncation. For example, the dispersion curves of the sixth-order continuum agree well with the discrete counterparts about in the range of the range of (0, 0.5); while those of the tenth-order do roughly in the range of (0, 0.7). On the other hand, with increasing SG orders, the continuum dispersion becomes more and more instable in the remining range of the first Brillouin zone. Here, the term "instability" means that the dispersive curves abruptly deviate from its discrete counterparts and can "penetrate" the band gaps predicted by the discrete model. Notably, such a kind of instability suggests that the corresponding continua fail to correctly predict band gap structures.
(II) When the wavelength of dependent Taylor expansion is employed, the continuum performance can be greatly enhanced. With increasing SG orders, the accuracy increases quickly with no instability issue. The fifth-and sixth-order continua already agree well with the discrete model except in a small range around ka 2π = 0.5. Such a continuum-discrete disparity can be explained as follows. According to Eq. (20), the maximum of d I noted as d max I is a. When ka 2π = 0.5 or equivalently λ = 2a, three positions x, x + a, and x + λ are equally spaced. Thus, calculating the displacement at x + d max I by Taylor expansion with reference to x or x + λ has an equal accuracy. Now the expansion accuracy is the lowest, since the distance between x + a and the reference position is the maximum among all wavelengths. Nevertheless, it is easy to overcome by increasing SG orders. When the order is increased to 8, the deviation can be eliminated, and the results of the present and discrete models are in much better agreement over the whole first Brillouin zone. Therefore, the eighth-order SG continuum is mainly adopted.

Lattices with various nonlocal interactions
In the above lattices studied, we have focused on behaviors of lattices with only one-neighbor interactions. In this section, we further investigate the SG model for lattices with both local and nonlocal interactions. Lattices with up to two-and three-neighbor interactions are analyzed to show the ability of the present model in dealing with nonlocal problems.
According to the derivation in Section 2, we obtain the dispersion equation of the discrete model with up to two-neighbor interactions, i.e., Eq. (10). To establish the equivalent SG continuum model, we select the unit cell 2 shown in Fig. 1(c). The periodic cell includes the whole atom located at 2n, two half atoms located at 2n−1 and 2n+1, respectively, two internal springs, two one-neighbor springs, two half two-neighbor springs, and one two-neighbor spring.
Once the kinetic and potential energy densities are figured out respectively, the EOMs of the corresponding SG continuum can be obtained by following the procedure in Section 3. By inserting Eq. (30) into these EOMs, we write the corresponding dispersion equation as where η 3 = (ikd I ) 2 + (ikdI) 4 12 + (ikdI) 6 360 + (ikdI) 8 20160 . Here, the value of d I in η 1 , η 2 , and η 3 is determined by Eq. (20).
The nonlocal interactions strongly affect dispersive characteristics and band gap distributions. Setting parameters α = 4, β = 2, α 0 = 2, β 0 = 0.5, γ = 2, and γ 0 = 1, we obtain the corresponding dispersion curves in Fig. 2(b). Dispersion curves of the present continuum model are still perfectly coincident with those of the corresponding discrete model. This indicates that the proposed model can effectively capture the dispersion behavior of periodic structures with nonlocal interactions when the wavelength-dependent Taylor expansion is adopted. Now curves become very different compared with the lattice simply with local interactions in Section 2. The fourth curve is no longer concave but convex, and the third curve is more convex, leading to the disappearance of the band gap between the third and fourth branches. Besides, the gap between the acoustic and optical branches also vanishes.
The existence of nonlocal interactions provides a new approach for designing desired band gaps. Figure 11 shows the effects of two-neighbor interactions on the band gap width and lower bound of each band gap. The lower bounds of the first and second band gaps vary mildly. The lower bound of the third band gap changes significantly and forms a maximum point with the increase in γ and γ 0 . From Fig. 11(d), within a wide range of γ and γ 0 , the first band gap tends to have a zero width, suggesting its disappearance. The maximum value of the first band gap width is far less than 1, while the values of the second band gap width are always close to 1. When γ > 1, the third gap width remains zero no matter how γ 0 changes. When γ < 1, the band gap becomes wider with decreasing γ. In brief, the nonlocal interactions play an important role in band gap structures of metamaterials.
To further explore the influence of nonlocal interactions, lattices with up to three-neighbor interactions are discussed. Following the derivation procedure in Section 2, we can get the where θ = K4 K1 , and K 4 is the stiffness coefficient of the third nearest neighbor spring. We select the unit cell 3 as shown in Fig. 1(d) to establish the SG continuum model by adopting the same procedure in Section 3. The corresponding dispersion relation can be finally expressed analytically. Notably, according to Fig. 5(b), the wavelength-dependent rule now becomes when λ ∈ a, 4a 3 , d I = a − λ, d II = 2(a − λ), when λ ∈ (2a, 4a], d I = a, d II = 2a − λ, (43) when λ ∈ (4a, +∞), d I = a, d II = 2a. (44) Inclusion of three-neighbor interactions makes the workable range of the long wave approximation narrow and makes the dispersion diagrams complex by causing more local stationary values. Setting parameters α = 4, β = 2, α 0 = 2, β 0 = 0.5, γ = 2, γ 0 = 1, and θ = 2, we can get the dispersion curves in Fig. 2(c). In the first Brillouin zone, the values of d I and d II are determined in a wavelength-dependent way, resulting in perfect agreement between the dispersion diagrams of the discrete and continuum models. The nonlocal interactions have greater effects on the third and fourth branches than the acoustic and optical ones. The third and fourth branches tend to be flattened [36] . Furthermore, band gap structures vary obviously compared with the case without nonlocal interactions. The continuum model is believed to have the capability of dealing with lattices with 1D periodic lattices with arbitrary nonlocal interactions.

Conclusions
With the Taylor expansion series of displacement modified, an exact SG model has been established for diatomic lattice metamaterials. The dispersion diagrams of continua coincide with those of discrete models in the first Brillouin zone. The wavelength-dependent Taylor expansion is adopted according to the periodicity of wave motion. To the best of our knowledge, the present continuum model is the first to successfully capture dispersive characteristics throughout the whole first Brillouin zone in various periodic lattices.
Some fundamental topics related to the present model have been studied, that is, the influence of different parameter ratios, appropriate truncation of SG orders, and local/nonlocal interactions. (I) The band gap can be designed by adjusting the parameters, and at this moment, the SG continuum model can still effectively predict the band gap characteristics and dispersion behavior of diatomic lattice metamaterials. The produced dispersion curves are essentially in agreement with those of the discrete model when the wavelength is close to the material scale. (II) The discussion also demonstrates that, by modifying the wavelengthdependent Taylor expansion, the prediction ability of the continuum model is enhanced, and by adopting higher-order terms, the accuracy of results is improved. Obviously, the addition of higher SG orders cannot improve in essence the capability for continua to capture dispersion relation, as the ZWLT models. In this work, it can be found that an eighth-order continuum is able to provide satisfactory results. (III) With the enhancement of nonlocal interactions, the practical range of long-wave approximation is narrowed, and the model becomes more complex. The suitable unit cell selection should also be considered, and thus the proposed continuum model can still approximate the mechanical behavior of metamaterials, presenting the consistent prediction. The results show that the shape of dispersion curves and band gap width are affected by nonlocal interactions.
In this study, we have focused on the dispersive characteristic of continua which is one key standard to judge the effectiveness of a continuum model. The assumption of infinite lattices spares us from discussing various kinds of boundary conditions which can be conveniently presented from the Hamilton formulations like Eq. (25). With a higher order included, the proposed SG model will be more accurate, while more complex boundary conditions arise. Thus, a balance between accuracy and complexity is important.
Although the 1D case is discussed here, the proposed methodology is also helpful in investigating features of wave propagations in plane and bulk media with periodic microstructures.