Nonlocal interaction engineering of 2D roton-like dispersion relations in acoustic and mechanical metamaterials

Abstract

The interior of the synthetic unit cells and their interactions determine the wave properties of metamaterials composed of periodic lattices of these cells. While local interactions with the nearest neighbors are well appreciated, nonlocal beyond-nearest-neighbor interactions are often considered as a nuisance. Here, by introducing a versatile effectively two-dimensional metamaterial platform for airborne sound and elastic waves, we exploit nonlocal effects as a powerful design tool. Within a simplified discrete model, we analytically show that the lowest band can be engineered by Fourier synthesis, where the $N$-th order Fourier coefficient is determined by the $N$-th nearest-neighbor interaction strength. Roton-like dispersion relations are an example. The results of the discrete model agree well with a refined model and with numerical calculations. In addition, we engineer the passage of waves from a local metamaterial into a nonlocal metamaterial by carefully tailoring the transition region between the two.

Introduction

In the solid-state physics of ordinary crystals, interactions way beyond the nearest atomic neighbors can be of crucial importance, for example for the crystal’s ionic binding energy¹. The Madelung constant² summarizes the effects of this long-range Coulomb interaction. In electromagnetism and optics, for artificial crystals called metamaterials^3,4, long-range interactions are well known to play an important role as well^5,6. For example, the dynamic electric dipole-dipole interaction asymptotically decays inversely with the distance between dipoles just like the static Coulomb potential^1,7. So far, however, long-range interactions in electromagnetic or optical metamaterials have mainly been considered a nuisance that complicates the physics rather than a useful design feature⁵.

In contrast, when conceptually introducing acoustical and optical phonons in ordinary crystals by mass-and-spring type models¹, the Hooke’s springs emerging from a given mass are usually only connected to the immediate neighbors of this mass^8,9,10,11. In few cases, springs between next-nearest neighbors have been considered as a correction^12,13. For many years, this only-nearest-neighbor spirit has also been taken for designing phonon dispersion relations in mechanical or acoustical metamaterials^{3,10,14,15,16,17}. Few exceptions have discussed the influence of long-range interactions in linear elastic materials^18,19,20 or in phononic crystals^21,22,23, albeit not in the spirit of a design tool but again rather as a correction to ordinary behavior.

Recently, we showed that suitable metamaterial designs using sufficiently strong $N$-th nearest-neighbor interactions with $N=3$, in addition to the usual nearest-neighbors interactions ($N=1$), allow for obtaining unusual acoustical-wave dispersion relations^24,25, which mimic the famous rotons in superfluid helium^26,27,28,29. However, the three-dimensional (3D) architectures discussed therein^24,25,30 were restricted in three regards, i)–iii). i) They were limited to tailoring the dispersion relation along only a single propagation direction. ii) They were not flexible enough to provide interactions way beyond $N=3$. iii) They could not easily be generalized to support more than just two types of interactions (nearest-neighbor and third-nearest-neighbor) simultaneously.

We note in passing that another recent publication³¹ used beyond-nearest-neighbor interactions in mechanical metamaterials to design and tailor the properties of higher phonon bands and states with the aim of obtaining topological band gaps. This work is not of immediate importance here because we focus on the lowest acoustical band.

Herein, we introduce and analyze by numerical and analytical calculations a flexible effectively two-dimensional metamaterial platform that allows for tailoring the lowest acoustical phonon dispersion relation for airborne sound by using beyond-nearest-neighbor interactions along two orthogonal directions. We design the long-range interactions by a network of cylindrical tubes, which connect cuboid compartments in a two-dimensional square array. The numerically calculated phonon bands reveal multiple roton-like minima along multiple directions. Through a simplified discrete analytical model for the lowest (acoustical) band that is mathematically equivalent to a mass-and-spring model, we show that the long-range interactions directly determine the Fourier coefficients of the acoustical-wave dispersion relation. Based on the simple model, we further demonstrate interesting negative refraction^32,33 and triple refraction at an interface between two metamaterials with and without beyond-nearest-neighbor interactions, respectively. We emphasize the importance of the transition region between an ordinary medium and the metamaterial comprising beyond-nearest-neighbor interactions as to which mode in the metamaterial the incident wave couples to. These calculations highlight the aspect of nonlocality^20,31, which is due to the beyond-nearest-neighbor interactions. We also propose a refined analytical model that can capture the higher bands in the band structure as well. The combination of these aspects shows that nonlocal effects are a powerful tool for designing acoustic and elastic metamaterials.

Results and discussion

Metamaterial design and numerical calculations

To implement the Fourier-synthesis idea in the simplest possible yet practical way, we consider airborne acoustical sound. Here, only longitudinally polarized pressure waves occur³⁴. We consider the air in tubes with rigid walls as mediator of the interactions. In acoustics, such hollow tubes do not lead to a finite minimum (“cut-off”) frequency, below which propagating waves do not occur. The behavior of airborne waves can easily be transferred to waterborne acoustical waves. In contrast, for waves in general elastic structures, two transverse modes occur in addition to the longitudinal modes. The situation would also be more complex for transversely polarized electromagnetic waves, for which tubes with walls made of a perfect electrical conductor do lead to a finite cut-off frequency³⁵. Such finite-frequency cut-off obviously inhibits an “acoustical” mode starting from zero frequency at zero wavenumber.

Figure 1a shows one unit cell of the suggested 3D structure of the metamaterial for airborne sound for the example of $N=3$. Panel b illustrates the resulting two-dimensional lattice. Panel c shows the unit cell for the fix cases $N=4,5,6,7,8$. In all of these cases, we consider a two-dimensional square lattice of hollow compartments with lattice constant $a$. The walls of these compartments are treated as rigid bodies. Mathematically, this corresponds to Neumann boundary conditions for the pressure field³⁴. Intuitively, the compartments can be seen as the “atoms”. Tubes connecting the compartments can be seen as mediators of the interactions among the atoms. Four thin tubes with inner radius ${r}_{1}$ directly connect each of these compartments with its four immediate neighbors. The walls of all tubes are treated as Neumann boundaries, too. These tubes mediate the nearest-neighbor interactions between the compartments. The thicker tubes with inner radius ${r}_{N} \, > \, {r}_{1}$ mediate the beyond-nearest-neighbor interactions with $N\ge 2$. To avoid collision of the tubes along the $x$- and the $y$-direction, respectively, one set of tubes is located above the ${xy}$-plane spanned by the thin tubes and the other one below. Therefore, the metamaterial structure itself has only two-fold rotational symmetry with respect to the $z$-axis. However, the structure exhibits an additional rotation-reflection symmetry with respect to the ${xy}$-plane. This symmetry ensures that the phonon dispersion relation ${\omega }_{n}({k}_{x},{k}_{y})$, which connects the angular frequency ${\omega }_{n}$ of the band with band index $n$ and wave vector ${{{{{\bf{k}}}}}}=({k}_{x},{k}_{y})$, has four-fold rotational symmetry, i.e., wave propagation along the $x$- and $y$-direction are strictly equivalent by symmetry. For a given fixed radius of the tubes, the maximum possible integer $N$ is clearly limited by geometrical constrains (cf. Fig. 1c). For the radii chosen in Fig. 1, the maximum possible value is $N=8$ as the tubes would overlap for $N=9$. However, conceptually, one can consider the limit of ${r}_{1}\to 0$ at fixed finite ratio ${r}_{N}/{r}_{1}$. In this limit, integers up to $N\to \infty$ are possible mathematically.

Fig. 1: Considered metamaterial platform supporting beyond-nearest-neighbor interactions along two orthogonal directions.

a Scheme of a single unit cell for $N=3$. b 2D square lattice with lattice constant $a$ built from this unit cell. The colors are for illustration only. The colored regions form a channel system that is bounded by rigid walls (not shown for clarity). The cuboid compartments (yellow) are connected to their four immediate neighbors along the $x$- and $y$-direction by tubes (yellow) with radius ${r}_{1}$. They are further connected to their $N$th nearest neighbors by tubes (blue) with radius ${r}_{N}$. The tubes along the $x$-direction ($y$-direction) are located above (below) the ${xy}$-plane spanned by the yellow tubes. c gallery of unit cells for $N=4$ to $N=8$ as indicated. For our calculations, we use the following fixed geometrical parameters: ${r}_{1}/a=0.03$, ${r}_{N}/a=0.05$, $w/a=0.30,h/a=0.80,{h}_{0}/a=0.24,{h}_{N}/a=0.25,$ and $d/a=0.86$. We choose $a=10{{{{{\rm{cm}}}}}}$ for the metamaterials supporting airborne sound and $a=100{{{{{\rm{\mu }}}}}}{{{{{\rm{m}}}}}}$ for the elastic metamaterials.

So far, we have only considered structures supporting nearest-neighbor interactions and $N$-th nearest-neighbor interactions with a single specific value of $N$. However, the metamaterial platform in Fig. 1b can easily be generalized to support several different beyond-nearest-neighbor interactions simultaneously. Figure 2 shows as an example a metamaterial structure supporting the three interactions with $N=1$, $N=3$, and $N=5$. The $N=3$ interactions lie in planes parallel to the ${xy}$-plane at different $z$-positions than the $N=5$ interactions. Clearly, for the chosen metamaterial platform, any further values of $N$ can be realized at yet different $z$-positions without imposing any geometrical constrains or difficulty (not depicted).

Fig. 2: Same as Fig. 1, but tubes (red) mediating the fifth-nearest-neighbor interactions (N = 5) are added as an example.

a Single metamaterial unit cell. b 2D square lattice composed thereof. Clearly, yet further values of $N$ could easily be implemented in other planes parallel to the ${xy}$-plane at other $z$-positions below and above the plane spanned by the yellow tubes ($N=1$). Again, due a rotation-reflection symmetry of the overall structure, wave propagation along the $x$- and $y$-direction is equivalent.

To compute the phonon dispersion relation of airborne acoustical waves propagating in the channel system shown in Figs. 1, 2, we consider the scalar wave equation for the spatial air-pressure modulation ${\widetilde{P}}_{{{{{{\bf{k}}}}}},i}({{{{{\bf{r}}}}}})$ on top of a constant background air pressure ${P}_{0}\gg {\widetilde{P}}_{{{{{{\bf{k}}}}}},i}({{{{{\bf{r}}}}}})$ at spatial position ${{{{{\bf{r}}}}}}=(x,y,z)$ and two-dimensional wave vector ${{{{{\bf{k}}}}}}=({k}_{x},{k}_{y})$¹$.$ We neglect damping. This assumption is expected to be reasonable if all absolute dimensions are sufficiently large³⁶. Furthermore, we have also performed calculations including damping (see Supplementary Note 1). The relative differences of the (real part of the) eigenfrequencies with respect to the lossless case are smaller than 4% for all cases (cf. Fig. 3 and Supplementary Figs. 1 and 2). The qualitative behavior is not changed at all. In the Fourier-domain, the resulting eigenvalue problem³⁴ reads

$${{{{{\boldsymbol{\nabla }}}}}}\cdot \left({{{{{\boldsymbol{\nabla }}}}}}{\widetilde{P}}_{{{{{{\bf{k}}}}}},n}\left({{{{{\bf{r}}}}}}\right)\right)=-\frac{{\omega }_{n}^{2}\left({{{{{\bf{k}}}}}}\right)}{{v}_{{{{{{\rm{air}}}}}}}^{2}}{\widetilde{P}}_{{{{{{\bf{k}}}}}},n}\left({{{{{\bf{r}}}}}}\right).$$

Fig. 3: Calculated phonon dispersion relation for the metamaterial structure shown in Fig. 1, i.e., for N = 3 beyond-nearest-neighbor interactions only.

a Tour through the first Brillouin zone (see inset) of the 2D square lattice with lattice constant $a$. The gray curve is obtained by numerical finite-element-method (FEM) calculations. The analytical results of the simple discrete model are shown by the red curve. Here, only the lowest band occurs. The results of the refined analytical model are shown by the blue dots. These almost completely overlap with the numerical calculations. The agreement with the simple discrete model is qualitatively very good, especially for the principal directions. The parameters for the gray and blue data are given in Fig. 1. The used parameters for the simple discrete model are: ${M=\rho {V}_{{{{{{\rm{c}}}}}}}=9.3\times {10}^{-5}{{{{{\rm{g}}}}}},{K}}_{1}=27.9{{{{{\rm{N}}}}}}\,{{{{{{\rm{m}}}}}}}^{-1}$, and ${K}_{3}=17.7{{{{{\rm{N}}}}}}\,{{{{{{\rm{m}}}}}}}^{-1}$. b Frequency as obtained from the numerical finite-element calculations represented on a false-color scale versus wavenumber ${k}_{x}$ and ${k}_{y}$ within the first Brillouin zone of the square lattice.

The band with index $n=1$ corresponds to an acoustical branch, the bands with $n\ge 2$ can be seen as “optical” phonon bands. We will focus our below design and discussion on the lowest band with $n=1$, but we will graphically also show several higher-frequency bands for completeness. We choose ${v}_{{{{{{\rm{air}}}}}}}=343{{{{{\rm{m}}}}}}/{{{{{\rm{s}}}}}}$ as the constant speed of sound in air. We solve this eigenvalue equation by a finite-element method (FEM) implemented in the commercial software Comsol Multiphysics. We assume Bloch periodic boundary conditions along the $x$- and the $y$-direction and treat the walls of all compartments and tubes as rigid immovable bodies via Neumann boundary conditions³⁵. All geometrical parameters of the metamaterial architecture are defined and given in Figs. 1, 2. As usual, the first 2D Brillouin zone of the 2D square lattice is given by the conditions $\left|{k}_{x}\right|\le \pi /a$ and $|{k}_{y}|\le \pi /a$.

Examples of calculated band structures are given in Figs. 3–5. Figure 3 shows the band structure of the metamaterial depicted in Fig. 1. It contains beyond-nearest-neighbor interactions with $N=3$. These 2D band structures can be compared with our previous quasi-1D band structures²⁴. Figure 3a shows the angular frequency for the usual tour through the 2D Brillouin zone (also see inset). Figure 3b represents the same results in 2D ${{{{{\bf{k}}}}}}$-space. It can be seen that we obtain local minima of the dispersion relation of the lowest acoustical band along the $({k}_{x},0)$ direction and the equivalent $(0,{k}_{y})$ direction. These local minima resemble the roton-like dispersion relations discussed in detail previously²⁴. In addition, we find four further (roton-like) local minima along the diagonals. As expected from the symmetry of the metamaterial structure (see above), the dispersion relation ${\omega }_{1}({k}_{x},{k}_{y})$ exhibits four-fold rotational symmetry. However, it is clearly not isotropic. The depth of the minima can be tailored by the effective strength of the beyond-nearest-neighbor interaction compared to the nearest-neighbor interactions, i.e., by the ratio of the tube radii ${r}_{3}/{r}_{1}$ (not depicted). The minimum is absent for ${r}_{3}/{r}_{1} \, = \, 0$, becomes deeper for increasing ratio ${r}_{3}/{r}_{1}\, > \,0$, and touches zero angular frequency, ${\omega }_{1}=0$, in the limit of ${r}_{3}/{r}_{1}\to \infty$.

Fig. 4: Dispersion relations as in Fig. 3a, but for different orders N of beyond-nearest-neighbor interactions in addition to the nearest-neighbor interactions.

a $N=3$. b $N=4$. c $N=5$. d $N=6$. e $N=7$. f $N=8$. The corresponding metamaterial structures are depicted in Fig. 1c.

Fig. 5: Calculated phonon dispersion relations for metamaterial structures as shown in Fig. 2.

Compared to the dispersion relations shown in Fig. 3, two different beyond-nearest-neighbor interactions are present simultaneously here. a–c correspond to $N=3$ and $N=4$ simultaneously, $N=3$ and $N=5$ simultaneously, and $N=3$ and $N=6$ simultaneously, respectively. Results from the numerical calculations (gray), the simple discrete model (red), and the refined analytical model (blue) are compared. d–f False-color representations of frequency versus wavenumber ${k}_{x}$ and ${k}_{y}$ from the numerical calculations. The behavior is due to the combined action of two different beyond-nearest-neighbor interactions.

Figure 4 shows the same as Fig. 3a, but for different values of $N$ instead of $N=3$. These results correspond to the metamaterial unit cells shown in Fig. 1. We find that the number of oscillations versus wave number within the first Brillouin zone increases with increasing $N$. This behavior already suggests a connection between the integer $N$ and the corresponding Fourier component. We will come back to this aspect in more detail in the following section.

Figure 5 is as Fig. 4 (where $N=3$ only), but for two different orders of long-range interaction in parallel, namely (a) $N=4$ in addition to $N=3$, (b) $N=5$ in addition to $N=3$, and (c) $N=6$ in addition to $N=3$. The corresponding metamaterial unit cells have already been shown in Fig. 2. Again, the angular frequency of the lowest acoustical band exhibits an oscillatory behavior versus wave number ${{{{{\bf{k}}}}}}$, however, now with two different Fourier components superimposed. We will come back to this aspect in more detail in the following section.

Mass-and-spring model

To obtain a more intuitive as well as an approximate analytical understanding of the findings of the previous section, we now consider a simple discrete analytical model. This model shows the acoustical dispersion relation of the metamaterial structures shown in Figs. 1, 2 can be tailored in the sense of Fourier synthesis. The beyond-nearest-neighbor interactions determine the Fourier coefficients.

We consider the mass-and-spring model shown in Fig. 6. We show that Newton’s equation of motion for the masses in this 2D lattice is mathematically equivalent to the approximate equation of motion for the air mass in the cuboid compartments shown in yellow in Fig. 1b. The hollow tubes connecting these compartments correspond to Hooke’s springs. This ansatz follows our previous reasoning^24,25 and is consistent with the numerically calculated air-pressure fields. An example referring to the parameters used in Fig. 3 is depicted in Supplementary Fig 3. We emphasize that this simple model accounts for the lowest (acoustical) band with $n=1$ only, as the model contains only one degree of freedom (the position of the single mass in the unit cell). We will expand the analytical modeling to also include the higher bands with $n\ge 2$ in the following section. Let us now derive the discrete model with a discrete set of equations of motion for the metamaterial structures illustrated in Figs. 1, 2.

Fig. 6: The shown 2D mass-and-spring model is mathematically equivalent to an approximate discrete model of the 2D metamaterial lattice shown in Fig. 1b.

To allow for comparison, the color of the masses (yellow) corresponds to that of the cuboid compartments in Fig. 1b. Likewise, the nearest-neighbor Hooke’s springs (yellow) correspond to the nearest-neighbor tubes in Fig. 1b. The beyond-nearest-neighbor Hooke’s springs with spring constants ${K}_{N}$, exemplified here for the case of $N=3$, are schematically illustrated by the blue lines. For clarity, they are only depicted for a single mass in the lattice. All masses $M$ in the square lattice with lattice constant $a$ are identical. Analytical dispersion relations obtained from this model are shown by the red dots in Figs. 3 and 4 (also see Supplementary Fig. 4). By virtue of the good agreement to the complete numerical calculations, the simple discrete model allows for discussing the highly unusual wave behavior at interfaces with reasonable effort (see Figs. 8–10).

Within the cuboid air compartment with volume ${V}_{{{{{{\rm{c}}}}}}}$ at the 2D lattice site defined by the pair of integers $(m,n)$, the air pressure ${P}_{{mn}}={P}_{0}+{\widetilde{P}}_{{mn}}$ shall be approximated by the constant mean pressure in that compartment. The air pressure directly translates into the number of air molecules ${N}_{{mn}}={N}_{0}+{\widetilde{N}}_{{mn}}$ in one compartment. ${N}_{0}$ is the number of air molecules in the compartment at fixed room temperature $T$, corresponding to the background pressure ${P}_{0}$. The ideal-gas equation becomes ${P}_{{mn}}={N}_{{mn}}{k}_{{{{{{\rm{B}}}}}}}T/{V}_{{{{{{\rm{c}}}}}}}$ or ${\widetilde{P}}_{{mn}}={\widetilde{N}}_{{mn}}{k}_{{{{{{\rm{B}}}}}}}T/{V}_{{{{{{\rm{c}}}}}}}$, with the Boltzmann constant ${k}_{{{{{{\rm{B}}}}}}}$. Furthermore, in each $N$th-nearest-neighbor horizontal and vertical tube with radius ${r}_{N}$ (see above), hence cross section $\pi {r}_{N}^{2}$, and length ${L}_{N}$, we approximate the air velocity along the tube axis as being constant throughout that tube. We define the corresponding velocities in the horizontal ($x$-direction) tubes, ${h}_{{mn}}^{(N)}$, and those in the vertical ($y$-direction) tubes, ${v}_{{ij}}^{(N)}$, mediating the $N$-th order interaction. The mass density, ${\rho }_{0}$, within all tubes is approximated as being constant, with ${\rho }_{0}={{m}_{0}P}_{0}/({k}_{{{{{{\rm{B}}}}}}}T)$, where ${m}_{0}$ is the mass of one air molecule. With these definitions, the continuity equation applied to compartment $(m,n)$ describes the in- and out-flux of air molecules from the tubes into the compartment and reads

$$\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}t}{{m}_{0}N}_{{mn}}=\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}t}{m}_{0}{\widetilde{N}}_{{mn}}=\mathop{\sum }\limits_{N=1}^{{{{{\infty }}}}}{-{\rho }_{0}\pi r}_{N}^{2}\left(\left({h}_{{mn}}^{(N)}-{h}_{m-N,n}^{(N)}\right)+\left({v}_{{mn}}^{(N)}-{v}_{m,n-N}^{(N)}\right)\right).$$

(1)

The acceleration $\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}t}{h}_{{mn}}^{(N)}$ in the $N$th-order horizontal tubes results from the net force corresponding to the pressure difference between the two ends of the tube with length ${L}_{N}$, i.e., from

$${\rho }_{0}\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}t}{h}_{{mn}}^{(N)}=-\frac{{\widetilde{P}}_{m+N,n}-{\widetilde{P}}_{{mn}}}{{L}_{N}},$$

(2)

and likewise for the $N$th-order vertical tubes

$${\rho }_{0}\frac{{{{{{\rm{d}}}}}}}{{{{{{\rm{d}}}}}}t}{v}_{{mn}}^{(N)}=-\frac{{\widetilde{P}}_{m,n+N}-{\widetilde{P}}_{{mn}}}{{L}_{N}}.$$

(3)

Taking the time derivative of Eq. (1), inserting Eqs. (2) and (3) into Eq. (1), and replacing ${\widetilde{N}}_{{mn}}=\frac{{\rho }_{0}{V}_{{{{{{\rm{c}}}}}}}}{M{P}_{0}}{\widetilde{P}}_{{mn}}$ leads to.

$$M\frac{{{{{{{\rm{d}}}}}}}^{2}}{{{{{{\rm{d}}}}}}{t}^{2}}{\widetilde{P}}_{{mn}}=\mathop{\sum }\limits_{N=1}^{{{{{\infty }}}}}{K}_{N}\left(\left({\widetilde{P}}_{m-N,n}-2{\widetilde{P}}_{{mn}}+{\widetilde{P}}_{m+N,n}\right)+\left({\widetilde{P}}_{m,n-N}-2{\widetilde{P}}_{{mn}}+{\widetilde{P}}_{m,n+N}\right)\right),$$

(4)

where we have introduced the two abbreviations

$$M={\rho }_{0}{V}_{{{{{{\rm{c}}}}}}}$$

(5)

and

$${K}_{N}={P}_{0}\frac{\pi {r}_{N}^{2}}{{L}_{N}}.$$

(6)

Equation (4) is identical to Newton’s equation for the mass-and-spring model shown in Fig. 6 if the pressure variation ${\widetilde{P}}_{m,n}$ in Eq. (4) is replaced by the out-of-plane displacement. $M$ in Eq. (4) is the air mass in one compartment. It is the same for all lattice sites at this point. Below, when addressing interfaces, we will consider different values of $M$ at different locations. According to Eq. (5), the mass $M$ can be varied in practice by changing the volume of the compartment ${V}_{{{{{{\rm{c}}}}}}}$. ${K}_{N}$ is the effective Hooke’s spring constant for the $N$th order interaction. According to Eq. (6), it can be tailored by the radius of the $N$th order tube, ${r}_{N}$. Obviously, at fixed tube radius ${r}_{N}$, the effective spring constant ${K}_{N}$ decreases inversely versus increasing order $N$ of the interaction, hence with increasing length ${L}_{N}\approx N{L}_{1}$ of the tube for $N\ge 2$. This behavior is perfectly analogous to that of an ordinary elastic Hooke’s spring. As the tube diameter $2{r}_{N}$ is geometrically constrained in size from above (cf. Fig. 1), and must eventually decrease with increasing $N$, the effective relative strength of the $N$th-order interaction eventually tends to zero inversely proportional to the order $N$, or stronger than that. In the previous section, we have seen that values up to $N=8$ can be achieved in practice.

Making a plane-wave ansatz for Eq. (4) according to

$${\widetilde{P}}_{{mn}}=\widetilde{P}{{\exp }}\left({{{{{\rm{i}}}}}}({{mk}}_{x}a+n{k}_{y}a-\omega t)\right),$$

(7)

with constant amplitude prefactor $\widetilde{P}$, we obtain the dispersion relation for acoustical pressure waves with band index $n=1$ in the discrete model

$${\omega }_{1}^{2}\left({k}_{x},{k}_{y}\right)=\mathop{\sum }\limits_{N=0}^{{{\infty }}}{F}_{N}{{\cos }}\left(N{k}_{x}a\right)+\mathop{\sum }\limits_{N=0}^{{{\infty }}}{F}_{N}{{\cos }}\left(N{k}_{y}a\right).$$

(8)

This expression is central to our paper. We see that, for $N\ge 1$, the Fourier coefficient ${F}_{N}$ in the two Fourier series on the right-hand side is directly connected to the $N$th nearest-neighbor effective Hooke’s spring constant ${K}_{N}$ and is given by

$${{F}_{N}=-2K}_{N}/M.$$

(9)

The constant $0$th-order term of the Fourier series with

$${F}_{0}=\mathop{\sum }\limits_{{N}^{{\prime} }=1}^{{{{{\infty }}}}}2{K}_{{N}^{{\prime} }}/M$$

(10)

guarantees that the dispersion relation Eq. (7) starts according to

$${\omega }_{1}({k}_{x},\,{k}_{y})\propto \left|{{{{{\bf{k}}}}}}\right|$$

(11)

in the limit $\left|{{{{{\bf{k}}}}}}\right|\to 0$. This is equivalent to saying that we restrict ourselves to “acoustical” dispersion relations. In the Fourier series on the right-hand side of Eq. (8), only cosine terms appear as sine terms would generally violate the reciprocity condition $\omega (-{{{{{\bf{k}}}}}})=\omega ({{{{{\bf{k}}}}}})$.

For suitable parameters, the above simple 2D discrete model again exhibits roton-like dispersion relations along multiple directions (see Supplementary Note 2 and Supplementary Fig 4). To allow for a direct comparison with the above numerical calculations, we plot the calculated dispersion relations of this model (red dots) together with those of the numerical calculations (gray curves) in Figs. 3 and  4. For the lowest band, the overall agreement is excellent for all values of $N$. This especially holds true for the principal direction ${{{{{\mathbf{\Gamma }}}}}}{{{{{\bf{M}}}}}}$. Quantitative discrepancies occur for other directions, but the agreement is still qualitatively good. Note, however, that we have taken the freedom to use fitted ratios ${K}_{N}/M$. In principle, these ratios could be calculated from Eqs. (5) and (6) and from the known geometrical parameters therein. The fitting procedure gives us freedom to correct for the simplicity of the mass-and-spring model. We will come back to a more stringent treatment in the following section.

Before we get there, we now apply the mass-and-spring model to treating metamaterials that are not infinitely extended, but that rather contain interfaces. The treatment of such configurations would be computationally very expensive on the level of the complete numerical calculations.

As a representative example, we consider only nearest-neighbor interactions and third-nearest-neighbor interactions ($N=3$). In such a system²⁴, the dispersion relation does not monotonically increase but features a region of backward waves, with negative scalar product of the phase and group velocity vectors, and a roton-like^26,28 local minimum in the dispersion relation. We expect interesting wave behavior, including negative refraction and triple refraction, which we shall study in the remainder of this section.

The considered interface and excitation configuration is illustrated in Fig. 7. In the region above the interface, we consider a discrete system simultaneously with nearest-neighbor interactions and third-nearest-neighbor interactions ($N=3$). In the region below the interface, we assume only nearest-neighbor interactions. We force displacements of masses located on the indicated black line according to ${{\cos }}\left(2\pi {f}_{{{{{{\rm{c}}}}}}}t\right){{\exp }}(-{\left(\pi {f}_{{{{{{\rm{c}}}}}}}t/50\right)}^{2}){{\exp }}(-({\left(x-{x}_{{{{{{\rm{c}}}}}}}\right)}^{2}+{\left(y-{y}_{{{{{{\rm{c}}}}}}}\right)}^{2})/{37.5}^{2})$, with ${x}_{{{{{{\rm{c}}}}}}}=-100$, ${y}_{{{{{{\rm{c}}}}}}}=-130$ being the coordinate of the middle point of the line in units of the lattice constant $a$, and ${f}_{{{{{{\rm{c}}}}}}}$ being the carrier frequency of the temporally Gaussian excitation signal. This excitation launches a spatially Gaussian wave packet towards the interface at an oblique angle of ${45}^{{{{{{\rm{o}}}}}}}$. We have chosen this particular angle just as an example to illustrate the principle. The temporal Gaussian envelope contains about 50 oscillation cycles (see Fig. 7c), corresponding to a fairly narrow spread in the frequency domain (see Fig. 7d).

Fig. 7: Simulation setup of wave refraction at interface.

a Wave transmission at an interface between a local system with only nearest-neighbor interactions (bottom, gray) and a nonlocal system additionally including third-nearest-neighbor interactions, ${K}_{3}\,\ne\, 0$ (top, yellow). The incident wave impinges under an angle of ${45}^{^\circ }$ onto the interface. b Due to the mismatch of interaction orders between the two regions, some of the third-order Hooke’s springs are cut away near the interface (dashed horizontal line) within the nonlocal region. To control the transmission of the incident waves, we fine-tune the mass parameters, ${M}_{01}$ and ${M}_{02}$, and three Hooke’s spring constants, ${K}_{01},$ ${K}_{02}$, and ${K}_{03}$, at the interface. To avoid reflection of the partial waves at the boundaries of the simulation domain, we choose a much larger simulation domain than the one shown. A wave packet that is Gaussian in space and time is launched by prescribing the pressure in Eq. (4) on the black lines according to ${{\cos }}\left(2\pi {f}_{{{{{{\rm{c}}}}}}}t\right){{\exp }}(-{\left(\pi {f}_{{{{{{\rm{c}}}}}}}t/50\right)}^{2}){{\exp }}(-({\left(x-{x}_{{{{{{\rm{c}}}}}}}\right)}^{2}+{\left(y-{y}_{{{{{{\rm{c}}}}}}}\right)}^{2})/{37.5}^{2})$. c Temporal profile of the launched pulse, containing about $50$ oscillation cycles. d Corresponding Fourier transform (absolute value).

Due to the mismatch of interaction orders between the two regions at the interface, one of the Hooke’s springs that mediate the third-nearest-neighbor interactions for the first three masses above the interface (see Fig. 7b) is truncated by the interface. Therefore, the connections at the interface are modified slightly (see Fig. 7b). These modifications in the transition region between the two media will turn out to be crucial for controlling the transmitted and reflected partial waves. As shown in the below simulations, which are based on a trial-and error procedure for the design and for the parameters of the transition layer, we can even achieve nearly $100 \%$ transmission for each of the three possible modes individually by fine-tuning the indicated interface parameters, including the two mass parameters, ${M}_{01}$ and ${M}_{02}$, and three spring constants, ${K}_{01}$, ${K}_{02}$ and ${K}_{03}$ in Fig. 7. Alternatively, the parameters for the transition region can be chosen such that an incident wave couples to all three refracted modes simultaneously (see Supplementary Fig 7).

We note in passing that an explicit and rational construction procedure for the transition layer between a local and a non-local medium and its parameters is presently elusive to our knowledge. Defining such a procedure beyond a trial-and-error approach is beyond the scope of our present paper.

In Fig. 8, we first demonstrate negative refraction at the interface. For the incidence region, we choose the parameters $M/{M}_{0}=1.0,{K}_{1}/{K}_{0}=1.2$, and ${K}_{3}/{K}_{0}=0$. Here, ${M}_{0}$ and ${K}_{0}$ are not relevant to the results and can take any constant reference values. In Fig. 8a, red (blue) solids lines represent iso-frequency contours for the discrete system in the incidence (transmission) region at carrier frequency (= center frequency) $\omega /{\omega }_{0}=0.7$ of the Gaussian wave packed used in the numerical simulation. Dashed lines correspond to a slightly larger frequency. From two nearby iso-frequency contours, we can identify the group velocity vector ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}}$. From Snell’s law, one negatively refracted mode with wave vector ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}}$ and group velocity ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}}$ is expected. A snapshot of the simulated displacement fields at one instant in time, $t/{T}_{{{{{{\rm{c}}}}}}}=150$, is shown in Fig. 8b. The interface parameters, as indicated in the figure, are tuned to achieve a nearly zero reflected partial wave. As expected, the refracted waves exactly move along the direction of the predicted group velocity ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}}$. The phase fronts of the refracted waves move along the wave vector ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}}$. The dashed lines are guides to the eye for the wave propagation trajectory. The negative refraction can be seen yet more clearly in Supplementary Movie 1.

Fig. 8: Negative refraction at an interface as illustrated in Fig. 7.

a Iso-frequency curves for the nonlocal system, with the same parameters as in Supplementary Fig 4, and the local system, with parameters $M/{M}_{0}=1,{K}_{1}/{K}_{0}=1.2$, and ${K}_{3}/{K}_{0}=0$. The solid lines correspond to the carrier frequency ${\omega }_{{{{{{\rm{c}}}}}}}/{\omega }_{0}=0.7$ of the Gaussian pulse. The blue arrow represents the incident wave vector, ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{i}}}}}}}$. From Snell’s law, a single negatively refracted mode with wave vector, ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}}$, and group velocity vector, ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}}$, is possible. b Snapshot of simulated pressure fields at $t/{t}_{{{{{{\rm{c}}}}}}}=150$, with ${t}_{{{{{{\rm{c}}}}}}}=2\pi /{\omega }_{{{{{{\rm{c}}}}}}}$. The interface parameters (cf. Fig. 7) are ${M}_{01}/{M}_{0}=0.5,{M}_{02}/{M}_{0}=1.0,{K}_{01}/{K}_{0}=1.0,{K}_{02}/{K}_{0}=1.0$, and ${K}_{03}/{K}_{0}=1.0$. The black dashed straight lines indicate the propagation path of the incident and the refracted wave. In the nonlocal region, the wave propagation direction agrees well with the group velocity vector as derived in a. An animated version of b is provided in Supplementary Movie 1. The displacement component ${u}_{z}$ for the two cuts perpendicular to the propagation direction at positions A (brown) and B (blue), respectively, is shown in Supplementary Fig 5.

Next, we demonstrate triple refraction. In this case, we choose different parameters for the incidence region, namely $M/{M}_{0}=0.44,{K}_{1}/{K}_{0}=2.3$, and ${K}_{3}/{K}_{0}=0$. The carrier frequency of the Gaussian packet is chosen as ${\omega }_{{{{{{\rm{c}}}}}}}/{\omega }_{0}=0.9$. Iso-frequency contours are depicted in Fig. 9 the same manner as in Fig. 8a. For this case, there are three possible refracted modes, with the wave vectors and the group velocities, ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}1}$ and ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}1}$, ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}2}$ and ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}2}$, and ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}3}$ and ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}3}$, respectively, instead of a single mode as in the previous example. The third mode is again a backward wave with ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}3}\cdot {{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}3} \, < \, 0$. The importance of the interface parameters becomes evident from the three different calculations shown in Fig. 10a–c. In each of the three cases, one of the three different modes is excited almost exclusively. This means that one can select a wanted mode by appropriately choosing the parameters of the transition region between the two media. For all three cases, we further illustrate the behavior by the Supplementary Movies 2–4.

Fig. 9: Triple refraction at an interface as illustrated in Fig. 7.

Iso-frequency contours for the nonlocal system, with the same parameters as in Supplementary Fig 4, and the local system, with parameters $M/{M}_{0}=0.44,{K}_{1}/{K}_{0}=2.3$, and ${K}_{3}/{K}_{0}=0$. The solid lines correspond to the carrier frequency$\,{{\omega }_{{{{{{\rm{c}}}}}}}/\omega }_{0}=0.9$ of the excited Gaussian wave packet. By using Snell’s law, triple refraction with wave vectors, ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}1}$, ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}2}$ and ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}3}$, and corresponding group velocity vectors, ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}1}$, ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}2}$ and ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}3}$, becomes possible. The partial wave with wave vector ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}3}$ and group velocity vector ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}3}$ is a backward wave with ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}3}\cdot {{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}3} \, < \, 0$.

Fig. 10: Relative amplitude control of the three refracted modes at the interface between a nonlocal system and a local system by tuning the interface parameters illustrated in Fig. 7b.

a Snapshot of simulated pressure fields with the interface parameters, ${M}_{01}/{M}_{0}=0.5$, ${M}_{02}/{M}_{0}=1.0$, ${K}_{01}/{K}_{0}=6.0$, ${K}_{02}/{K}_{0}=1.0$, and ${K}_{03}/{K}_{0}=1.0$, respectively. In this case, only the refracted mode with the wave vector ${{{{{{\bf{k}}}}}}}_{{{{{{\rm{t}}}}}}1}$ and group velocity vector ${{{{{{\bf{v}}}}}}}_{{{{{{\rm{g}}}}}}1}$ occurs. The dashed lines indicate the propagation path of the incident and refracted waves. b, c Same as a, but with interface parameters ${M}_{01}/{M}_{0}=0.8$, ${M}_{02}/{M}_{0}=0.5$, ${K}_{01}/{K}_{0}=2.0$, ${K}_{02}/{K}_{0}=0.2$, ${K}_{03}/{K}_{0}=6.0$, and ${M}_{01}/{M}_{0}=1.0$, ${M}_{02}/{M}_{0}=0.6$, ${K}_{01}/{K}_{0}=1.0$, ${K}_{02}/{K}_{0}=1.0$, and ${K}_{03}/{K}_{0}=4.0$, respectively. Animated versions of the three scenarios are provided in the Supplementary Movies 2–4. The displacement component ${u}_{z}$ for the cuts perpendicular to the propagation direction at the positions A (brown) and B (blue) in a–c is shown in Supplementary Fig 6.

We note in passing that the rich physics of the transition region can likely also be captured by generalizing the effective-medium approximation of the mass-and-spring model in terms of a generalized wave equation for homogeneous non-local media, as introduced in ref. ²⁵, to the case of heterogeneous non-local elastic media. However, this aspect goes well beyond the scope of the present paper.

Refined analytical model

In this section, we present a refined analytical model, which captures not only the behavior of the lowest band but also that of the higher bands for the proposed metamaterials. We make the following two assumptions: i) The acoustic pressure in the cuboid compartments in Fig. 1 is approximated as being constant within. ii) The pressure modes in all tubes are the fundamental waveguide mode, i.e., at a given position along the tubes, the pressure is constant over the cross section³⁴. These two assumptions are expected to be valid for low-frequency sound, for which the wavelength in air is much larger than the unit cell size.

On this basis, we derive the metamaterial dispersion relations analytically based on the Floquet-Bloch theorem¹. As an ansatz, the acoustic pressure modulation at the lattice site defined by the pair of integers $(m,n)$ is given by ${\widetilde{P}}_{{mn}}=\widetilde{P}\,{{\exp }}({{{{{\rm{i}}}}}}({k}_{x}{ma}+{k}_{y}{na}-\omega t))$, with $\widetilde{P}$ being a constant prefactor as in Eq. (7). The key point to derive the dispersion relation is to analyze acoustic pressure wave propagation in the tubes. As one example, we consider the tube that connects the cuboid compartment at site $\left(m,{n}\right)$ to its $N$th nearest neighbor along $x$-direction at site $(m+N,{n})$. The pressure field inside the tube is composed of a forward waveguide mode, from the cuboid at site $\left(m,{n}\right)$ to that at site $\left(m+N,{n}\right)$, and a corresponding backward waveguide mode

$$p(s)=\left({A}^{+}{{\exp }}\left({{{{{\rm{i}}}}}}\frac{\omega }{{v}_{{{{{{\rm{air}}}}}}}}s\right)+{A}^{-}{{\exp }}\left(-{{{{{\rm{i}}}}}}\frac{\omega }{{v}_{{{{{{\rm{air}}}}}}}}s\right)\right){{\exp }}\left(-{{{{{\rm{i}}}}}}\omega t\right).$$

(12)

Herein, $s$ indicates the distance along the central axis of the tube, with $s=0$ corresponding to one end of the tube, where it connects the cuboid at site $\left(m,{n}\right)$. We further have $s={L}_{N}$, with ${L}_{N}=d-w+2{h}_{N}+\sqrt{{d}^{2}+{N}^{2}{a}^{2}}$ being the tube length, representing the other end. ${A}^{+}$ and ${A}^{-}$ are two unknown amplitude coefficients. The corresponding particle velocity is derived as

$$v\left(s\right)=\frac{1}{{{{{{\rm{i}}}}}}\omega {\rho }_{0}}\frac{\partial p}{\partial s}=\frac{1}{{\rho }_{0}{v}_{{{{{{\rm{air}}}}}}}}\left({A}^{+}{{\exp }}\left({{{{{\rm{i}}}}}}\frac{\omega }{{v}_{{{{{{\rm{air}}}}}}}}s\right)-{A}^{-}{{\exp }}\left(-{{{{{\rm{i}}}}}}\frac{\omega }{{v}_{{{{{{\rm{air}}}}}}}}s\right)\right){{\exp }}\left(-{{{{{\rm{i}}}}}}\omega t\right),$$

(13)

with ${\rho }_{0}$ representing the average air density as previously. Due to continuity, the acoustic pressure at the two ends of the tube must be the same as in the corresponding cuboids³⁴

$$p\left(0\right)=\widetilde{P}\,{{\exp }}\left({{{{{\rm{i}}}}}}{k}_{x}{ma}+{{{{{\rm{i}}}}}}{k}_{y}{na}-{{{{{\rm{i}}}}}}\omega t\right),$$

(14)

$$p\left({L}_{N}\right)=\widetilde{P}\,{{\exp }}({{{{{\rm{i}}}}}}{k}_{x}(m+N)a+{{{{{\rm{i}}}}}}{k}_{y}{na}-i\omega t).$$

(15)

The two amplitude coefficients, ${A}^{+}$ and ${A}^{-}$, result from

$${A}^{+}=\frac{\widetilde{P}}{2i}{{\exp }}\left({{{{{\rm{i}}}}}}{k}_{x}{ma}+{{{{{\rm{i}}}}}}{k}_{y}{na}\right){{\csc }}\left(\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)\left({{\exp }}({{{{{{\rm{i}}}}}}k}_{x}{Na})-{{\exp }}\left(-{{{{{\rm{i}}}}}}\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)\right),$$

(16)

$${A}^{-}=\frac{\widetilde{P}}{2i}{{\exp }}\left({{{{{\rm{i}}}}}}{k}_{x}{ma}+{{{{{\rm{i}}}}}}{k}_{y}{na}\right){{\csc }}\left(\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)\left(-{{\exp }}\left({{{{{{\rm{i}}}}}}k}_{x}{Na}\right)+{{\exp }}\left({{{{{\rm{i}}}}}}\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)\right).$$

(17)

Substitution of the above two equations into Eq. (13) leads to the particle velocity

$$v(s)=\frac{{{{{{\rm{i}}}}}}\,\widetilde{P}}{{\rho }_{0}{v}_{{{{{{\rm{air}}}}}}}}{{\csc }}\left(\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right){{\exp }}\left({{{{{\rm{i}}}}}}{k}_{x}{ma}+{{{{{\rm{i}}}}}}{k}_{y}{na}-i\omega t\right)\\ \left({{\cos }}\left(\frac{\omega s-\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)-{{\cos }}\left(\frac{\omega s}{{v}_{{{{{{\rm{air}}}}}}}}\right){{\exp }}\left({{{{{{\rm{i}}}}}}k}_{x}{Na}\right)\right).$$

(18)

From this expression, the total air mass flowing away from the cuboid at lattice site $\left(m,{n}\right)$ through the above tube is given by

$${Q}_{m+N,n}={\rho }_{0}{S}_{N}v\left(0\right)={{{{{\rm{i}}}}}}\frac{{S}_{N}{\widetilde{P}}_{{mn}}}{{v}_{{{{{{\rm{air}}}}}}}}{{\csc }}\left(\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)\,\left({{\cos }}\left(\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)-{{\exp }}\left({{{{{{\rm{i}}}}}}k}_{x}{Na}\right)\right).$$

(19)

Here, the area ${S}_{N}=\pi {R}_{N}^{2}$ represents the cross section of the tube and $v(0)$ the particle velocity at one end of the tube. Similarly, we derive the average mass flow away from the cuboid at lattice site $\left(m,{n}\right)$ through the tube that connects the two cuboids at sites $\left(m,{n}\right)$ and $\left(m,{n}+N\right)$ as

$${Q}_{m,n+N}={{{{{\rm{i}}}}}}\frac{{S}_{N}{\widetilde{P}}_{{mn}}}{{v}_{{{{{{\rm{air}}}}}}}}{{\csc }}\left(\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)\,\left({{\cos }}\left(\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)-{{\exp }}\left({{{{{{\rm{i}}}}}}k}_{y}{Na}\right)\right).$$

(20)

With the conservation law for the air mass inside the cuboid compartment at site $\left(m,{n}\right)$, we have³⁴

$${{{{{\rm{i}}}}}}\frac{\omega }{{v}_{{{{{{\rm{air}}}}}}}^{2}}{V}_{{{{{{\rm{c}}}}}}}{\widetilde{P}}_{{mn}}=\mathop{\sum}\limits_{N=1}^{{{{{\infty }}}}}\left({Q}_{m+N,n}+{Q}_{m-N,n}+{Q}_{m,n+N}+{Q}_{m,n-N}\right).$$

(21)

${Q}_{m-N,{n}}$ is the mass flow through the tube that connects the two cuboids at sites $(m,n)$ and $(m-N,n)$. ${Q}_{m,{n}-N}$ is defined analogously. After some mathematical simplifications, we arrive at the following expression for the dispersion relation

$$\frac{\omega }{{v}_{{{{{{\rm{air}}}}}}}}{V}_{c}=\mathop{\sum}\limits_{N=1}^{{{\infty }}}\left(\pi {R}_{N}^{2}{{\csc }}\left(\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)\left(4{{\cos }}\left(\frac{\omega {L}_{N}}{{v}_{{{{{{\rm{air}}}}}}}}\right)-2{{\cos }}\left({k}_{x}{Na}\right)-2{{\cos }}({k}_{y}{Na})\right)\right).$$

(22)

${V}_{{{{{{\rm{c}}}}}}}$ represents the volume of the compartment. This equation connects the angular frequency $\omega$ with ${k}_{x}$ and ${k}_{y}$. For a given Bloch wave vector ${{{{{\bf{k}}}}}}\,{{{{{\boldsymbol{=}}}}}}{{{{{\boldsymbol{(}}}}}}{k}_{x},\,{k}_{y}{{{{{\boldsymbol{)}}}}}}$, the implicit formula Eq. (22) provides multiple solutions for $\omega$. These solutions correspond to the different bands. Compared to the previous section, the price we pay for the refined modeling in the present section is that we do not obtain a closed explicit expression for the angular frequency $\omega$ as a function of the Bloch wave vector. However, the implicit Eq. (22) can easily be solved numerically. Results are shown by the blue dots in Figs. 3 and 4. They agree well with the numerical FEM calculations (gray curves) regarding the lowest bands as well as the higher bands. For the case of two different orders of long-range interaction in parallel shown in Fig. 5a–c, the qualitative agreement for the lowest bands is again good, but quantitative differences occur, especially along the ${{{{{\bf{MK}}}}}}$ direction. The differences are yet more pronounced for the higher bands. We assign these differences to the fact that the vertical parts of the tubes for different $N$ partially overlap (cf. Fig. 2, part with height ${h}_{3}$), leading to interference of the respective partial waves. The refined model neglects these interferences. However, we recall that the focus of this study lies on the lowest bands. For these, the refined model generally performs better than the simple discrete model in Figs. 3–5.

Elastic waves instead of airborne sound

So far, we have exclusively discussed longitudinally polarized airborne sound waves. It is interesting to investigate whether we can translate the overall approach of the metamaterial platform illustrated in Figs. 1 and 2 to elastic waves, for which three modes emerge from the ${{{{{\boldsymbol{\Gamma }}}}}}$ point rather than just a single mode for airborne sound. For elastic waves, the structures in Figs. 1 and 2 should be interpreted as being composed of a single ordinary elastic material (e.g., a polymer, see Supplementary Note 3) rather than as air channels. Supplementary Fig 8 provides numerically calculated example band structures for $N=3,4$, and $5$, that can be compared with Fig. 4a–c. Clearly, the overall behavior of the lowest elastic band is closely similar to that of the lowest band for airborne sound. Specifically, the number of extrema of $\omega({{{{{\bf{k}}}}}})$, with ${{{{{\bf{k}}}}}}=(k,0)$ or ${{{{{\bf{k}}}}}}=(0,k)$, versus $k$ in the interval $k\in[0,\pi/a]$ is equal to $N$—in agreement with the Fourier synthesis idea presented above.

It is presently not clear whether this idea can also be translated to experimentally accessible metamaterial structures for yet other types of waves, such as, e.g., surface water waves, microwaves or light waves. We hope that our work stimulates future experiments and further design studies in this direction.