# Closed-form solutions for wave propagation in hexagonal diatomic non-local lattices

F. Ongaro, P.H. Beoletto, F. Bosia, M. Miniaci, N.M. Pugno

► To cite this version:

F. Ongaro, P.H. Beoletto, F. Bosia, M. Miniaci, N.M. Pugno. Closed-form solutions for wave propagation in hexagonal diatomic non-local lattices. International Journal of Mechanical Sciences, 2025, 291-292, pp.110095. 10.1016/j.ijmecsci.2025.110095. hal-05014235

HAL Id: hal-05014235

https://hal.science/hal-05014235v1

Submitted on 20 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Closed-form solutions for wave propagation in hexagonal diatomic non-local lattices

F. Ongaro<sup>a,*</sup>, P.H. Beoletto<sup>b</sup>, F. Bosia<sup>b</sup>, M. Miniaci<sup>c</sup>, N.M. Pugno<sup>a,d,**</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Mass-spring systems

Elastic waves

Non-neighbouring interactions

Roton-like dispersion diagram

Topologically-protected waveguiding

Bound modes

## ABSTRACT

Periodic mass-spring lattices are commonly used to investigate the propagation of waves in elastic systems, including wave localisation and topological protection in phononic crystals and metamaterials. Recent studies have shown that introducing non-neighbouring (i.e., beyond nearest neighbour) connections in these chains leads to multiple topologically localised modes, while generating roton-like dispersion relations. This paper focuses on the theoretical analysis of elastic wave propagation in hexagonal diatom mass-spring systems in which both neighbouring and non-neighbouring interactions occur through linear elastic springs. Closed-form expression for the dispersion equations are derived, up to an arbitrary order of beyond-the-nearest connections for both in-plane and out-of-plane mass displacements. This allows to explicitly determine the influence of the order of non-neighbouring interactions on the band gaps, the local minima and the slope inversions in the first Brillouin zone for the considered unit cell. All analytical solutions are numerically verified. Finally, examples are provided on how non-neighbouring connections can be exploited to enhance the localisation of topologically-protected edge modes in waveguides constructed using mirror symmetric diatomic lattices constituted by two regions with different unit cell orientations. The study provides further insight on how to design phononic crystals generating roton-like behaviour and to exploit them for topologically protected waveguiding.

## 1. Introduction

Over the past twenty years, there has been a surge of interest in the unique wave control properties of photonic and phononic crystals (PCs) and metamaterials (MMs) [1–8]. Originally explored in electromagnetism, their working principles were quickly extended to the field of elasticity [9–11], demonstrating great potential in a wide range of applications spanning different spatial and frequency ranges, i.e., from micro-scale surface acoustic waves [12] and microelectromechanical devices [13] to macro-scale systems for noise and vibration reduction [14,15], as well as large-scale seismic metamaterials for vibration shielding of buildings [16–19], and infrastructures [20].

In more recent years, topological MMs [21,22] have emerged as a topic of particular interest, allowing to replicate in macroscopic systems well-known effects in condensed matter theory such as Quantum Hall [23], Quantum Spin Hall [24], and Quantum Valley Hall effects [25]. These topologically protected modes emerge from symmetry-breaking in correspondence of Dirac cones and have been employed for applications such as waveguides with minimal backscattering through sharp corners or with resilience against defects [26], or for rainbow trapping of elastic waves [27,28].

In elasticity, the study of mass-spring lattices can provide fundamental insights into the dynamic behaviour of both discrete and continuous mechanical systems. Due to the simplicity of their building blocks, mass-spring lattices represent an ideal playground to explore

<footnote><sup>a</sup> Laboratory for Bio-inspired, Bionic, Nano, Meta Materials & Mechanics, Department of Civil, Environmental and Mechanical Engineering, University of Trento, 38123 Trento, Italy</footnote>

<footnote><sup>b</sup> Department of Applied Science and Technology, Politecnico di Torino, 10129 Torino, Italy</footnote>

<footnote><sup>c</sup> Univ. Lille, CNRS, Centrale Lille, Junia, Univ. Polytechnique Hauts-de-France, UMR 8520 - IEMN - Institut d' Electronique de Microelectronique et de Nanotechnologie, F-59000 Lille, France</footnote>

<footnote><sup>d</sup> School of Engineering and Materials Science, Queen Mary University of London, Mile End Road, London E1 4NS, United Kingdom</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>** Corresponding author at: Laboratory for Bio-inspired, Bionic, Nano, Meta Materials & Mechanics, Department of Civil, Environmental and Mechanical Engineering, University of Trento, 38123 Trento, Italy.</footnote>

<footnote>E-mail addresses: federica.ongaro@unitn.it (F. Ongaro), nicola.pugno@unitn.it (N.M. Pugno).</footnote>
the unconventional dynamic properties of ideal PCs or MMs [29]. Based on the correspondence between continuous systems and their discrete counterparts, a bulk of literature has been devoted to examining a variety of dynamic behaviours of one-, two-, and three-dimensional lattices, including band gaps in mono-atomic and diatomic chains [30], damping [31], negative refraction [32], waveguiding [33], cloaking [34], nonlinearity [35,36] and topological effects [37-43], to cite a few. Regarding the latter, mass-spring systems with a hexagonal configuration constitute a model system for invoking topological phases. For example, based on this model, the Quantum Valley Hall effects have been realised in [44] by alternating the masses of the unit cell, while the Quantum Spin Hall effects have been explored in [45] in mechanical granular graphene systems. Starting from a hexagonal mass-spring lattice, Quantum Spin Hall effects have also been analysed in [46] by varying the masses and the relative spring constants. This strategy, in particular, has been used to create, tune and invert the band structure, while preserving time reversal symmetry and thus emulating the Quantum Spin Hall effect. However, in most of these works, the interaction between masses or unit cells is limited to the nearest neighbours, i.e., only interactions between adjacent cells are considered to derive the dispersion relations.

On the contrary, recent studies have shown that considering interactions beyond the nearest neighbours leads to richer mechanics, both from the static and dynamic point of view.

In dynamics, non-nearest neighbour interactions can lead to the existence of multiple localised modes with a topological origin, identified by the number of Dirac points in the Brillouin zone or going beyond the conventional winding number in topological mechanics [37,47,48]. In [49], third-nearest neighbouring interactions were theoretically and numerically shown to induce topological phase transitions in one-dimensional mass-spring chains, while in [50,51], active methods to manipulate topological states via non-local connections were proposed.

Moreover, beyond nearest neighbour interactions lead to roton-like dispersion relations [52-54]. Originally predicted by Landau [55] to explain the unusual thermodynamic properties of superfluid helium, rotons are quasiparticles which can be observed in a quantum condensate with long-range dipolar interactions or spin-orbit coupling [56]. Their dispersion curves display both increasing and decreasing trends, included between a maximum (the maxon) and a minimum (the roton).

In the context of PCs and MMs, the slope inversion of the dispersion curves due to non-nearest neighbour interaction has been experimentally demonstrated in [57], where an acoustic MM including cube resonators and different types of connecting tubes mimicking non-neighbouring connections was manufactured: when increasing the order on long-range interactions, a larger number of slope inversions can be obtained [57]. These findings, in accordance with the works in [56,58], were also numerically and experimentally confirmed in [59] for the case of a fully reconfigurable elastic MMs built with the popular construction kit Meccano. Starting with a simple one-dimensional mass-spring model, non-neighbouring effects were also triggered in [60,61] by connecting the masses to their immediate neighbours and to an increasing order of non-nearest ones: again, dispersion relations characterised by an increased number of frequency regions with alternating positive and negative slopes were found. In [60], the same concept was extended to three-dimensional MMs, where the masses and springs of the one-dimensional mass-spring chain were replaced, respectively, by small cubes connected by cylindrical rods with a variable radius: a three-dimensional tube-based metamaterial for airborne sound is then proposed, allowing the authors to explore the effects of adding non-neighbouring interactions in the context of acoustical wave dispersion. More recently, Kazemi et al. [58] proposed a detailed method based on Fourier series coefficients for tailoring the dispersion relations of PCs with non-neighbouring connections. Similarly, in [62], beyond-nearest-neighbour interactions are used to tailor the properties of higher phonon bands, and obtain topological bandgaps. The effects of adding nonlinearity in conjunction with

non-neighbouring connections is explored in [63]: by focusing on a monoatomic mass-spring chain with up to third-neighbour interactions, a roton-like configuration emerges in the dispersion curves, together with an enhancement of the tunability provided by non-locality. Additional contributions were discussed in [64,65], where authors considered second-neighbouring interactions in addition to the nearest-neighbour ones. Farzbod et al. [64], in particular, established a force analysis for one-dimensional mass-spring chains, developing a theoretical framework for vibrational analysis of structures in which interactions are not limited to the nearest neighbours. An application of the developed method to a square mass-spring lattice was also illustrated, with special attention to structures with electromagnetic and/or electrostatic forces. By focusing on a square mass-spring lattice with equal masses, Wang et al. [65] observed negative refraction and directional radiation of elastic waves due to the non-neighbouring connection of masses. Similarly, in [60,66], adding beyond-the-nearest neighbour interactions to the classical first-nearest ones leads to unusual acoustic wave dispersion relations, mimicking roton-like behaviour. Finally, non-neighbouring interactions up to the third order in two-dimensional mass-spring hexagonal lattices were numerically investigated in [67], where the topological modes of a two-dimensional mass-spring hexagonal lattice with connections between both nearest and third nearest neighbouring masses was examined.

From this overview, it appears that most of the works currently available in literature are based on numerical or experimental investigations or, when an analytical approach is envisaged, it is often limited to one-dimensional mass-spring systems. To bridge this research gap, this paper investigates both analytically and numerically how elastic wave propagation in two-dimensional hexagonal diatomic mass-spring systems is affected by non-neighbouring interactions up to an arbitrary order, both for in-plane and out-of-plane displacements, and examines how this can be exploited for topologically protected waveguiding. Here, it is shown that non-nearest connections allow for (i) the formation of additional Dirac cones and (ii) a migration in their location in the reciprocal space as a function of the relative stiffness between nearest and third nearest neighbouring connections. These additional Dirac cones are linked to a corresponding increase in the number of topological edge modes, which hybridise and result in bound modes at interfaces between lattices that are inverted copies of each other.

## 2. Theoretical model: lattice description and dispersion relations

In this section, we briefly describe the main steps of Bloch's analysis allowing us to derive a closed-form expression for the dispersion equation of a 2D mass-spring system. A "local" (first-neighbouring) system is initially considered in Section 2.1. The theory is then extended to the case of not-neighbouring connections in Section 2.2, where a scenario with a general number of nth connections is proposed.

### 2.1. First nearest connections (local interactions)

Let us focus on the infinite two-dimensional hexagonal diatom lattice, as the one illustrated in Fig. 1(a). The unit cell, identified by the light blue parallelogram, is made of two different masses, $m_a$ and $m_b$, periodically arranged in a hexagonal configuration of lattice constant $\ell$. The periodicity and infinite extent of the structure allow us to restrict the analysis to the basic unit cell (0,0) tessellating the entire structure when repeated along the lattice direct vectors $\mathbf{d}_1 = \ell e_1$ and $\mathbf{d}_2 = (\ell/2)e_1 + (\ell/\sqrt{3}/2)e_2$. With reference to the chosen base cell, any other cell obtained by $n_1$ integer translations along $\mathbf{d}_1$ and $n_2$ integer translations along $\mathbf{d}_2$, is defined by the two coordinates $(n_1, n_2)$. Similarly, denoting with $\mathbf{r}_j^{(0,0)}$ the position vector of the general jth physical lattice point in the base cell (0,0), the corresponding point in the $(n_1, n_2)$ cell is identified by the position vector

$$ \mathbf{r}_j^{(n_1, n_2)} = \mathbf{r}_j^{(0,0)} + n_1 \mathbf{d}_1 + n_2 \mathbf{d}_2. \quad (1) $$
<figure><img src="image_4.png" /><figcaption>Fig. 1. The non-local mass-spring system: (a) the two-dimensional hexagonal configuration described by the lattice vectors (d₁, d₂) and reference unit cell (highlighted) composed by the pair of masses m_a and m_b; (b), (c) non-local interactions up to the third order for mass (b) m_a and (c) m_b. In both scenarios, elastic springs are represented by continuous lines while the coloured parallelograms denote the unit cells in the case of first- (blue), second- (yellow) and third- (magenta) level connections.</figcaption></figure>

Finally, each mass of the lattice is connected to its immediate
neighbouring mass by a linear elastic spring of stiffness k₁ indicated,
for simplicity, by the blue continuous lines in Figs. 1(b), 1(c). Note
that, for sake of clarity, connections are only illustrated for the unit
cell (0,0). Finally, each mass is assumed to have only one degree of
freedom, corresponding to an out-of-plane displacement.

Imagine now to apply a time-harmonic vibration to the considered
system. Indicating with $u_a^{(0,0)}$ and $u_b^{(0,0)}$ the displacements of the $m_a$ and
$m_b$ masses contained in the cell (0,0), the governing equations of motion
for infinitesimal deformations read

$$
\left\{
\begin{aligned}
m_a \ddot{u}_a^{(0,0)} &= k_1 (\dot{u}_b^{(0,0)} - u_a^{(0,0)}) + k_1 (\dot{u}_b^{(-1,0)} - u_a^{(0,0)}) + k_1 (\dot{u}_b^{(0,-1)} - u_a^{(0,0)}) \\
m_b \ddot{u}_b^{(0,0)} &= k_1 (\dot{u}_a^{(0,0)} - u_b^{(0,0)}) + k_1 (\dot{u}_a^{(1,0)} - u_b^{(0,0)}) + k_1 (\dot{u}_a^{(0,1)} - u_b^{(0,0)}),
\end{aligned}
\right.
\qquad (2)
$$

with

$$
\ddot{u}_{(-)}^{(0,0)} := \frac{\partial^2 u_{(-)}^{(0,0)}}{\partial t^2} \qquad (3)
$$

and $u_i^{(n_1, n_2)}$ the displacement of the $i = (a, b)$ mass pertaining to the unit cell with coordinates $(n_1, n_2)$. For example, $u_b^{(-1,0)}$ stands for the displacement of the $m_b$ mass contained in the cell $(-1,0)$, which is connected to the mass $m_a$ in $(0,0)$ according to the first neighbouring connection scheme (Figs. 1(b), 1(c)).

If a plane wave solution is assumed, the generalised displacement
of the two masses in the reference cell (0,0) can be written in terms of

their position and time dependence as

$$
u_j^{(0,0)} \equiv u_j^{(0,0)}(\mathbf{r}_j^{(0,0)}, t) = \hat{u}_j^{(0,0)}(\mathbf{r}_j^{(0,0)}) e^{-i\omega t} = \hat{u}_j^{(0,0)} e^{i\mathbf{q}\cdot\mathbf{r}_j^{(0,0)} - \omega t}, \quad j = a, b, \qquad (4)
$$

with $i := \sqrt{-1}$ the imaginary unit, $\omega$ the circular frequency of the time-
harmonic wave propagating along the system, $\hat{u}_j^{(0,0)}$ the wave amplitude
and $\mathbf{q} := q_x \mathbf{e}_1 + q_y \mathbf{e}_2$ the wave vector that defines the direction of
propagation. Also, given the periodicity of the examined lattice, the
displacement $u_j^{(n_1, n_2)}$ can be obtained by applying Bloch's theorem [68]

$$
\begin{align}
u_j^{(n_1, n_2)} &\equiv u_j^{(n_1, n_2)}(\mathbf{r}_j^{(n_1, n_2)}, t) = u_j^{(0,0)}(\mathbf{r}_j^{(0,0)}, t) e^{i\mathbf{q} \cdot (\mathbf{r}_j^{(n_1, n_2)} - \mathbf{r}_j^{(0,0)})} \nonumber \\
&= \hat{u}_j^{(0,0)}(\mathbf{r}_j^{(0,0)}) e^{i(\mathbf{q} \cdot (n_1 \mathbf{d}_1 + n_2 \mathbf{d}_2) - \omega t)}, \quad j = a, b. \tag{5}
\end{align}
$$

Substituting Eqs. (4), (5) into Eq. (2) and omitting, for simplic-
ity, the dependence on the time-harmonic factor $e^{-i\omega t}$, provides the
homogeneous system

$$
\begin{equation}
\begin{aligned}
& \begin{cases}
-m_a \omega^2 \hat{u}_a = \hat{u}_a (-3k_1) + \hat{u}_b (1 + e^{-i\mathbf{q} \cdot \mathbf{d}_1} + e^{-i\mathbf{q} \cdot \mathbf{d}_2}) k_1 \\
-m_b \omega^2 \hat{u}_b = \hat{u}_a (1 + e^{-i\mathbf{q} \cdot \mathbf{d}_1} + e^{i\mathbf{q} \cdot \mathbf{d}_2}) k_1 + \hat{u}_b (-3k_1),
\end{cases}
\end{aligned}
\tag{6}
\end{equation}
$$

where the notation $\hat{u}_{(-)} := \hat{u}_{(-)}^{(0,0)}(\mathbf{r}_{\mathbf{j}}^{(0,0)})$ is used to improve readability.
Expressed in matrix form, Eq. (6) leads to the classical eigenproblem

$$
(k_1 \mathbf{K}_1(\mathbf{q}) - \omega^2 \mathbf{M}) \mathbf{U} = \mathbf{0}, \qquad (7)
$$
with

$$
\mathbf{K}_1(\mathbf{q}) = \begin{bmatrix} 3 & -(1+e^{-i\mathbf{q}\cdot\mathbf{d}_1}+e^{-i\mathbf{q}\cdot\mathbf{d}_2}) \\ -(1+e^{i\mathbf{q}\cdot\mathbf{d}_1}+e^{i\mathbf{q}\cdot\mathbf{d}_2}) & 3 \end{bmatrix}, \\
\mathbf{M} = \begin{bmatrix} m_a & 0 \\ 0 & m_b \end{bmatrix}, \quad \mathbf{U} = \begin{bmatrix} \hat{u}_a \\ \hat{u}_b \end{bmatrix}. \tag{8}
$$

The non-trivial solution of Eq. (7), obtained by imposing the determinant of the matrix of coefficients to vanish, i.e.,

$$
\det (k_1 \mathbf{K}_1(\mathbf{q}) - \omega^2 \mathbf{M}) = 0, \quad (9)
$$

provides the dispersion relation $\omega(\mathbf{q})$ of the investigated lattice system:

$$
\omega^4 m_a m_b - 3\omega^2 k_1 (m_a + m_b) + k_1^2 (9 - \Psi'_1(\mathbf{q})) = 0, \quad (10)
$$

from which

$$
\omega^2(\mathbf{q}) = k_1 \frac{3(m_a + m_b) \pm \sqrt{9(m_a + m_b)^2 - 4m_a m_b (9 - \Psi'_1(\mathbf{q}))}}{2m_a m_b}. \quad (11)
$$

Here, for conciseness,

$$
\Psi'_1(\mathbf{q}) := 3 + e^{i\mathbf{q}\cdot\mathbf{d}_1} + e^{-i\mathbf{q}\cdot\mathbf{d}_1} + e^{i\mathbf{q}\cdot\mathbf{d}_2} + e^{-i\mathbf{q}\cdot\mathbf{d}_2} + e^{i\mathbf{q}\cdot(\mathbf{d}_1 - \mathbf{d}_2)} + e^{-i\mathbf{q}\cdot(\mathbf{d}_1 - \mathbf{d}_2)} \quad (12)
$$

or, expressing the exponential terms in trigonometric functions by the
Euler's equation

$$
e^{\pm i\alpha} = \cos \alpha \pm i \sin \alpha \quad \text{and} \quad e^{i\alpha} + e^{-i\alpha} = 2 \cos \alpha, \tag{13}
$$

$$
\Psi'_1(\mathbf{q}) = 3 + 2\cos(q_x \ell) + 2\cos\left(q_x \frac{\ell}{2} + q_y \frac{\sqrt{3}}{2} \ell\right) + 2\cos\left(q_x \frac{\ell}{2} - q_y \frac{\sqrt{3}}{2} \ell\right). \quad (14)
$$

Clearly, when taking the square root on both sides of Eq. (11), we
obtain two signs for $\omega(\mathbf{q})$. As usual, we adopt the convention to consider
positive angular frequencies.

Before concluding, it is important to highlight the difference existing between the position vector and the displacement field. The position vector, which in our analysis is denoted with $\mathbf{r}_j^{(n_1,n_2)}$, defines the position of a given mass pertaining to the cell $(n_1,n_2)$ as a function of the vector $\mathbf{r}_j^{(0,0)}$, defining the position of the corresponding mass in the reference cell (0,0). $\mathbf{r}_j^{(\cdot,\cdot)}$, in particular, is a plane vector expressed in terms of the lattice vectors $\mathbf{d}_1$ and $\mathbf{d}_2$. According the Bloch theory, this allows us to analyse the entire lattice by focusing on a single unit cell. We have then introduced the displacement field $u_j^{(\cdot,\cdot)}$, corresponding to the out-of-plane displacement of the mass located at the position $\mathbf{r}_j^{(\cdot,\cdot)}$, and wrote the equilibrium equations as a function of it. $\mathbf{r}_j^{(\cdot,\cdot)}$ and $u_j^{(\cdot,\cdot)}$ are thus very different quantities: the first defines the position of the masses in the lattice plane, the second the displacement field (out-of-plane in this case) associated to them.

2.2. Beyond nearest-neighbour connections (non-local interactions): the
general case of nth nearest connections

Consider again the hexagonal lattice in Fig. 1 and imagine to con-
nect each $m_a$ and $m_b$ mass to its $n$th nearest masses, with $n = 1, 2, 3, ...$ ($n = 1$ coincides with the nearest interactions in Section 1). Connections are realised by linear springs of stiffness $k_n$.

The reasoning of Section 2.1 can be clearly extended to this general
case. In particular, writing the equilibrium equations for the $m_a$ and
$m_b$ masses in the base cell (0,0) and solving the corresponding eigen-
problem, provides the dispersion relation for a time-harmonic vibration
propagating within the system. This leads to

$$
\omega^2(\mathbf{q}) = k_n \frac{3(m_a + m_b) \pm \sqrt{9(m_a + m_b)^2 - 4m_a m_b (9 - \Psi'_n(\mathbf{q}))}}{2m_a m_b}, \quad (15)
$$

with

$$
\Psi_n(\mathbf{q}) \equiv \Psi'_1(2^{n-1}\mathbf{q}), \quad n = 1, 2, \dots \tag{16}
$$

such that $\Psi'_n(\mathbf{q}) \equiv \Psi_{n-1}(2\mathbf{q})$.

It is interesting to note that the dispersion relations obtained in the case of first- (Eq. (11)) and nth (Eq. (15)) neighbouring interactions have the same form. The only difference is the replacement of the quantities $k_1$ and $\Psi'_1(\mathbf{q})$ with $k_n$ and $\Psi'_n(\mathbf{q})$. This determines the possibility of establishing a rule of thumb to easily derive the dispersion equation, and corresponding properties like bandgaps, in the general case of nth nearest connections by just solving the simplest case of first nearest interactions. For the interested reader, this is demonstrated in Appendix A, where the extended calculations for deriving the dispersion relation for the cases of second- ($n = 2$) and third- ($n = 3$) nearest connections are reported.

An important comment concerns the configuration of the non-local interactions illustrated in Fig. 1. We know that, when non-local interactions are assumed periodic, two-dimensional regular lattices with a given topology (hexagonal in the examined case), are characterised by a non-univocal choice in passing from local to non-local configurations. That is, a prescribed level of non-locality does not imply one possible configuration of the periodic non-local interactions. In the examined case, in particular, non-local connections copy the hexagonal pattern of the fundamental lattice, whereas both neighbouring and non-neighbouring masses could be connected without following the hexagonal symmetry as well. Even if sharing the same level of non-locality, any differently arranged non-local interactions mounted on the same fundamental hexagonal lattice would be associated to completely different dynamics and dispersion relations. However, the family of generated non-local lattices could be not so trivial to describe from a dynamical point of view, even with the application of the Bloch's theory. It is thus important to highlight that the dynamic equations presented in this and the following sections are referred to the particular choice of tracing a hexagonal pattern for the non-local connections.

2.3. Local plus non-local interactions: simultaneous nearest and nth beyond nearest neighbour connections

The model presented in the previous sections can be easily ex-
tended to contain not only nearest (or n-th nearest) connections but
also multiple orders of beyond-nearest connections in addition to the
nearest ones. In particular, let us consider the case of nearest-neighbour
connections (n = 1) together with n-orders of beyond-nearest-neighbour
interactions. In this configuration, obtained by “superimposing” the
configurations described in Sections 2.1 and 2.2, the m_a and m_b masses
interact with nearest neighbours via linear springs of stiffness k_1 and
with the n-th non-nearest neighbours via linear springs of stiffness k_n.

As usual, the dispersion relation for a time-harmonic vibration prop-
agating within the lattice can be obtained by writing the equilibrium
equations for the m_a and m_b masses in the (0,0) cell and solving the
eigenproblem that emerges

$$
\left( \sum_{i=1}^{n} k_i \mathbf{K}_j(\mathbf{q}) - \omega^2 \mathbf{M} \right) \mathbf{U} = \mathbf{0}, \quad (17)
$$

with **M** and **U** derived from Eq. (8) and **K**<sub>i</sub>(**q**) the stiffness matrix corresponding to the *i*th neighbouring connections (cf. Eq. (8) and Appendix A). Finally, imposing

$$
\det \left( \sum_{i=1}^{n} k_i \mathbf{K}_j(\mathbf{q}) - \omega^2 \mathbf{M} \right) = 0, \quad (18)
$$

provides

$$
\omega^4 m_a m_b - 3\omega^2 \sum_{i=1}^{n} k_i (m_a + m_b) + \sum_{i=1}^{n} k_i^2 (9 - \Psi'_i(\mathbf{q})) \\
+ 2 \sum_{i=1}^{n-1} \left( \sum_{j=i+1}^{n} k_i k_j (9 - \Psi'_{i,j}(\mathbf{q})) \right) = 0, \quad (19)
$$
where $\Psi_i(q) \equiv \Psi'_i(2^{n-1}q)$, $i = 1, 2, \dots, n$ (cf. Eq. (16)),

$$ \Psi'_{ij}(q) = \frac{K_i(2, 1)K_j(1, 2) + K_j(2, 1)K_i(1, 2)}{2}, \quad i, j = 1, 2, \dots, n \quad (20) $$

and $K_e(m, n)$ the $(m, n)$th component of the stiffness matrix describing the eth order of neighbouring connections. The latter are obtained by considering the $m_a$ and $m_b$ masses connected with their eth neighbours only. This provides a quicker way to obtain the coupling terms $\Psi_{ij}(q)$ in Eq. (19).

As can be expected for linear elastic systems, from Eq. (19) it also emerges that the examined configuration can be seen as a superimposition of multiple nth neighbouring connections (second and third terms) plus their mutual interactions (last term), expressed as a function of the corresponding spring stiffness.

More details about the derivation of the listed equations can be found in Appendix B, where the calculations for the case of first- ($n=1$) plus second- ($n=2$) and third- ($n=3$) neighbouring connections are reported.

### 3. Generalising the theory to the case of in-plane displacements

We now wish to extend the analysis presented in the previous sections, based on the assumption of transverse displacements of the masses, to the case of in-plane mass motion. In particular, the equilibrium equations and dispersion relations that we obtain in the two situations have exactly the same form, thus revealing the equivalence existing between the two cases in terms of wave dispersion characteristics.

This can be demonstrated by focusing on Fig. 2(a), illustrating the simplest case of two masses, $m_a$ and $m_b$, connected by an elastic spring of stiffness $k$.

Denoting with $\Delta x$ and $\Delta y$ the elongation of the spring in the $e_1$ and $e_2$ direction and with $u_a = u_{x,a}e_1 + u_{y,a}e_2$ and $u_b = u_{x,b}e_1 + u_{y,b}e_2$ the in-plane displacements of the two masses, the elastic force $F_a$ acting on the mass $m_a$ is given by

$$ |\mathbf{F}_a| = k \sqrt{\Delta x^2 + \Delta y^2} = k \sqrt{(u_{x,b} - u_{x,a})^2 + (u_{y,b} - u_{y,a})^2}. \quad (21) $$

Its components along the $x$ and $y$ directions are, respectively,

$$ F_{a,x} = |\mathbf{F}_a|\cos\theta, \quad F_{a,y} = |\mathbf{F}_a|\sin\theta, \quad (22) $$

being $\theta$, the spring inclination angle illustrated in Fig. 2(a), such that

$$ \sin\theta = \frac{\Delta y}{\sqrt{\Delta x^2 + \Delta y^2}}, \quad \cos\theta = \frac{\Delta x}{\sqrt{\Delta x^2 + \Delta y^2}}. \quad (23) $$

Substituting Eq. (23) into Eq. (22) leads to

$$ F_{a,x} = k\Delta x = k(u_{x,b} - u_{x,a}), \quad F_{a,y} = k\Delta y = k(u_{y,b} - u_{y,a}). \quad (24) $$

Finally, writing the equilibrium equations for the mass $m_a$ in the $x$ and $y$ directions, provides

$$ \left\{ \begin{array}{l} m_a \ddot{u}_{x,a} = k (u_{x,b} - u_{x,a}) \\ m_a \ddot{u}_{y,a} = k (u_{y,b} - u_{y,a}). \end{array} \right. \quad (25) $$

Consider now Fig. 2(b), which more closely resembles the theoretical model described in Section 2: a central mass $m_a$ connected by elastic springs of stiffness $k$ to the three $m_b$ masses $m_{b1}$, $m_{b2}$, $m_{b3}$. Again, it emerges that the $x$ and $y$ components of the elastic forces $F_{a1}$, $F_{a2}$, $F_{a3}$ acting on $m_a$ due to the springs elongation, take the form

$$ F_{a1,x} = k(u_{x,bi} - u_{x,a}), \quad F_{a1,y} = k(u_{y,bi} - u_{y,a}), \quad i = 1, 2, 3, \quad (26) $$

from which the following equilibrium equations can be established

$$ \left\{ \begin{array}{l} m_a \ddot{u}_{x,a} = k(u_{x,b1} - u_{x,a}) + k(u_{x,b2} - u_{x,a}) + k(u_{x,b3} - u_{x,a}) \\ m_a \ddot{u}_{y,a} = k(u_{y,b1} - u_{y,a}) + k(u_{y,b2} - u_{y,a}) + k(u_{y,b3} - u_{y,a}). \end{array} \right. \quad (27) $$

Let us now apply these general concepts to the hexagonal lattice presented in Section 2 and, in particular, let us focus on the two masses $m_a$ and $m_b$ pertaining to the reference cell (0,0) highlighted in Fig. 1. Even if, for simplicity, the case of first-neighbouring connections

is considered in the following (Figs. 1(b), 1(c)), analogous considerations apply for the nth neighbouring ones.

By adopting the same notation of Section 2, writing the equilibrium equations for the $m_a$ and $m_b$ masses in the $x$ and $y$ directions provides

$$ \left\{ \begin{aligned} m_a \ddot{u}_{j,a}^{(0,0)} &= k_1 (u_{j,b}^{(0,0)} - u_{j,a}^{(0,0)}) + k_1 (u_{j,b}^{(-1,0)} - u_{j,a}^{(0,0)}) + k_1 (u_{j,b}^{(0,-1)} - u_{j,a}^{(0,0)}) \\ m_b \ddot{u}_{j,b}^{(0,0)} &= k_1 (u_{j,a}^{(0,0)} - u_{j,b}^{(0,0)}) + k_1 (u_{j,a}^{(1,0)} - u_{j,b}^{(0,0)}) + k_1 (u_{j,a}^{(0,1)} - u_{j,b}^{(0,0)}), \quad j = x, y, \end{aligned} \right. \quad (28) $$

with $u_{x,i}^{(n_1, n_2)}$, $u_{y,i}^{(n_1, n_2)}$ the in-plane displacements of the $i = (a, b)$ mass pertaining to the $(n_1, n_2)$ cell. Again, if a plane wave is assumed, applying the Bloch's theorem (cf. Eq. (5)) leads to the homogeneous system

$$ \left\{ \begin{aligned} -m_a \omega^2 \ddot{u}_{j,a} &= \hat{u}_{j,a} (-3k_1) + \hat{u}_{j,b} (1 + e^{-i\mathbf{q}\cdot\mathbf{d}_1} + e^{-i\mathbf{q}\cdot\mathbf{d}_2}) k_1 \\ -m_b \omega^2 \ddot{u}_{j,b} &= \hat{u}_{j,a} (1 + e^{-i\mathbf{q}\cdot\mathbf{d}_1} + e^{i\mathbf{q}\cdot\mathbf{d}_2}) k_1 + \hat{u}_b (-3k_1), \quad j = x, y, \end{aligned} \right. \quad (29) $$

where $\hat{u}_{(c)} := \hat{u}_{(c)}^{(0,0)}(\mathbf{r}_j^{(0,0)})$ (cf. Section 2.1).

Expressed in matrix form, Eq. (29) yields

$$ (k_1 \mathbf{K}_1(\mathbf{q}) - \omega^2 \mathbf{M}) \mathbf{U}_j = \mathbf{0}, \quad j = x, y, \quad (30) $$

with $\mathbf{K}_1(\mathbf{q})$ and $\mathbf{M}$ the matrices defined in Eq. (8) and

$$ \mathbf{U}_j = \begin{bmatrix} \hat{u}_{j,a} \\ \hat{u}_{j,b} \end{bmatrix}, \quad j = x, y. \quad (31) $$

Finally, imposing the determinant of the matrix of coefficients to vanish, provides the corresponding dispersion relation. It emerges that the latter has exactly the same form of Eq. (11), given the equivalence between Eqs. (7) and (30). This result, revealing the equivalence existing between the in-plane and out-of-plane wave dispersion characteristics, is numerically confirmed in Section 5.

In the above calculations note that no particular values of the angle $\theta$ and/or lattice configuration have been assumed. The proposed method has thus general applicability.

### 4. Band diagrams

The dispersion equations obtained in the previous Sections allow the computation of the so-called band diagrams, illustrating the relation between wavevectors and wave propagation frequencies and highlighting negative refraction effects. Bandgap opening and closing frequencies can also be determined.

In this section, we investigate how the introduction of non-local effects influences the dispersion curves and frequency bandgaps of a traditional locally interacting mass-spring system. Different scenarios are considered, involving a hexagonal mass-spring system with first-neighbour connections having springs of stiffness $k_1$, second-neighbour connections with springs stiffness $k_2$, third-neighbour connections with springs stiffness $k_3$, first- plus second- plus third-neighbour connections with spring stiffnesses of respectively $k_1$, $k_2$, $k_3$ (Figs. 1(b), 1(c)). In all cases, $m_a = 0.001$ kg and $k_1 = 10^6$ N/m while, to understand how the system characteristics, i.e., mass and stiffness, affect the dispersion behaviour, the parameters $m_b$, $k_2$ and $k_3$ are taken as variable. Specifically, denoting with $\alpha := m_b/m_a$, $\beta := k_2/k_1$ and $\gamma := k_3/k_1$, the following values are considered: $\alpha = (0.5, 1, 1.5)$, $\beta = (0.1, 1, 10)$, $\gamma = (0.1, 1, 10)$.

It should be noted that different possibilities exist to define the springs stiffness of the various order connections. One solution, adopted in the present work, is to introduce two scaling parameters, $\beta$ and $\gamma$, and express the stiffness of non-local connections as a function of the stiffness of the local ones: $k_2 = \beta k_1$, $k_3 = \gamma k_1$. Changing the value of
<figure><img src="image_5.png" /><figcaption>Fig. 2. Theory for the case of in-plane displacements: (a) masses $m_a$ and $m_b$ connected by a linear elastic spring of stiffness $k$ experiencing the elongation $\Delta x$ ($e_1$ direction) and $\Delta y$ ($e_2$ direction) and corresponding elastic force $F_a$ acting on $m_a$, (b) central mass $m_a$ connected to the three masses $m_{b1}, m_{b2}, m_{b3}$ by linear elastic springs of stiffness $k$ and corresponding elastic forces $F_{a1}, F_{a2}, F_{a3}$.</figcaption></figure>

$\beta$ and $\gamma$ allows us to tailor the strength of the ith order connection. Clearly, this can also be done when higher order connections are introduced, like the 4th, 5th, ... A second possibility, probably more convenient to use when a larger number of non-local interactions is present, is to define the spring stiffness according to a law of the type $k_i = k_1 i^{u_s}$ or $k_i = k_1 / i^{u_s}$, with $u_s$ a scaling parameter.

Given the periodicity of the examined lattice, the usual procedure for the analysis of wave propagation consists in defining the wavevector $q$ constrained to the first Brillouin zone (FBZ) [68], illustrated in Fig. 3(d) together with the components of the wavevector $q_x, q_y$. However, once symmetry is maintained, to determine the band diagrams, it is sufficient only to scan the contour of the irreducible Brillouin zone (IBZ) instead of considering every wavevector inside the FBZ [69]. If a bandgap is detected on the contour of the IBZ, it is also present in its interior and, by symmetry, in the whole FBZ. This can be seen in Fig. 3(d) where the IBZ, defined by the points $K = (2\pi/3\ell, 2\pi/\sqrt{3}\ell)$, $\Gamma = (0, 0)$ and $M = (0, 2\pi/\sqrt{3}\ell)$ is highlighted.

The study can be clearly extended to n-orders of beyond-nearest neighbour connections (cf. Sections 3, 4). Here, for simplicity, non-local interactions only up to $n=3$ are considered since, as shown below, the resulting behaviour is already rich and complex.

## 4.1. Nearest and beyond-nearest neighbour connections: analysis of the influencing parameters on the bandgap formation

Let us consider the case of nearest-neighbour connections ($n=1$). The corresponding dispersion diagrams are illustrated in Fig. 3 as a function of the non-dimensional frequency $\hat{\omega} := \omega / \sqrt{k_1 / m_a}$ and wavenumber $\hat{q} := q \ell / \pi$.

Different values of the ratio $\alpha$ between masses $m_a$ and $m_b$, are considered: $\alpha = 0.5$ in Fig. 3a, $\alpha = 1$ in Fig. 3b, $\alpha = 1.5$ in Fig. 3(c). As can be seen, for $\alpha = 1$, i.e., when the two masses $m_a$ and $m_b$ are equal, the investigated configuration exhibits no complete bandgaps since, at any given frequency, there is always at least one propagating wave. Conversely, for $\alpha = 0.5$ and $\alpha = 1.5$ one complete bandgap appears, spanning the range $\hat{\omega} = [1.73, 2.45]$ and $\hat{\omega} = [1.41, 1.73]$, respectively in the two cases. Similar considerations apply for second- and third-neighbouring connections, as reported in Fig. 3. Again, there are no complete bandgaps for $\alpha = 1$ and one complete bandgap for $\alpha = 0.5$ and $\alpha = 1.5$, located, respectively, at $\hat{\omega} = [1.73, 2.45]$ and $\hat{\omega} = [1.41, 1.73]$ for second-neighbouring connections and at $\hat{\omega} = [1.73, 2.45]$ and $\hat{\omega} = [1.41, 1.73]$ for the third-neighbouring ones.

Considering non-local effects provides a more complex behaviour of the dispersion curves, which display multiple slope inversions inside the $\Gamma - M, M - K$ and $K - \Gamma$ boundaries. More and more local minima are also observed by increasing the order of non-locality $n$. For second-order connections, the lowest dispersion curve goes to zero

in the symmetry point as well as in $\Gamma$. For third-order connections, other zeros occur between $\Gamma - M$ and $K - \Gamma$. Conversely, the Dirac cone in $K$ is always conserved, independently of the order of non-nearest interactions. Interestingly, by tuning the order of non-local effects, the position of the negative-slope regions can be controlled, so that they occur at a desired frequency and wavenumber. This is particularly evident by comparing the yellow and green curves in Fig. 3, corresponding, in turn, to the second- and third-neighbour connections. Note that this slope inversion does not occur without beyond-nearest neighbour interactions, which are enabled by the addition of the $k_2$ and $k_3$ springs. Regarding this aspect, varying the parameters $\beta$ and $\gamma$, i.e., the stiffness of the $k_2$ and $k_3$ springs, for a fixed value of $\alpha = 0.5$, i.e., $m_b = 0.5 m_a$, dramatically changes the frequency width and position of the dispersion curves, as reported in Fig. 4. For example, in the case of second-neighbour connections, increasing $\beta$ from 0.1 (Fig. 4(a)) to 1 (Fig. 4(b)) and 10 (Fig. 4(c)) leads to an increase of the bandgap width of approximately 300% and 1000%. A shift to higher frequencies also emerges. A comparison with the dispersion curves of first-nearest connections is also established in Fig. 4, from which it emerges that for $\beta = 1$ (Fig. 4(b)) and $\gamma = 1$ (Fig. 4(e)) both the bandgap amplitude and position are exactly the same, regardless of the order of non-locality $n$.

## 4.2. Local plus non-local interactions: effects of changing parameters on the bandgap formation

To further characterise the dynamic response of the proposed mass-spring system, we extend the study in Section 4.1 to the combined action of local plus non-local interactions.

In the investigated scenario, the $m_a$ and $m_b$ masses, disposed in a hexagonal configuration, are connected to their first-, second- and third-neighbours by springs having stiffness $k_1, k_2$ and $k_3$, respectively. As in Section 4.1, $m_a = 0.001$ kg and $k_1 = 10^6$ N/m while, to explore the effects of changing the mass and stiffness characteristics of the system, the parameters $\alpha := m_b / m_a$, $\beta := k_2 / k_1$ and $\gamma := k_3 / k_1$ are set to vary so that $\alpha = (0.5, 1, 1.5)$, $\beta = (0.1, 1, 10)$, $\gamma = (0.1, 1, 10)$.

Results are presented in Figs. 5 and 6 where, in accordance with Section 4.1, non-dimensional frequency $\hat{\omega} := \omega / \sqrt{k_1 / m_a}$ and wavenumber $\hat{q} := q \ell / \pi$ are used.

Let us initially focus on Fig. 5, where the dispersion curves corresponding to first-, second-, third- and first- plus second- plus third-neighbouring interactions are compared for a fixed value of $\alpha = 0.5$ and varying values of the parameters $\beta$ and $\gamma$. In general, a complex configuration of the dispersion curves emerges that, similarly to Section 4.1, display multiple slope inversions inside the $\Gamma - M, M - K$ and $K - \Gamma$ boundaries. However, some peculiarities can be observed. The first can be seen in Fig. 5(a), where $\beta = \gamma = 0.1$, i.e., springs $k_1$ stiff, springs $k_2$ and $k_3$ soft. The dispersion curves corresponding to the first- plus second- plus third-order connections (the red curve)
<figure><img src="image_6.png" /><figcaption>Fig. 3. Dispersion curves corresponding to first-, second- and third-neighbouring interactions for varying values of α. Comparison of the obtained dispersion curves in the case of $m_a = 0.001$ kg, $k_1 = 10^6$ N/m, $\beta = \gamma = 1$ and: (a) $\alpha = 0.5$, (b) $\alpha = 1$, (c) $\alpha = 1.5$. (d) The irreducible Brillouin zone for a hexagonal lattice highlighting the high-symmetry points K, M, Γ.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Dispersion curves obtained in the case of local and non-local effects for varying values of $\beta$ and $\gamma$. Comparison of the dispersion curves corresponding to $m_a = 0.001$ kg, $k_1 = 10^6$ N/m, $\alpha = 0.5$ and: (a) $\beta = 0.1$, (b) $\beta = 1$, (c) $\beta = 10$; (d) $\gamma = 0.1$, (e) $\gamma = 1$, (f) $\gamma = 10$.</figcaption></figure>
closely resemble those corresponding to the first-order ones (the blue line). Conversely, when $\beta = 0.1$ and $\gamma = 10$, i.e., spring $k_2$ is soft and spring $k_3$ is stiff (Fig. 5(c)), the dispersion curves of the first-plus second-plus third-neighbouring connections resemble those of the third-neighbouring ones (green line). Finally, for $\beta = 0.1$ and $\gamma = 1$, a sort of 'mixed' behaviour can be observed, in which the configuration of the dispersion curves corresponding to the first-plus second-plus third-order interactions are a sort of superposition of those corresponding to the first-(blue curve), second-(yellow curve) and third-order (green curve) ones. That is, the red curves in Fig. 5(b) are more complex than those in Fig. 5(a), i.e., multiple slope inversions and local minima occur, but they are smoother than those of Fig. 5(c), with a higher number of slope inversions and an higher difference between local minima and maxima.

Clearly, the stiffer the springs $k_i$, the stronger the influence of the ith non-local effect will be.

These results are confirmed in Figs. 5(d-f), based on a fixed value of $\beta = 1$ and a variable value of $\gamma$: $\gamma = 0.1$ in Fig. 5(d), $\gamma = 1$ in Fig. 5(e), $\gamma = 10$ in Fig. 5(f). Again, a sort of mitigation effect can be seen for the dispersion curves of the first-plus second-plus third-order connections even if, in this case, the influence of the second-neighbour ones if higher than that in Figs. 5(a-c).

Consider now Fig. 6, where the dispersion curves corresponding to the first-plus second-plus third-order interactions are depicted for a fixed value of $\beta = 10$ and variable values of $\gamma$ and $\alpha$. Here, no complete bandgaps emerge for $\alpha = 1$, that is $m_a = m_b$ (Fig. 6(b)), regardless of the values of $k_2$ and $k_3$, i.e., the parameters $\beta$ and $\gamma$. Conversely, for $\alpha = 0.5$ (Fig. 6(a)) and $\alpha = 1.5$ (Fig. 6(c)), one complete bandgap appears. However, a decreasing bandgap width of approximately 55% can be identified for increasing values of $\alpha$, together with a shift to lower frequencies of $\tilde{\omega}$.

Similar results corresponding to the cases with $\beta = 0.1$ and $\beta = 1$ are reported in Appendix C, for which similar considerations still apply. The illustrated examples, in particular, confirm that, when $\alpha = 1$, adjusting the parameters $\beta$ and $\gamma$ does not open the bandgap. Specifically, no bandgaps are opened by changing the parameters $\beta$, from 0.1 to 1, and $\gamma$, taken in the range [0.1, 1, 10], when $\alpha = 1$. Conversely, a more significant influence of $\beta$ and $\gamma$ is observed in the bandgap width and position for $\alpha \neq 1$.

## 4.3. Discussion

The results presented above show that, in terms of dispersion properties, including non-local interactions in addition to local ones leads to roton-like dispersion curves displaying multiple slope inversions inside the first Brillouin zone. A shift in the bandgaps is also observed. From our analysis, it emerges that these effects are more and more significant, and thus not negligible, for large values of the parameter $n$, i.e., the non-locality order. In this case, in particular, an increasing number of slope inversions is observed, independently of the parameters $\alpha := m_b/m_a$, $\beta := k_2/k_1$ and $\gamma := k_3/k_1$. We also found that non-local effects are not negligible in terms of bandgap width and position when the stiffness of the non-local springs, $k_2$ and $k_3$, becomes larger, i.e., for increasing values of the parameters $\gamma$ and $\beta$. Conversely, in the particular case $\beta = \gamma = 1$, configurations with 2nd and 3rd order connections provide a bandgap width coinciding with that of the 1st order one. Regarding the parameter $\alpha$, we observed that in the limit situation of $\alpha = 1$, i.e., equal masses, no frequency bandgaps occur by adding non-local effects, independently of the strength of the interactions. For $\alpha \neq 1$, non-local effects are not negligible, especially for large values of the non-locality order $n$.

# 5. Numerical validation of the model

Finite Element simulations using Comsol Multiphysics 6.2 are now presented to verify the analytical relations previously obtained. Different scenarios are investigated in which the planar system is simulated in an eigenfrequency study in three dimensions. Bloch-Floquet boundary conditions are applied to the unit cells with first-, second- and third-neighbouring interactions, illustrated in Fig. 7. The mass-spring model is reproduced using the solid mechanics interface. The masses are defined as circles with rigid material properties, avoiding compression and rotation in order to consider only their inertia. In particular, the value of $m_1$ is set to be equal to 1 g, while $m_2$ depends on the value of the parameter $\alpha$. The circles have depth $d = 0.1$ mm and the radius is set so that $r = \sqrt{\frac{m_1}{\rho_1 \pi d}}$, with $\rho_1 = 1$ kg m$^{-3}$. The lattice constant is $\ell = 20r_1$. The springs are modelled with a spring-damper in matrix form, with translational isotropic coupling $k_n$, depending on the order $n$ of the neighbours. The constant $k_1 = 1$ kN m$^{-1}$ defines the values of $k_2 = \beta k_2$ and $k_3 = \gamma k_1$. The boundary conditions allow to define the wavevector and obtain the desired eigenfrequencies of the systems. The first Brillouin zone is explored along the same path chosen in the previous sections (i.e., $\Gamma - M - K - \Gamma$) in order to verify the closed form solution. The three orders of coupling are studied with values of $\alpha$ equal to 0.5, 1 and 1.5.

Results of the eigenfrequency study are shown in Fig. 8, showing the relation between non-dimensional wavenumber and frequency. All the conditions lead to results that agree perfectly with the analytical predictions of the dispersion diagrams, showing that the closed-form solutions for this system are exact. The analytical prediction in Section 3 highlights an equivalence between in-plane and out-of-plane wave dispersion: this is confirmed numerically as the eigenfrequency study indicates triple degeneracy for each mode. The polarisation of the wave is well defined close to $\Gamma$, allowing discrimination between the degenerate solutions in two in-plane and one out-of-plane eigenmodes. On the contrary, the modes hybridise when moving away from $\Gamma$, leading to coupled in-plane and out-of-plane polarisation. The same curves are obtained for a planar two-dimensional system: removing the out-of-plane degree of freedom simply results in reducing the level of degeneracy to two.

# 6. Edge modes in the presence of beyond-nearest neighbour connections

As mentioned in Section 1, the excellent wave-control properties of MMs make them ideal candidates in the field of topological protection [70–73]. Originally introduced in quantum mechanics [74], topological protection provides new possibilities for wave guiding through symmetry and topology properties of dispersion surfaces [75–77]. Two possibilities exist for achieving topological protection in elasticity. The first consists in using active components, like fluids [78] or rotating gyroscopes [79,80], to break the time-reversal symmetry, while the second employs possible strategies in which the spatial symmetry of the unit cell needs to be broken [81,82]. In the literature, these phenomena have been largely investigated by using local mass-spring systems [83–89] even though the great potential of non-local lattices in this field has been recently demonstrated (cf. Section 1).

It is well known that two-dimensional elastic hexagonal lattices support topologically protected edge modes at the interface between regions of the lattice realised with mirror symmetry. This has been demonstrated in [90] by adding unequal masses at the vertices of an hexagonal lattice and in [44] in the case of a plate-like structure with a hexagonal configuration of mass-spring resonators. A 2D hexagonal lattice model including non-local interactions among masses is investigated in [91] to capture the aforementioned topological effects. Here, similarly to [90], symmetry is broken by setting unequal values of the vertex masses, resulting in the localisation of waves at the corner
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Dispersion curves corresponding to first-, second-, third- and first- plus second- plus third-neighbouring connections for varying values of $\beta$ and $\gamma$. Comparison of the dispersion curves obtained in the case of $m_a = 0.001$ kg, $k_1 = 10^6$ N/m, $\alpha = 0.5$ and: (a) $\gamma = 0.1$, $\beta = 0.1$; (b) $\gamma = 1$, $\beta = 0.1$; (c) $\gamma = 10$, $\beta = 0.1$; (d) $\gamma = 0.1$, $\beta = 1$; (e) $\gamma = 1$, $\beta = 1$; (f) $\gamma = 10$, $\beta = 1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Dispersion curves for the case of first- plus second- plus third-neighbouring connections for varying values of $\alpha$. Comparison of the dispersion curves obtained for $m_a = 0.001$ kg, $k_1 = 10^6$ N/m, $\beta = 10$ and: (a) $\alpha = 0.5$, (b) $\alpha = 1$, (c) $\alpha = 1.5$.</figcaption></figure>

regions. The existence of this anti-symmetric corner state is experi-
mentally confirmed via local pump-probe measurements. Non-trivial
topologically protected bandgaps are finally developed in [92] after
breaking the mirror symmetry of a piezoelectric metamaterial compris-
ing a reconfigurable hexagonal-like honeycomb structure with periodic
bonded piezoelectric beams by means of a negative capacitance circuit.
Specifically, when all piezoelectric circuits are consistent, a symmetric
configuration is obtained, while when the external circuits of two
adjacent piezoelectric beams become inconsistent, the spatial symmetry
is broken. After calculating the Chern number in the system, authors

discover the existence of the interface between different topological
microstructures.

As a final demonstration of how the derived closed-form equations
can be useful to describe wave propagation in more complex systems,
we consider the case of a mass-spring lattice with non-local connections
based on that considered in previous sections, introducing a mirror
symmetry. A first numerical study is related to the band diagram of
supercells of systems with first-, second- and third-neighbour interac-
tions. In particular, all the systems are studied both in standard single
crystal form and in presence of an edge between regions with mirror
<figure><img src="image_10.png" /><figcaption>Fig. 7. Unit cells used to model the system in Comsol Multiphysics FEM simulations. Bloch-Floquet boundary conditions are used to obtain the dispersion diagrams of the system with (a) first-, (b) second- and (c) third-order neighbour interactions. Rigid material conditions are applied to the masses and the connections are modelled as spring-dampers.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Comparison of the dispersion curves obtained with the closed-form analytical solutions (solid line) and the triply degenerate ones obtained with eigenfrequency study with Comsol in FEM simulations (black squares). $\beta = 1$ and $\gamma = 1$ for all the graphs reported in the figure. Plots are displayed for first connections (solid line in blue) and with (a) $\alpha = 0.5$, (b) $\alpha = 1$, (c) $\alpha = 1.5$; second connections (solid line in yellow) and with (d) $\alpha = 0.5$, (e) $\alpha = 1$, (f) $\alpha = 1.5$; third connections (solid line in green) and with (g) $\alpha = 0.5$, (h) $\alpha = 1$, (i) $\alpha = 1.5$.</figcaption></figure>

symmetry. As shown in Fig. 9.a-b-c, all three systems are characterised by the presence of two modes inside the band-gap.

The advantage of the model developed here is that it allows to investigate of the impact of beyond-nearest-neighbour connections on the localisation of edge waves. By introducing second and third levels of interaction within the lattice, we add more degrees of freedom compared to the standard “local” configuration, which could potentially be optimised for specific objectives. Here, the approach is numerical, based on the band analysis presented previously. The FEM model consists in a finite lattice composed of two regions with different orientations (Fig. 10(a)). Numerical band diagrams of a supercell, with $\alpha = 0.5$ in the presence of symmetry breaking are compared to the ones for a single crystal structure in Fig. 9, for networks with first-, second- and third-neighbour connections: the presence of the interface between the two phases leads to the appearance of edge modes inside the band-gap, characterised by an oscillation localised at the edge.

The interface, characterised by symmetry breaking at the unit cell level, creates an interface along which a topologically protected mode can propagate at frequencies that are forbidden for the original structure. We compare results for a lattice with local interactions only and one with first, second, and third-order springs (mixed configuration) in frequency domain simulations. A signal with unit amplitude applied to the left side of the lattice, with frequencies within the band gap for both structures, allows us to analyse the displacement of masses on the opposite end, and to assess localisation properties. Fig. 10(b) shows the obtained waveguiding effect at the interface for a lattice with $\alpha = 0.5$, $\beta = 0.1$, $\gamma = 0.1$: the topologically protected mode propagates with normalised frequency $\omega_0 = 2.21$ (in the band gap for this configuration), from left to right in the finite system. Analysis of the displacement profile of masses on the right end reveals that in the presence of non-local springs, with the chosen parameters, the elastic wave decays more
<figure><img src="image_12.png" /><figcaption>Fig. 9. Numerical band diagrams of supercells with single crystal symmetry (black dots) compared to crystals with an edge between two regions with inverted symmetry (red dots). The band gap opening is obtained with the parameter $α = 0.5$. (a) First-neighbour interactions. (b) Second-neighbour interaction. (c) Third-neighbour interaction.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. (a) Lattice configuration: two realisations of the same lattice with a rotation of 30 deg are joined, and a topologically waveguiding interface is formed. (b) Displacement magnitude distribution obtained with a FEM frequency domain simulation. The normalised frequency $ω_0 = 2.21$ occurs in the band gap for the mixed spring configuration (1st, 2nd and 3rd neighbour connections) with $α = 0.5$, $β = 0.1$, $γ = 0.1$. The load is imposed on the left end of the waveguide. (c) Comparison of the displacement field at the right end of the waveguide for the local and the mixed configurations: the localisation effect is stronger for the considered mixed configuration. (d) Lattice configuration with z-like interface. (e) Displacement field for the topologically protected mode in presence of a z-like interface.</figcaption></figure>

rapidly away from the interface, resulting in improved localisation at the edge (Fig. 10(c)).

Additionally, we examine the topologically protected mode for an angled waveguide at the interface between the two lattice sections. Fig. 10(d) illustrates the test configuration, featuring a z-shaped interface between sublattices. Fig. 10(e) confirms the robustness of the topological waveguiding, since scatter-free wave propagation is observed.

## 7. Conclusions

In conclusion, elastic wave propagation in hexagonal diatomic mass-spring lattices, including non-local interactions up to a generic nth

order has been investigated. Closed-form analytical solutions for the dispersion curves has been obtained as a function of the model parameters, i.e., non-locality order, mass and stiffness characteristics. Results confirm the presence of roton-like behaviour when non-local interactions are present, with the number of local minima reflected in the order of non-locality. Combining first-, second- and third-nearest neighbouring interactions, it is possible to tune the system to obtain the desired complex dispersion characteristics. Numerical FEM simulations are used to validate the analytical solutions and to study the influence of nonlocal interactions on the presence and localisation of edge modes in a mirror symmetric diatomic lattice constituted by two regions with different unit cell orientations. Inclusion of beyond-nearest neighbour
interaction with appropriately chosen spring constants allows to improve the localisation of the edge modes compared to the case where only nearest neighbour interactions are present.

Despite the ideal nature of a mass-spring system, results obtained in our study can be of use to guide the design of more realistic and complex 2D and 3D mechanical structures based on resonating elements, as has been shown in previous studies in the past. In particular, lattices with non-local interactions have been realised experimentally in various systems, for example in an elastic waveguide using simple connecting beams (Meccano elements) for beyond nearest neighbour connections [59], or in an acoustic system using an array of cavities connected beyond nearest neighbour [93]. Our approach provides a simple tool to choose appropriate parameters for the achievement of such experimental systems. Finally, the study of edge mode localisation and topological protection in lattices endowed with nonlocal, roton-like interactions can add a rich new phenomenology to currently extensively studied systems.

### CRediT authorship contribution statement

F. Ongaro: Writing – review & editing, Writing – original draft, Methodology, Formal analysis, Conceptualization. P.H. Beoletto: Writing – review & editing, Software, Validation. F. Bosia: Writing – review & editing, Supervision. M. Miniaci: Writing – review & editing, Conceptualization. N.M. Pugno: Writing – review & editing, Supervision, Conceptualization, Funding acquisition.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

PHB, MM, FB and NMP are supported by the European Union Horizon 2020 FET Open (Boheme) project under grant agreement No. 863179. PHB, FB and NMP acknowledge the financial support of the European Union – Next Generation EU - Piano Nazionale di Ripresa e Resilienza (PNRR) – MISSIONE 4 COMPONENTE 2, INVESTIMENTO N. 1.1, BANDO PRIN 2022 D.D. 104/02-02-2022 (PRIN 2022 2022ATZCJN AMPHYBIA) CUP N.E53D23003040006. MM is supported by the European Union’s Horizon Europe programme in the framework of the ERC StG POSEIDON under Grant Agreement No. 101039576.

## Appendix A

### A.1. Second nearest connections

Consider the hexagonal lattice depicted in Fig. 1 and imagine to connect each mass to its second neighbouring masses by linear springs of constant $k_2$. The resulting system is illustrated in Figs. 1 (b), 1(c), where, for simplicity, the springs are represented by straight yellow lines. Again, connections are shown only for the unit cell (0,0).

As in Section 1, if a time-harmonic vibration is applied to the considered model, writing the equations governing the small-amplitude motion of masses $m_a$ and $m_b$ in the reference cell (0,0) provides

$$ \left\{ \begin{aligned} m_a \ddot{u}_a^{(0,0)} &= k_2 (u_b^{(-1,1)} - u_a^{(0,0)}) + k_2 (u_b^{(-1,-1)} - u_a^{(0,0)}) + k_2 (u_b^{(1,-1)} - u_a^{(0,0)}) \\ m_b \ddot{u}_b^{(0,0)} &= k_2 (u_a^{(1,-1)} - u_b^{(0,0)}) + k_2 (u_a^{(1,1)} - u_b^{(0,0)}) + k_2 (u_a^{(-1,1)} - u_b^{(0,0)}). \end{aligned} \right. \quad (32) $$

Replacing Eqs. (4), (5) into Eq. (32) yields

$$ \left\{ \begin{aligned} -m_a \omega^2 \dot{u}_a &= \hat{u}_a (-3k_2) + \hat{u}_b (e^{i\mathbf{q} \cdot (\mathbf{d}_2 - \mathbf{d}_1)} + e^{-i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (\mathbf{d}_1 - \mathbf{d}_2)}) k_2 \\ -m_b \omega^2 \dot{u}_b &= \hat{u}_a (e^{i\mathbf{q} \cdot (\mathbf{d}_1 - \mathbf{d}_2)} + e^{i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (\mathbf{d}_2 - \mathbf{d}_1)}) k_2 + \hat{u}_b (-3k_2) \end{aligned} \right. \quad (33) $$

that, in matrix form, leads to the eigenproblem stated as

$$ (k_2 \mathbf{K}_2(\mathbf{q}) - \omega^2 \mathbf{M}) \mathbf{U} = \mathbf{0}, \quad (34) $$

with (See Box I) with M and U defined in Eq. (8). Again, in Eq. (33), the factor $e^{-i\omega t}$ is omitted and the more compact notation $\hat{u}_{(-)} := \hat{u}^{(0,0)}(\mathbf{r}_j^{(0,0)})$ is used for simplicity.

The dispersion equation for the case of second-neighbouring connections can be obtained by imposing

$$ \det (k_2 \mathbf{K}_2(\mathbf{q}) - \omega^2 \mathbf{M}) = 0, \quad (36) $$

from which

$$ \omega^4 m_a m_b - 3\omega^2 k_2 (m_a + m_b) + k_2^2 (9 - \Psi_2(\mathbf{q})) = 0, \quad (37) $$

leading to

$$ \omega^2(\mathbf{q}) = k_2 \frac{3(m_a + m_b) \pm \sqrt{9(m_a + m_b)^2 - 4m_a m_b (9 - \Psi_2(\mathbf{q}))}}{2m_a m_b}, \quad (38) $$

with

$$ \Psi_2(\mathbf{q}) := 3 + e^{i\mathbf{q} \cdot 2\mathbf{d}_1} + e^{-i\mathbf{q} \cdot 2\mathbf{d}_1} + e^{i\mathbf{q} \cdot 2\mathbf{d}_2} + e^{-i\mathbf{q} \cdot 2\mathbf{d}_2} + e^{i\mathbf{q} \cdot 2(\mathbf{d}_1 - \mathbf{d}_2)} + e^{-i\mathbf{q} \cdot 2(\mathbf{d}_1 - \mathbf{d}_2)} \quad (39) $$

or, in trigonometric form,

$$ \Psi_2(\mathbf{q}) = 3 + 2\cos(q_x\ell + q_y\ell\sqrt{3}) + 2\cos(2q_x\ell) + 2\cos(q_x\ell - q_y\ell\sqrt{3}). \quad (40) $$

It clearly emerges that the dispersion relations obtained in the case of first- (Eq. (11)) and second- (Eq. (38)) neighbouring interactions have the same form. The only difference is the replacement of the quantities $k_1$ and $\Psi_1(\mathbf{q})$ with $k_2$ and $\Psi_2(\mathbf{q})$. Regarding the latter, note that

$$ \Psi_2(\mathbf{q}) \equiv \Psi_1(2\mathbf{q}), \quad (41) $$

which opened the way to the possibility of establishing a rule of thumb to easily derive the dispersion equation in the general case of nth neighbouring connections. However, before doing so, let us proceed to the case of third neighbouring connections and let us verify if the previous considerations still apply.

### A.2. Third nearest connections

As in section A.1, let us start with the hexagonal lattice in Fig. 1 and let us connect each $m_a$ and $m_b$ mass to its third neighbouring masses by linear springs with stiffness $k_3$. This leads to the mass-spring system represented in Figs. 1(b), 1(c), where the $k_3$ springs are only illustrated for the unit cell (0,0) by blue straight lines.

Applying a time-harmonic vibration to the mass-spring model with third neighbouring interactions and writing the corresponding governing equations of motion for the $m_a$ and $m_b$ masses at the cell (0,0) yields

$$ \left\{ \begin{aligned} m_a \ddot{u}_a^{(0,0)} &= k_3 (u_b^{(1,-3)} - u_a^{(0,0)}) + k_3 (u_b^{(1,1)} - u_a^{(0,0)}) + k_3 (u_b^{(-3,1)} - u_a^{(0,0)}) \\ m_b \ddot{u}_b^{(0,0)} &= k_3 (u_a^{(-1,3)} - u_b^{(0,0)}) + k_3 (u_a^{(-1,-1)} - u_b^{(0,0)}) + k_3 (u_a^{(3,-1)} - u_b^{(0,0)}). \end{aligned} \right. \quad (42) $$
$$ K_2(\mathbf{q}) = \begin{bmatrix} 3 & - (e^{i\mathbf{q} \cdot (\mathbf{d}_2 - \mathbf{d}_1)} + e^{-i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (\mathbf{d}_1 - \mathbf{d}_2)}) \\ -(e^{i\mathbf{q} \cdot (\mathbf{d}_2 - \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (\mathbf{d}_1 - \mathbf{d}_2)}) & 3 \end{bmatrix}, \quad (35) $$

Box I.

$$ K_3(\mathbf{q}) = \begin{bmatrix} 3 & - (e^{i\mathbf{q} \cdot (\mathbf{d}_1 - 3\mathbf{d}_2)} + e^{i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{-i\mathbf{q} \cdot (3\mathbf{d}_1 - \mathbf{d}_2)}) \\ -(e^{-i\mathbf{q} \cdot (\mathbf{d}_1 - 3\mathbf{d}_2)} + e^{-i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (3\mathbf{d}_1 - \mathbf{d}_2)}) & 3 \end{bmatrix}, \quad (45) $$

Box II.

Similarly to the previous sections, introducing Eqs. (4), (5) into Eq. (42) yields the homogeneous system

$$ \left\{ \begin{aligned} -m_a \omega^2 \hat{u}_a &= \hat{u}_a (-3k_3) + \hat{u}_b (e^{i\mathbf{q} \cdot (\mathbf{d}_1 - 3\mathbf{d}_2)} + e^{i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{-i\mathbf{q} \cdot (3\mathbf{d}_1 - \mathbf{d}_2)}) k_3 \\ -m_b \omega^2 \hat{u}_b &= \hat{u}_a (e^{-i\mathbf{q} \cdot (\mathbf{d}_1 - 3\mathbf{d}_2)} + e^{-i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (3\mathbf{d}_1 - \mathbf{d}_2)}) k_3 + \hat{u}_b (-3k_3), \end{aligned} \right. \quad (43) $$

which can be expressed in matrix form by the following eigenproblem (See Box II)

$$ (k_3 K_3(\mathbf{q}) - \omega^2 \mathbf{M}) \mathbf{U} = \mathbf{0}, \quad (44) $$

being **M** and **U** the matrices in Eq. (8) and $\hat{u}_{(\cdot)} := \hat{u}_{(\cdot)}^{(0,0)}(\mathbf{r}_{\mathbf{j}}^{(0,0)})$ for brevity.

Imposing

$$ \det (k_3 K_3(\mathbf{q}) - \omega^2 \mathbf{M}) = 0, \quad (46) $$

gives the dispersion equation

$$ \omega^4 m_a m_b - 3\omega^2 k_3 (m_a + m_b) + k_3^2 (9 - \Psi_3(\mathbf{q})) = 0 \quad (47) $$

resulting in

$$ \omega^2(\mathbf{q}) = k_3 \frac{3(m_a + m_b) \pm \sqrt{9(m_a + m_b)^2 - 4m_a m_b (9 - \Psi_3(\mathbf{q}))}}{2m_a m_b}, \quad (48) $$

with

$$ \begin{aligned} \Psi_3(\mathbf{q}) &:= 3 + e^{i\mathbf{q} \cdot 4a_1} + e^{-i\mathbf{q} \cdot 4a_1} + e^{i\mathbf{q} \cdot 4a_2} + e^{-i\mathbf{q} \cdot 4a_2} + e^{i\mathbf{q} \cdot (4a_1 - a_2)} + e^{-i\mathbf{q} \cdot (4a_1 - a_2)} \\ &= 3 + 2\cos(2q_x \ell + 2q_y \ell \sqrt{3}) + 2\cos(4q_x \ell) + 2\cos(4q_y \ell - 2q_y \ell \sqrt{3}). \end{aligned} \quad (49) $$

Thus, we obtain a dispersion relation having the same form of Eqs. (11) and (38) describing, respectively, the case of first- and second-nearest neighbouring connections. Again, Eq. (48) can be derived by substituting $k_1$ and $\Psi'_1(\mathbf{q})$ with $k_3$ and $\Psi'_3(\mathbf{q})$ in Eq. (11). Also, by comparing Eqs. (39) and (49), we have

$$ \Psi_3(\mathbf{q}) \equiv \Psi_2(2\mathbf{q}). \quad (50) $$

This relation, coupled with the observation $\Psi_2(\mathbf{q}) \equiv \Psi'_1(2\mathbf{q})$ (cf. Section A.1), allowed us to naturally generalise the discussion to the general case of n neighbouring connections. Regarding the latter, it can be said that solving the simplest case of first nearest interactions and deriving, in particular, the corresponding $\Psi'_1(\mathbf{q})$ function, is thus sufficient to quickly obtain the dispersion properties of the lattice with general n-th nearest connections.

Appendix B

Nearest and beyond-the-nearest interactions: simultaneous first-, second- and third-neighbour connections

Let us consider the case of nearest neighbours ($n=1$) together with neighbours with $n=2$ and $n=3$, i.e., first-, second- and third-neighbour interactions. In this configuration, obtained by 'superimposing' the three configurations described in Sections 2 and 3, connections between nearest neighbours are realised by linear springs of stiffness $k_1$ and between non-nearest neighbours by linear springs of stiffness $k_2$ (second neighbours) and $k_3$ (third neighbours).

Applying a time-harmonic vibration to the system and writing the small-amplitude motion equations for the $m_a$ and $m_b$ masses in the cell (0,0) (Fig. 1) leads to

$$ \left\{ \begin{aligned} m_a \ddot{u}_a^{(0,0)} &= k_1(u_b^{(0,0)} - u_a^{(0,0)}) + k_1(u_b^{(-1,0)} - u_a^{(0,0)}) + k_1(u_b^{(0,-1)} - u_a^{(0,0)}) + \\ & \quad + k_2(u_b^{(-1,1)} - u_a^{(0,0)}) + k_2(u_b^{(-1,-1)} - u_a^{(0,0)}) + k_2(u_b^{(1,-1)} - u_a^{(0,0)}) + \\ & \quad + k_3(u_b^{(1,-3)} - u_a^{(0,0)}) + k_3(u_b^{(1,1)} - u_a^{(0,0)}) + k_3(u_b^{(-3,1)} - u_a^{(0,0)}) \\ m_b \ddot{u}_b^{(0,0)} &= k_1(u_a^{(0,0)} - u_b^{(0,0)}) + k_1(u_a^{(1,0)} - u_b^{(0,0)}) + k_1(u_a^{(0,1)} - u_b^{(0,0)}) + \\ & \quad + k_2(u_a^{(1,-1)} - u_b^{(0,0)}) + k_2(u_a^{(1,1)} - u_b^{(0,0)}) + k_2(u_a^{(-1,1)} - u_b^{(0,0)}) + \\ & \quad + k_3(u_a^{(-1,3)} - u_b^{(0,0)}) + k_3(u_a^{(-1,-1)} - u_b^{(0,0)}) + k_3(u_a^{(3,-1)} - u_b^{(0,0)}). \end{aligned} \right. \quad (51) $$

The generalised displacement $\ddot{u}_a^{(\cdot, \cdot)}$ and $\ddot{u}_b^{(\cdot, \cdot)}$ of the two masses can be expressed as a Bloch solution having the form of Eqs. (4) and (5). The latter, substituted into Eq. (51), gives

$$ \left\{ \begin{aligned} -m_a \omega^2 \hat{u}_a &= \hat{u}_a (-3k_1) + \hat{u}_b (1 + e^{-i\mathbf{q} \cdot \mathbf{d}_1} + e^{-i\mathbf{q} \cdot \mathbf{d}_2}) k_1 + \\ & \quad + \hat{u}_a (-3k_2) + \hat{u}_b (e^{i\mathbf{q} \cdot (\mathbf{d}_2 - \mathbf{d}_1)} + e^{-i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (\mathbf{d}_1 - \mathbf{d}_2)}) k_2 + \\ & \quad + \hat{u}_a (-3k_3) + \hat{u}_b (e^{i\mathbf{q} \cdot (\mathbf{d}_1 - 3\mathbf{d}_2)} + e^{i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{-i\mathbf{q} \cdot (3\mathbf{d}_1 - \mathbf{d}_2)}) k_3 \\ -m_b \omega^2 \hat{u}_b &= \hat{u}_a (1 + e^{i\mathbf{q} \cdot \mathbf{d}_1} + e^{i\mathbf{q} \cdot \mathbf{d}_2}) k_1 + \hat{u}_b (-3k_1) + \\ & \quad + \hat{u}_a (e^{i\mathbf{q} \cdot (\mathbf{d}_1 - \mathbf{d}_2)} + e^{i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (\mathbf{d}_2 - \mathbf{d}_1)}) k_2 + \hat{u}_b (-3k_2) + \\ & \quad + \hat{u}_a (e^{-i\mathbf{q} \cdot (\mathbf{d}_1 - 3\mathbf{d}_2)} + e^{-i\mathbf{q} \cdot (\mathbf{d}_2 + \mathbf{d}_1)} + e^{i\mathbf{q} \cdot (3\mathbf{d}_1 - \mathbf{d}_2)}) k_3 + \hat{u}_b (-3k_3), \end{aligned} \right. \quad (52) $$

with $\hat{u}_{(\cdot, \cdot)}$, $\mathbf{q}$ and $\mathbf{d}_j$, $j=1,2$, the quantities previously introduced (cf. Section 1).
<figure><img src="image_14.png" /><figcaption>Fig. 11. Local plus non-local interactions: dispersion curves for first- plus second- plus third-neighbouring connections with $m_a = 0.001$ kg, $k_1 = 10^6$ N/m, $\beta = 0.1$ and varying values of $\alpha$: (a) $\alpha = 0.5$, (b) $\alpha = 1$, (c) $\alpha = 1.5$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Comparison of the dispersion curves obtained for first- plus second- plus third-neighbouring connections with $m_a = 0.001$ kg, $k_1 = 10^6$ N/m, $\beta = 1$ and varying values of $\alpha$ and $\gamma$: (a) $\alpha = 0.5$, (b) $\alpha = 1$, (c) $\alpha = 1.5$.</figcaption></figure>

Eq. (52) can be written in matrix form, thus leading to the eigenproblem

$$ (k_1 \mathbf{K}_1(\mathbf{q}) + k_2 \mathbf{K}_2(\mathbf{q}) + k_3 \mathbf{K}_3(\mathbf{q}) - \omega^2 \mathbf{M}) \mathbf{U} = \mathbf{0}, \quad (53) $$

being the matrices $\mathbf{M}$, $\mathbf{U}$, $\mathbf{K}_i(\mathbf{q})$, $i = 1, 2, 3$, defined in Eqs. (8), (35) and (45). Imposing

$$ \det (k_1 \mathbf{K}_1(\mathbf{q}) + k_2 \mathbf{K}_2(\mathbf{q}) + k_3 \mathbf{K}_3(\mathbf{q}) - \omega^2 \mathbf{M}) = 0, \quad (54) $$

yields the dispersion relation of the investigated system

$$ \omega^4 m_a m_b - 3\omega^2 (k_1 + k_2 + k_3) (m_a + m_b) + k_1^2 (9 - \Psi_1(\mathbf{q})) + k_2^2 (9 - \Psi_2(\mathbf{q})) \\ + k_3^2 (9 - \Psi_3(\mathbf{q})) + 2k_1 k_2 (9 - \Psi_{12}(\mathbf{q})) + 2k_1 k_3 (9 - \Psi_{13}(\mathbf{q})) + 2k_2 k_3 (9 - \Psi_{23}(\mathbf{q})) = 0, \quad (55) $$

with $\Psi_1, \Psi_2, \Psi_3$ the functions in Eqs. (12), (39), (49),

$$ \Psi_{12}(\mathbf{q}) := \frac{eb + af}{2}, \quad \Psi_{13}(\mathbf{q}) := \frac{ed + cf}{2}, \quad \Psi_{23}(\mathbf{q}) := \frac{ad + bc}{2}, \quad (56) $$

and

$$ a := a(\mathbf{q}) = - (e^{i\mathbf{q} \cdot (\mathbf{a}_1 - \mathbf{a}_2)} + e^{i\mathbf{q} \cdot (\mathbf{a}_1 + \mathbf{a}_2)} + e^{i\mathbf{q} \cdot (\mathbf{a}_2 - \mathbf{a}_1)}) = \mathbf{K}_2(2, 1), \quad (57) $$

$$ b := b(\mathbf{q}) = - (e^{i\mathbf{q} \cdot (\mathbf{a}_2 - \mathbf{a}_1)} + e^{-i\mathbf{q} \cdot (\mathbf{a}_1 + \mathbf{a}_2)} + e^{i\mathbf{q} \cdot (\mathbf{a}_1 - \mathbf{a}_2)}) = \mathbf{K}_2(1, 2), \quad (58) $$

$$ c := c(\mathbf{q}) = - (e^{-i\mathbf{q} \cdot (\mathbf{a}_1 - 3\mathbf{a}_2)} + e^{-i\mathbf{q} \cdot (\mathbf{a}_1 + \mathbf{a}_2)} + e^{i\mathbf{q} \cdot (3\mathbf{a}_1 - \mathbf{a}_2)}) = \mathbf{K}_3(2, 1), \quad (59) $$

$$ d := d(\mathbf{q}) = - (e^{i\mathbf{q} \cdot (\mathbf{a}_1 - 3\mathbf{a}_2)} + e^{i\mathbf{q} \cdot (\mathbf{a}_1 + \mathbf{a}_2)} + e^{-i\mathbf{q} \cdot (3\mathbf{a}_1 - \mathbf{a}_2)}) = \mathbf{K}_3(1, 2), \quad (60) $$

$$ e := e(\mathbf{q}) = - (1 + e^{i\mathbf{q} \cdot \mathbf{a}_1} + e^{i\mathbf{q} \cdot \mathbf{a}_2}) = \mathbf{K}_1(2, 1), \quad (61) $$

$$ f := f(\mathbf{q}) = - (1 + e^{-i\mathbf{q} \cdot \mathbf{a}_1} + e^{-i\mathbf{q} \cdot \mathbf{a}_2}) = \mathbf{K}_1(1, 2). \quad (62) $$

By adopting a more compact notation, Eq. (55) can be expressed by

$$ \omega^4 m_a m_b - 3\omega^2 \sum_{i=1}^{3} k_i (m_a + m_b) + \sum_{i=1}^{3} k_i^2 (9 - \Psi_i(\mathbf{q})) $$
<figure><img src="image_16.png" /><figcaption>Fig. 13. Study of the effect of β on crystals with simultaneous first- and second-neighbour interaction. (a-b-c) Analytical dispersion curves of the system with α = 1 and values of β = 0.1, 1, 10 respectively. The level ω<sub>max</sub> corresponds to the maximum frequency of the lowest eigenmodes in the Γ–M path while ω<sub>M</sub> is the frequency of this eigenmode in the Γ point. (d) Trend of the ratio ω<sub>M</sub>/ω<sub>max</sub> for increasing values of β. (e) Lowest eigenmode of the system for values of β ranging from 10<sup>-2</sup> to 10<sup>5</sup>.</figcaption></figure>

$$+ 2 \sum_{i=1}^{2} \left( \sum_{j=i+1}^{3} k_i k_j (9 - \Psi_{ij}(\mathbf{q})) \right) = 0, \qquad (63)$$

from which

$$\omega^2(\mathbf{q}) = \frac{3 \sum_{i=1}^{3} k_i (m_a + m_b) \pm \sqrt{9 \left(\sum_{i=1}^{3} k_i (m_a + m_b)\right)^2 - 4m_a m_b (\Lambda_i + 2\Lambda_{ij})}}{2m_a m_b}, \qquad (64)$$

where

$$\Lambda_i := \sum_{i=1}^{3} k_i^2 (9 - \Psi_i(\mathbf{q})), \quad \Lambda_{ij} := \sum_{i=1}^{2} \left( \sum_{j=i+1}^{3} k_i k_j (9 - \Psi_{ij}(\mathbf{q})) \right). \quad (65)$$

## Appendix C

See Figs. 11 and 12.

## Appendix D

Following the analytical approach presented in Appendix B and limiting the study to the simultaneous presence of first and second neighbour connections, it is possible to evaluate the effect of the parameter β = k₂/k₁ on the dispersion relation. In Fig. 13.a-b-c the dispersion curves for values of β = 0.1, 1, 10 are presented, for a crystal with α = 1: it is clear that for increasing value of β the non-local effects have a more relevant impact as the springs related to the second-neighbour connections prevail. While for β = 0.1 the curves resemble the ones of the crystal with local connections, when the ratio k₁/k₂ increases it is possible to observe the roton-like slope inversion at the M point together with the appearance of a second Dirac cone.

A quantitative analysis of the roton behaviour is proposed by examining the ratio ω<sub>M</sub>/ω<sub>max</sub>, which represents the frequency of the lowest eigenmode at the M point relative to the maximum frequency of the

same eigenmode along the Γ–M path. In a purely local system, this ratio equals 1, as the frequency trend in this region of the irreducible Brillouin zone (IBZ) is monotonic. However, the onset of slope inversion reduces the value of ω<sub>M</sub>, thereby lowering the ratio. As illustrated in Fig. 13.d, the roton emerges for values of β approaching 1, with the ratio decreasing further and approaching 0 as β tends towards infinity.

Fig. 13.e shows the lowest eigenmode for the whole set of values of β that were considered for this study: it is possible to observe that, while the system it is very similar to the local crystal for β ≪ 1, if β tends to infinity the dispersion curves reproduce the ones of the structure with only second-neighbour connections.

## Data availability

No data was used for the research described in the article.

## References

[1] Liu Z, Zhang X, Mao Y, Zhu YY, Yang Z, Chan CT, et al. Locally resonant sonic materials. Science 2000;289:1734–6.

[2] Lu MH, Feng L, Chen YF. Phononic crystals and acoustic metamaterials. Mater Today 2009;12:34–42.

[3] Maldovan M. Sound and heat revolutions in phonics. Nature 2013;503:209–17.

[4] Cummer SA, Christensen J, Alu A. Controlling sound with acoustic metamaterials. Nat Rev Mater 2016;1.

[5] Ma G, Sheng P. Acoustic metamaterials: from local resonances to broad horizons. Sci Adv 2016;1.

[6] Jiao P, Mueller J, Raney JR, Zheng X, Alavi AH. Mechanical metamaterials and beyond. Nat Commun 2023;14.

[7] Wu L, Wang Y, Chuang K, Wu F, Wang Q, Lin W, et al. A brief review of dynamic mechanical metamaterials for mechanical energy manipulation. Mater Today 2021;44:168–93.

[8] Craster RV, Guenneau SR, Muamer K, Wegener M. Mechanical metamaterials. Rep Progr Phys 2023.
[9] Lee JH, Singer JP, Thomas EL. Micro-nanostructured mechanical metamaterials. Adv Mater 2012;24(36):4782-810.

[10] Bertoldi K, Vitelli V, Christensen J, van Hecke M. Flexible mechanical metamaterials. Nat Rev Mater 2017;2.

[11] Kadic M, Milton GW, van Hecke M, Wegener M. 3D metamaterials. Nat Rev Phys 2019;1:198-210.

[12] Kahler H, Platz D, Schmid S. Surface acoustic wave coupling between micromechanical resonators. Commun Phys 2022;5.

[13] Zhu W, Liu A, Bourouina T, et al. Microelectromechanical Maltese-cross metamaterial with tunable terahertz anisotropy. Nat Commun 2012;3.

[14] Gao N, Zhang Z, Deng J, Guo X, Cheng B, Hou H. Acoustic metamaterials for noise reduction: a review. Adv Mater Technol 2022;7.

[15] Rui S, Zhang W, Yu R, Wang X, Ma F. A multi-band elastic metamaterial for low-frequency multi-polarization vibration absorption. Mech Syst Signal Process 2024;216.

[16] Miniaci M, Krushynska A, Bosia F, Pugno NM. Large scale mechanical metamaterials as seismic shields. New J Phys 2016;18.

[17] Colombi A, Colquitt D, Roux P, Guenneau S, Craster RV. A seismic metamaterial: the resonant metawedge. Sci Rep 2016;6:1-6.

[18] Palermo A, Kradel K, Marzani A, Daraio C. Engineered metabarrier as shield from seismic surface waves. Sci Rep 2016;6:1-10.

[19] Miniaci M, Kherraz N, Croenne C, Mazzotti M, Morvarid M, Gliozzi AS, et al. Hierarchical large-scale elastic metamaterials for passive seismic wave mitigation. EPJ Appl Metamaterials 2021;8:14.

[20] Nistri F, Bosia F, Gliozzi AS, D'Alessandro A, Caverni S, Charkaluk P, et al. Design and in field validation of a modular metamaterial for mitigation of railway induced vibrations. Soil Dyn Earthq Eng 2024;180.

[21] Huber SD. Topological mechanics. Nat Phys 2016;12:621-3.

[22] Huang H, Chen J, Huo S. Recent advances in topological elastic metamaterials. J Phys Condens Matter 2021;33:503002.

[23] Fleury R, Khanikaev AB, Alù A. Floquet topological insulators for sound. Nat Commun 2016;7:11744.

[24] He C, Ni X, Ge H, Sun XC, Bin Chen Y, Lu MH, et al. Acoustic topological insulator and robust one-way sound transport. Nat Phys 2016;12:1124-9.

[25] Lu J, Qiu C, Ye L, Fan X, Ke M, Zhang F, et al. Observation of topological valley transport of sound in sonic crystals. Nat Phys 2017;13:369-74.

[26] Mousavi SH, Khanikaev AB, Wang Z. Topologically protected elastic waves in phononic metamaterials. Nat Commun 2015;6:8682.

[27] Fang H, Xie G, Huang H, et al. Coupled topological rainbow trapping of elastic waves in two-dimensional phononic crystals. Sci Rep 2024;14:17011.

[28] Sui F, Chen J, Huang H. Tunable topological edge states and rainbow trapping in two dimensional magnetoelastic phononic crystal plates based on an external magnetostatic field. Int J Mech Sci 2022;225:107360.

[29] Hussein MI, Leamy MJ, Ruzzene M. Dynamics of phononic materials and structures: Historical origins. Recent Prog Futur Outlook Appl Mech Rev 2014;66:040802.

[30] Jensen JS. Phononic band gaps and vibrations in one- and two-dimensional mass-spring structures. J Sound Vib 2003;266:1053-78.

[31] Hussein MI. Theory of damped Bloch waves in elastic media. Phys Rev B 2009;80.

[32] Huang H, Huo S, Chen J. Subwavelength elastic topological negative refraction in ternary locally resonant phononic crystals. Int J Mech Sci 2021;198:106391.

[33] Laude V. Principles and properties of phononic crystal waveguides. APL Mater 2021;9.

[34] Farhat M, Guenneau S, Enoch S. Ultrabroadband elastic cloaking in thin plates. Phys Rev Lett 2009;103.

[35] Deng B, Wang P, He Q, Tournat V, Bertoldi K. Metamaterials with amplitude gaps for elastic solitons. Nat Commun 2018;9.

[36] Onorato MN, Lvov YV, Dematteis G, Chibbaro S. Wave turbulence and thermalization in one-dimensional chains. Phys Rep 2023;1040:1-36.

[37] Chen H, Nassar H, Huang GL. A study of topological effects in 1D and 2D mechanical lattices. J Mech Phys Solids 2018;117:22-36.

[38] Chaplain GJ, De Ponti JM, Aguzzi G, Colombi A, Craster RV. Topological rainbow trapping for elastic energy harvesting in graded su-schrieffer-heeger systems. Phys Rev Appl 2020;14.

[39] Lu L, Joannopoulos JD, Solja M. Topological photonics. Nat Photonics 2014;8:821-9.

[40] Xue H, Yang Y, Zhang B. Topological acoustics. Nat Rev Mater 2022;7.

[41] Sui F, Chen J, Huang H. Tunable topological edge states and rainbow trapping in two dimensional magnetoelastic phononic crystal plates based on an external magnetostatic field. Int J Mech Sci 2022;225:107360.

[42] Yu G, Wang W, Lai L, Peng T, Jiang C, Li Y. Thermally reconfigurable topological protection routes for the 2D phononic system. Int J Mech Sci 2023;241:107989.

[43] Veenstra J, Gamayun O, Guo X, et al. Non-reciprocal topological solitons in active metamaterials. Nature 2024;627:528-33.

[44] Pal RK, Ruzzene M. Edge waves in plates with resonators: an elastic analogue of the quantum valley hall effect. New J Phys 2017;19.

[45] Zheng LY, Theocharis G, Tournat V, Gusev V. Quasitopological rotational waves in mechanical granular graphene. Phys Rev 2018;97:060101.

[46] Zhou Y, Bandaru PR, Sievenpiper DF. Quantum-spin hall topological insulator in a spring-mass system. New J Phys 2018;20:123011.

[47] Grundmann M. Topological states due to third-neighbour coupling in diatomic linear elastic chains. Phys Status Solidi (B) 2020;257.

[48] Han L, Xuequin H, Mou Y, Jiuyang L, Weiyin D, Zhengyou L. Acoustic topological metamaterials of large winding number. Phys Rev Appl 2023;19.

[49] Betancur-Ocampo Y, Manjarrez-Montanez B, Martínez-Arguello A, Mendez-Sanchez RA. Twofold topological phase transitions induced by third-nearest-neighbour hopping in one-dimensional chains. Phys Rev B 2024;109.

[50] Wang S, Wang YZ. Nonlocal topological states in elastic wave metamaterials with active feedback control. Int J Solids Struct 2023;281.

[51] Wang S, Wang YZ. Active control on topological interface states of elastic wave metamaterials with double coupled chains. J Acoust Soc Am 2023;154:2440-52.

[52] Iglesias Martinez JA, Gross MF, Chen Y, Frenzel T, Laude V, Muamer K, et al. Experimental observation of roton-like dispersion relations in metamaterials. Sci Adv 2021;7.

[53] Wang K, Chen Y, Muamer K, Wang C, Wegener M. Nonlocal interaction engineering of 2D roton-like dispersion relations in acoustic and mechanical metamaterials. Commun Mater 2022;3:35.

[54] Zhu Z, Zhen G, Gui-Geng L, Yong G, Yin W, Xiang X, et al. Observation of multiple rotons and multidirectional roton-like dispersion relations in acoustic metamaterials. New J Phys 2022;24.

[55] Landau L. Theory of superfluidity of helium II. Phys Rev 1941;60:356-8.

[56] Fleury R. Non-local oddities. Nat Phys 2021;17:760-72.

[57] Zhu Z, Gao Z, Liu GG, Ge Y, Wang Y, Xi X, et al. Observation of multiple rotons and multidirectional rotor-like dispersion relations in acoustic metamaterials. New J Phys 2022;24.

[58] Kazemi A, Deshmukh KJ, Chen F, Liu Y, Deng B, Fu HC, et al. Non-local phononic crystals for dispersion customisation and undulation-point dynamics. Phys Rev Lett 2023;131.

[59] Chaplain GJ, Hooper IR, Hibbins AP, Starkey TA. Reconfigurable elastic metamaterials: Engineering dispersion with beyond nearest neighbours. Phys Rev Appl 2023;19.

[60] Chen Y, Kadic M, Wegener M. Roton-like acoustical dispersion relations in 3D metamaterials. Nat Commun 2021;12.

[61] Chen Y, Wang K, Kadic M, Guenneau S, Wang C, Wegener M. Phonon transmission through a nonlocal meta material slab. Commun Phys 2023;6.

[62] Peri V. Experimental characterisation of fragile topology in an acoustic meta material. Science 2020;367.

[63] Sepehri S, Mashhadi MM, Fakhrabadi MMS. Nonlinear nonlocal phononic crystals with rotten-like behavior. Nonlinear Dynam 2023;111:8591-610.

[64] Farzbod F, Scott-Emuakpor OE. Interactions beyond the nearest neighbours in a periodic structure: force analysis. Int J Solids Struct 2020;199:203-11.

[65] Wang J, Huang Y, Chen W, Zhu W. Abnormal wave propagation behaviours in two-dimensional mass-spring structures with nonlocal effects. Math Mech Solids 2019;24:3632-43.

[66] Martinez JAI, Grob MF, Chen Y, Frenzel T, Laude V, Kadic M, et al. Experimental observation of roton-like dispersion relations in metamaterials. Sci Adv 2021;7.

[67] Dal Poggetto VF, Pal RK, Pugno NM, Miniaci M. Topological bound modes in phononic lattices with nonlocal interactions. Int J Mech Sci 2024;281:109503.

[68] Bloch F. Über die Quantenmechanik der Elektronen in Kristallgittern. Z Phys 1929;52:555-600.

[69] Maurin F, Claeys C, Deckers E, Desmet W. Probability that a band-gap extremum is located on the irreducible Brillouin-zone contour for the 17 different plane crystallographic lattices. Int J Solids Struct 2018;35:26-36.

[70] Zangeneh-Nejad F, Alu A, Fleury R. Topological wave insulators: a review. Phys 2020;21:467-99.

[71] Khanikaev AB, Alu A. Topological photonics: robustness and beyond. Nat Commun 2024;151.

[72] Zarate Y, Babaee S, Kang SH, Neshev DN, Shadrivov IV, Bertoldi K, et al. Elastic metamaterials for tuning circular polarization of electromagnetic waves. Sci Rep 2016;6.

[73] Ma G, Fu C, Wang G, Del Hougne P, Christensen J, Lai Y, et al. Polarization bandgaps and fluid-like elasticity in fully solid elastic metamaterials. Nat Commun 2016;7.

[74] Hasan MZ, Kane CL. Colloquium: topological insulators. Rev Modern Phys 2010;82.

[75] Xia BZ, Liu TT, Huang GL, Dai HQ, Jiao JR, Zang XG, et al. Topological phononic insulator with robust pseudospin-dependent transport. Phys Rev B 2017;96.

[76] Jin Y, Torrent D, Djafari-Rouhani B. Robustness of conventional and topologically protected edge states in phononic crystal plates. Phys Rev B 2018;98.
[77] Wang G, Chen Z, Shi Z, Lim C. Three-dimensional acoustic metamaterials with topological states of different orders and multidirectional waveguiding. Phys Rev B 2024;109.

[78] Souslov A, Van Zuiden BC, Bartolo D, Vitelli V. Topological sound in active-liquid metamaterials. Nat Phys 2017;13:1091–4.

[79] Nash LM, Kleckner D, Read A, Vitelli V, Turner AM, Irvine WT. Topological mechanics of gyroscopic metamaterials. Proc Natl Acad Sci 2015;112.

[80] Wang P, Lu L, Bertoldi K. Topological phononic crystals with one-way elastic edge waves. Phys Rev Lett 2015;115.

[81] Ma G, Xiao M, Chan CT. Topological phases in acoustic and mechanical systems. Nat Rev Phys 2019;1:281–94.

[82] Miniaci M, Pal R. Design of topological elastic waveguides. J Appl Phys 2021;130.

[83] Kane CL, Lubensky TC. Topological boundary modes in isostatic lattices. Nat Phys 2014;10:39–45.

[84] Rocklin DZ, Chen BG, Falk M, Vitelli V, Lubensky T. Mechanical Weyl modes in topological Maxwell lattices. Phys Rev Lett 2013;116.

[85] Pal RK, Schaeffer M, Ruzzene M. Helical edge states and topological phase transitions in phononic systems using bi-layered lattices. J Appl Phys 2016;119.

[86] Lu J, Qiu C, Ye L, Fan X, Ke M, Zhang F, et al. Observation of topological valley transport of sound in sonic crystals. Nat Phys 2017;13:369–74.

[87] Du Z, Chen H, Huang G. Optimal quantum valley hall insulators by rationally engineering Berry curvature and band structure. J Mech Phys Solids 2020;135.

[88] Al Ba ba H, Yu K, Wang Q. Elastically-supported lattices for tunable mechanical topological insulators. Extrem Mech Lett 2020;38.

[89] Zhou W, Su Y, Chen W, Lim C, et al. Voltage-controlled quantum valley Hall effect in dielectric membrane-type acoustic metamaterials. Int J Mech Sci 2020;172.

[90] Vila J, Pal RK, Ruzzene M. Observation of topological valley modes in an elastic hexagonal lattice. Phys Rev B 2017;96.

[91] An S, Liu T, Fan H, Gao H, Gu Z, Liang S, et al. Second-order elastic topological insulator with valley-selective corner states. Int J Mech Sci 2022;224.

[92] Chen Y, Liu D, Wu Y, Yu P, Liu Y. Valley hall elastic topological insulator with large chern numbers. Int J Mech Sci 2023;239.

[93] Moore DB, Sambles JR, Hibbins AP, Starkey TA, Chaplain GJ. Acoustic surface modes on metasurfaces with embedded next-nearest-neighbor coupling. Phys Rev B 2023;107:144110.