# Tailoring dispersion and evanescent modes in multimodal nonlocal lattices using positive-only interactions

Lucas Rouhi, Christophe Droz

► To cite this version:

Lucas Rouhi, Christophe Droz. Tailoring dispersion and evanescent modes in multimodal nonlocal lattices using positive-only interactions. 2025. hal-05330430

HAL Id: hal-05330430

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

Preprint submitted on 24 Oct 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
# Tailoring dispersion and evanescent modes in multimodal nonlocal lattices using positive-only interactions

Lucas Rouhi¹ Christophe Droz¹

¹Univ. Gustave Eiffel, Inria, COSYS-SII, I4S, Rennes, France

lucas.rouhi@inria.fr

## Highlights

* Customization of band diagrams through nonlocal interactions is key to advanced phononic design.

* Versatility of the method enables selective customization of Rotons, ZGV and GVD indicators.

* The method applies to a broad range of lattices with nonlocal spring and beam interactions.

* We achieved passive and stable mechanical systems through real, positive-only nonlocal stiffnesses.

* We control multimodal evanescent decay inside bandgaps by tuning the complex wavenumbers.

## Abstract

Metamaterials derive their unconventional properties from engineered microstructures, with periodic lattices providing a versatile framework for modeling wave propagation. Dispersion relations, obtained from Bloch-Floquet theory, govern how waves propagate, attenuate, or localize within such systems. Extending interactions beyond nearest neighbors, through nonlocality, substantially enriches the design space of band diagrams, enabling phenomena such as negative or zero group velocities, roton-like extrema, and band-gap localization. However, existing approaches to dispersion tailoring often rely on analytical formulations or Fourier-based identifications, which become impractical for complex coupling mechanisms and offer limited control over physical constraints such as stiffness positivity.

This work introduces a general interpolation-based framework for customizing dispersion relations in uniform nonlocal lattices. Rather than reconstructing full dispersion curves, the method enforces prescribed frequency-wavenumber points as interpolation constraints, enabling localized and tunable control of wave behavior. The formulation is applied to both spring- and beam-interaction lattices, and demonstrated on an Euler-Bernoulli beam model with adjustable nonlocal couplings. Through systematic parameter tuning, the framework enables the creation of rotons, the adjustment of group-velocity dispersion, and the design of evanescent waves with controlled exponential decay within band gaps, all while ensuring real, positive-only stiffness parameters and passive mechanical behavior.

Altogether, this parametric interpolation strategy provides a physically consistent and computationally efficient route for engineering advanced phononic functionalities in periodic nonlocal systems.

**Keywords:** Nonlocal lattices, Dispersion engineering, Wave propagation control, Roton modes, Group velocity dispersion, Positive-only nonlocal stiffnesses, Evanescent decay tuning

## 1 Introduction

Metamaterials are artificially engineered composites whose microstructural design grants access to physical properties absent in conventional materials. Among them, phononic metamaterials stand out for their ability to manipulate mechanical waves through carefully designed periodic architectures. Such structural periodicity enables tailored control of wave propagation, with key applications in vibration damping, noise mitigation, and acoustic cloaking. As a result, metamaterials play an essential role across diverse industries, ranging from civil engineering and automotive design to aerospace and acoustic insulation systems.
A common mathematical abstraction for modeling such materials is the *periodic lattice*, obtained by infinitely repeating an elementary unit cell. These unit cells may range from simple monoatomic arrangements to more elaborate multi-mass structures such as diatomic lattices. Exploiting periodicity, the dynamics of the infinite lattice reduce to supercells governed by setsof ordinary differential equations [1, 2], providing an analytically tractable yet powerful framework for investigating wave phenomena.

At the heart of this analysis lies the Bloch-Floquet theory, which expresses solutions as harmonic waves mod-
ulated by spatially periodic functions. This formalism naturally yields dispersion relations connecting angular
frequency and wavenumber. Dispersion curves, in turn, are indispensable tools for identifying propagating,
evanescent, and attenuating modes, thereby offering deep insight into transmission and localization mecha-
nisms.

While conventional models typically assume nearest-neighbor couplings, nonlocal lattice models have attracted increasing attention since Brillouin's pioneering work in 1946 [1], whose archetypal example, a monoatomic chain with long-range spring interactions, is illustrated in Figure 1. His original formulation of nonlocal lattices has been recently revived and extended in various modern contexts, highlighting their potential for advanced control of wave propagation. By extending interactions beyond nearest neighbors, nonlocality introduces additional tunable parameters and thereby greatly enriches the design space of dispersion relations. This expanded framework allows access to remarkable wave phenomena, including negative or hypersonic group velocities [3, 4], roton-like dispersion profiles [5–8], dirac cones [9], frozen evanescent modes [10, 11], metadamping [12] and the emergence of additional transient modes [1, 13]. Such versatility has established nonlocality as a promising strategy for tailoring wave behavior in periodic media.

<figure><img src="image_2.png" /><figcaption>Figure 1: Monoatomic chain with non-local interactions up to order $P = 2$. Local interactions are shown in black, while higher-order couplings are highlighted in blue.</figcaption></figure>

The design of nonlocal lattices can be understood as a two-step process: selecting a geometry of nonlocality to be periodically repeated, and specifying the physical nature of the interactions. Several nonlocal lattice geometries have been explored in the literature. The most common configuration connects each node to its *P*-th neighbors [13–17], while other works have considered diatomic extensions introducing two distinct species of nodes [13, 18], or nonlocal architectures incorporating local resonators [19]. Different physical realizations of the coupling mechanisms have also been investigated. Most studies employ spring-based interactions [13, 15–17], but alternative implementations have been proposed, including networks of coaxial cables [20], reconfigurable assemblies using model-construction kits such as Meccano™ [21], electromagnetic forces [4], and Euler-Bernoulli beam interactions [14].

The central challenge lies in precisely shaping dispersion relations. Brillouin [1] already recognized that the dispersion relation of a generalized nonlocal spring lattice could be expressed as a Fourier series, with spring stiffnesses corresponding, up to scaling, to Fourier coefficients. In principle, this identification allows any target dispersion curve to be approximated by assigning stiffness values according to Fourier decomposition. However, not all dispersion profiles are physically realizable. Kazemi et al. [13] clarified the feasibility conditions of this approach, showing that while certain constraints exist, such as linear behavior at low wavenumbers and vanishing group velocity at high wavenumbers, they remain relatively mild.

Subsequent extensions included Kazemi's generalization of Brillouin's method to diatomic lattices with spring-based nonlocality, and Paul's formulation of dual nonlocal interactions [22]. More recently, Edge et al. [17] proposed an initial approach for parameterizing evanescent modes in nonlocal spring-interaction models. Yet, these strategies remain difficult to apply when the coupling mechanism departs from simple springs. They rely on closed-form analytical expressions of the dispersion relation, often unavailable in practice, and on linear mappings between stiffnesses and Fourier coefficients that break down for more complex interaction laws, leading instead to nonlinear identification problems.

In addition, nothing in the general identification process guarantees that the resulting stiffness coefficients remain positive, which becomes problematic when aiming to design passive and mechanically stable structures. Negative stiffness values would correspond to active components that inject energy into the system, leading
to loss of passivity and potential instability. Edge et al. [17] proposed sufficient criteria on admissible target dispersion curves to ensure positive stiffnesses, but these conditions appear overly restrictive and limit the generality of the approach. This limitation stems from the very objective of reconstructing an entire target dispersion curve, which tends to overconstrain the design problem. In many practical contexts, however, the interest lies not in reproducing a full curve but in tailoring specific points or regions that govern distinctive wave phenomena, such as rotons [5–8], frozen evanescent modes [10, 11], or Dirac cones [9], or in locally shaping the dispersion to control group-velocity dispersion or band-gap properties [17].

Motivated by these limitations, the present work introduces a complementary paradigm for dispersion tailoring. Rather than reconstructing an entire target curve, the proposed approach focuses on localized customization: specific frequency-wavenumber points are prescribed as interpolation constraints that the dispersion relation must satisfy. This formulation enables direct and flexible control over selected regions of the spectrum without relying on a global curve-fitting procedure. The general framework is first established and then applied to the Euler-Bernoulli beam lattice model introduced by Edge et al. [14], which is here reformulated to incorporate tunable nonlocal couplings and multiple degrees of freedom per node. The resulting multimodal dynamics allow simultaneous control of several interacting wave branches within a unified setting. Several representative applications are presented to illustrate the method’s capabilities, including the creation of rotons, the tuning of group-velocity dispersion, and the controlled design of evanescent modes within band gaps. In addition, the interpolation scheme naturally accommodates design constraints, such as enforcing positive stiffness parameters, to ensure passive and mechanically stable interactions at all stages of the tuning process. Altogether, the framework provides a unified and physically consistent foundation for engineering wave phenomena in nonlocal lattices: it applies seamlessly to both spring and beam systems, guarantees stability through real positive stiffness coefficients, and enables fine control of propagation and attenuation mechanisms via complex dispersion analysis.

## 2 Tuning of dispersion relations in uniform nonlocal lattices

The ability to tailor the dispersion characteristics of nonlocal lattices systems comes from a simple yet fundamental observation: each nonlocal interaction, regardless of its underlying physical nature, is governed by one or more tunable parameters. Introducing additional nonlocal interactions effectively increases the number of degrees of freedom in the system, which can in turn be exploited to shape the dispersion curves in a controlled manner.

The approach developed here follows a two-step rationale. First, a formal description of the class of *uniform nonlocal lattices* is introduced, providing a precise mathematical framework in which the dispersion relation can be expressed and analyzed. Second, a systematic method is proposed for selecting the interaction parameters so as to enforce prescribed points on the dispersion curves, thus achieving the desired wave-propagation properties.

### 2.1 Uniform nonlocal models

The lattices considered in this work are infinite, one-dimensional structures: their nodes are aligned along a single axis, and waves propagate exclusively in this direction. To introduce the concept of nonlocality, let the nodes be equally spaced at a distance $l$. A reference node is denoted by $x_0$, and the next node along the positive direction is $x_1$, located at a distance $l$ from $x_0$. This indexing extends such that $x_n$ denotes the node located at $nl$ from $x_0$ in the positive direction, and $x_{-n}$ in the opposite direction.

A nonlocal interaction of order $p$ couples two nodes separated by a distance $pl$, that is, between $x_n$ and $x_{n+p}$. A lattice is said to be nonlocal of order $P$ if interactions extend to at most $P$ neighbors. In the present study, the focus is on *uniform nonlocal interactions*, meaning that if an order-$p$ interaction exists between a given pair of nodes, then every pair of nodes separated by $pl$ is coupled by the same type of interaction, characterized by identical parameters.

Once the interaction model between nodes is fixed, a uniform nonlocal lattice of order $P$ is entirely characterized by its maximum order $P$ and the set of parameters $\beta = (\beta_p)$ configuring each interaction of order $p$. These parameters are generally real and positive, but may acquire a imaginary part in the presence of damping.

Each node possesses $d$ mechanical degrees of freedom, and the displacement of node $x_n$ from equilibrium is denoted by $\mathbf{u}_n \in \mathbb{R}^d$. For a uniform nonlocal lattice of order $P$, the equations of motion for each $\mathbf{u}_n$ take the
generic form

$$ \mathbf{M} \partial_t^2 \mathbf{u}_n + \mathbf{K} [\mathbf{u}_{n-P} \dots \mathbf{u}_n \dots \mathbf{u}_{n+P}]^T = 0 \quad (2.1) $$

where **M** and **K** are matrices of size $d \times d$ and $d \times (2P + 1)$, respectively.

To analyze wave propagation in such systems, Floquet-Bloch theory is employed [2], seeking solutions of the harmonic form

$$ \mathbf{u}_n(t) = \mathbf{\Phi} \lambda^n e^{i\omega t} \quad (2.2) $$

where $\lambda \in \mathbb{C}$ is the propagation constant, $\omega \in \mathbb{R}_+$ the angular frequency, and $\mathbf{\Phi} \in \mathbb{C}^d$ the mode shape. It is often convenient to set $\lambda = e^{-il\kappa}$, where $\kappa$ denotes the wavenumber.

The detailed derivation is provided in Appendix A.1, substituting (2.2) into (2.1) leads to the dispersion relation

$$ \left[ \sum_{p=1}^{P} \mathcal{K}_p(\lambda, \beta_p) - \omega^2 \mathbf{M} \right] \mathbf{\Phi} = 0 \quad (2.3) $$

with

$$ \mathcal{K}_p(\lambda, \beta_p) = \mathbf{A}_p^+(\beta_p) + \mathbf{B}_p^+(\beta_p)\lambda^p + \mathbf{A}_p^-(\beta_p) + \mathbf{B}_p^-(\beta_p)\lambda^{-p} \quad (2.4) $$

where $\mathbf{A}_p^+(\beta_p)$, $\mathbf{A}_p^-(\beta_p)$, $\mathbf{B}_p^+(\beta_p)$ and $\mathbf{B}_p^-(\beta_p)$ are square matrix of size $d \times d$. For a fixed $\lambda \in \mathbb{C}$ and set of parameters $\beta$, solving (2.3) yields $d$ admissible angular frequencies $\omega_j(\lambda, \beta)$, with $1 \le j \le d$.

This formal framework sets the stage for a systematic tuning strategy, where the parameters $\beta_p$ are chosen so as to enforce specific target points on the dispersion curves. The procedure is detailed in the next subsection.

## 2.2 Tuning method

Building upon the above framework, the aim is to parametrically control the dispersion curves of uniform nonlocal lattices. In this context, “parametric control” means prescribing specific points $(\lambda_i, \omega_i)$ that the dispersion relation of the model must satisfy. Formally, for each prescribed $\lambda_i$, one of the $d$ angular frequencies $\omega_j(\lambda_i, \beta)$ must satisfy

$$ \omega_j(\lambda_i, \beta) = \omega_i. $$

For the various parameterization objectives discussed in Section 4, it can be more convenient to operate directly on the dispersion curves, and thus to prescribe target points in the form $(\kappa_i, \omega_i)$. The procedure begins by selecting the maximum nonlocal order $P$ and fixing the interaction model between the nodes. The chosen model depends on $P$ physical parameters $\beta_p$ characterizing the nonlocal couplings. It is then possible to select $P$ target points $(\lambda_i, \omega_i)$ (or $(\kappa_i, \omega_i)$) to be interpolated by the dispersion relation. This leads to the nonlinear system

$$ \omega_j(\lambda_i, \beta)^2 - \omega_i^2 = 0, \quad i \in \{1, \dots, P\}, \quad (2.5) $$

with the $P$ unknowns $\beta = (\beta_p)$. In general, (2.5) is nonlinear in $\beta$. However, squaring the angular frequencies makes it possible, in certain cases, to obtain a linear system, for instance in the monoatomic chain with uniform nonlocal spring interactions. The applicability and versatility of the interpolation framework are best appreciated through concrete examples. Two representative cases are considered in this work: the fundamental monoatomic chain with uniform nonlocal spring interactions, presented in Appendix A.2, and the Euler-Bernoulli beam lattice, discussed in the next section, which serves to illustrate the method within a physically richer and structurally more complex context.

## 3 Application to an Euler-Bernoulli beam lattice

In this model, each node possesses two degrees of freedom ($d=2$), corresponding to transverse displacement and rotation, and is connected to neighboring nodes by Euler-Bernoulli beams. Each connecting beam is characterized by its Young’s modulus $E_p$ and second moment of inertia $I_p$, defining the stiffness parameter

$$ \beta_p = E_p I_p. $$
The nearest-neighbor configuration of such lattices has been studied by Madine and Colquitt [23], while uniform nonlocal extensions with identical beam properties for all interaction orders were introduced by Edge et al. [14]. The present formulation further generalizes this setting by allowing the mechanical properties of the connecting beams to depend explicitly on the interaction order p.

Following the formalism of Section 2 and as detailed in Section A.3.1, the dispersion relation of order-P lattices takes the form of Eq. (2.3), with explicit coefficient matrices

$$ \mathbf{A}_{p}^{+}(\beta) = \frac{\beta_{p}}{p^{3}} \begin{bmatrix} 12 & 6p \\ 6p & 4p^{2} \end{bmatrix}, \quad \mathbf{B}_{p}^{+}(\beta) = \frac{\beta_{p}}{p^{3}} \begin{bmatrix} -12 & 6p \\ -6p & 2p^{2} \end{bmatrix} \qquad (3.1) $$

$$ \mathbf{A}_{p}^{-}(\beta) = \frac{\beta_{p}}{p^{3}} \begin{bmatrix} 12 & -6p \\ -6p & 4p^{2} \end{bmatrix}, \quad \mathbf{B}_{p}^{-}(\beta) = \frac{\beta_{p}}{p^{3}} \begin{bmatrix} -12 & -6p \\ 6p & 2p^{2} \end{bmatrix} \qquad (3.2) $$

The corresponding lumped mass matrix is

$$ \mathbf{M} = \begin{bmatrix} 1.0 & 0.0 \\ 0.0 & \mu \end{bmatrix} \qquad (3.3) $$

with $0 < \mu \ll 1$. For each fixed wavenumber, the generalized eigenvalue problem (2.3) admits two distinct angular frequencies. In this work, only cases where $\Im(\omega) = 0$ are considered. The two positive roots of (2.3) are denoted by $\omega_{-} \le \omega_{+}$.

In Appendix A.3.2, analytical expressions for $\omega_{-}$ and $\omega_{+}$ are fully derived, along with a discussion on how these closed-form results can be used to simplify the interpolation system (2.5) for this model. However, such analytical formulations are not required for the present approach, and numerical computations are employed throughout this work.

The dispersion relation therefore consists of two branches. The interpolation method introduced earlier can be directly applied to this model. Interpolation points ($\lambda_{i}, \omega_{i}$) are selected with explicit specification of the branch at each point, leading to a system of the form (2.5). Depending on the design objective, either a single branch or both branches may be interpolated simultaneously.

As an illustrative example, whose result is shown in Figure 2a, consider the objective of flattening the lower branch $\omega_{-}$ over the interval $\kappa \in (1.5, \pi)$. This is achieved by selecting three target points:

$$ (\kappa_{1} = 1.5, \omega_{1} = 2.0), \quad (\kappa_{2} = 2.0, \omega_{2} = 2.0), \quad (\kappa_{3} = 2.5, \omega_{3} = 2.0), $$

leading to the nonlinear system

$$ \left\{ \begin{array}{l} \omega_{-}(\kappa_{1} = 1.5, \beta)^{2} - 2.0^{2} = 0 \\ \omega_{-}(\kappa_{2} = 2.0, \beta)^{2} - 2.0^{2} = 0 \\ \omega_{-}(\kappa_{3} = 2.5, \beta)^{2} - 2.0^{2} = 0 \end{array} \right. \qquad (3.4) $$

As explained in Appendix A.4.1, this system can be solved using numerical optimization methods, to which additional constraints, such as the positivity of the solution corresponding to the positivity of each parameter $\beta_{p}$, can be incorporated. Two additional examples of similar parameterizations are shown in Figure 2.

An important consequence of the nonlinear system (2.5) is the possible absence of solutions. In beam-interaction lattices, two main scenarios may lead to infeasibility: (i) attempting to impose excessively strong oscillations on a branch, especially at low frequencies; (ii) constraining both branches at the same wavenumber without ensuring sufficient separation between $\omega_{+}$ and $\omega_{-}$. In such cases, introducing additional nonlocal interactions increases the number of design parameters, thereby enlarging the feasible solution set. Finally, the analytical properties of the dispersion relation, established in Appendix A.3.2, provide useful validity criteria for admissible target curves, such as the vanishing derivatives of both branches at $\kappa = 0$ and $\kappa = \pi$.

In summary, the Euler-Bernoulli beam lattice provides a representative setting in which the interpolation framework can be explicitly implemented and tested. The model highlights both the potential and the limitations of the method, while offering practical criteria for feasible dispersion shaping. Building on this foundation, the next section explores how parametric interpolation can be employed to engineer a variety of distinctive wave phenomena.
<figure><img src="image_3.png" /><figcaption>Figure 2: Examples of dispersion curves for the Euler-Bernoulli beam lattice. In the first two cases, a single branch is interpolated, while in the third case both branches are simultaneously constrained (see Table 1).</figcaption></figure>

# 4 Engineering wave phenomena via parametric interpolation

The previous sections established the theoretical framework for uniform nonlocal lattices and introduced a systematic interpolation method for tailoring their dispersion relations. The Euler-Bernoulli beam lattice was then adopted as a representative model, providing explicit expressions for the interaction matrices and serving as a reference setting for practical implementation.

This section illustrates the potential of the proposed approach through a series of application examples. By prescribing selected points on the dispersion curves, distinct and physically relevant wave phenomena can be engineered, ranging from qualitative modifications of branch topology to quantitative control of wave-packet dynamics. Three representative cases are considered: the construction of rotons, the tuning of group velocity dispersion, and the design of evanescent modes within band gaps. Together, these examples demonstrate the versatility of parametric interpolation as a tool for engineering wave propagation in nonlocal lattices.

## 4.1 Constructing rotons via parametric dispersion control

In the Euler-Bernoulli beam lattice restricted to local interactions, both branches of the dispersion relation are injective over $\kappa \in [0, \pi]$, as can be observed from the analytical expressions derived in Appendix A.3.2. For a given angular frequency $\omega$, the Floquet-Bloch condition (2.2) admits a single propagation constant $\lambda \in [0, 1]$ in the purely local case. Introducing nonlocal interactions fundamentally changes this behavior: the dispersion branches may lose their injectivity, allowing multiple propagation constants $\lambda$ to correspond to the same frequency $\omega$.

The parametrization method enables the systematic construction of such non-injective dispersion curves. For a prescribed frequency $\omega_0$, several target points $(\lambda_i, \omega_0)$ with distinct propagation constants are selected. The resulting system (2.5) then yields a lattice model whose dispersion relation admits multiple propagation constants at the same angular frequency. The dispersion curve of such a system is presented in Figure 3b. When $\partial_{\kappa}\omega(\kappa_1) \neq \partial_{\kappa}\omega(\kappa_2)$, these solutions correspond to distinct group velocities associated with $\omega_0$.

The implications of this property can be highlighted through transient simulations. Consider a Gaussian wave packet imposed on the reference node $x_0$, centered around the frequency $\omega_0$:

$$u_0(t) = A \cos(\omega_0(t - t_0)) e^{-((t-t_0)/\tau)^2}. \quad (4.1)$$

As shown in Figure 3c, the initially localized excitation evolves into two distinct packets that propagate with
the two group velocities determined at $\omega_0$.

Loss of injectivity in the dispersion relation is accompanied by the appearance of local extrema along the branches. Such extrema are commonly referred to as rotons in the literature [5–8]. A roton corresponds to a vanishing group velocity at the associated frequency, leading to the localization of a wave packet around its emission point.

The interpolation method provides a direct strategy to enforce such conditions. Three interpolation points $(\lambda_{-1}, \omega_{-1})$, $(\lambda_0, \omega_0)$, and $(\lambda_{+1}, \omega_{+1})$ are chosen symmetrically around $\kappa_0$, with spacing $\Delta\kappa$:

$$ \omega_{-1} < \omega_0 < \omega_{+1}, \quad \kappa_{-1} < \kappa_0 < \kappa_{+1}, \quad \kappa_0 - \kappa_{-1} = \kappa_{+1} - \kappa_0 = \Delta\kappa. $$

A vanishing slope at $(\lambda_0, \omega_0)$ is then imposed through the central finite-difference criterion,

$$ \frac{\omega_{+1} - \omega_{-1}}{2\Delta\kappa} = 0. \qquad (4.2) $$

To ensure that the point indeed corresponds to a roton, the curvature must be nonzero, which is verified through the second-order finite difference:

$$ \omega_{-1} - 2\omega_0 + \omega_{+1} \neq 0, \qquad (4.3) $$

positive for a local minimum and negative for a local maximum. The dispersion curve of a model parametrized to exhibit a local minimum at $(\kappa = 2.0, \omega_+ = 8)$ is shown in Figure 3a.

<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Example of a roton created at the frequency $\omega_+ = 8$. Three interpolation points are selected to enforce a vanishing central finite difference at $\omega_+$ (see Table 2). (b) Model exhibiting two propagating solutions for the same frequency $\omega_- = 6$, obtained by prescribing two interpolation points at $(\kappa_1 = 1, \omega_-)$ and $(\kappa_2 = \pi, \omega_-)$. (c) Transient response associated with (b): the excitation at $\omega_- = 6$ splits into two wave packets with distinct group velocities.</figcaption></figure>

The emergence of rotons illustrates how parametric interpolation can modify the qualitative structure of dispersion relations. The same principle can be extended to tune more subtle features, such as the group velocity dispersion, which governs the spreading of wave packets over time.

## 4.2 Tailoring group velocity dispersion via parametric interpolation

In discrete Euler-Bernoulli beam lattices, the relation between angular frequency $\omega$ and wavenumber $\kappa$ is inherently nonlinear, reflecting their dispersive nature. As a consequence, a wave packet centered around a given frequency $\omega_0$ gradually spreads over time due to dispersion effects. This spreading is governed by the group velocity dispersion (GVD), defined as

$$ \mathbf{GVD}(\omega_0) = \partial_{\omega}^2 \kappa(\omega_0) \qquad (4.4) $$
By applying the local inversion theorem, the GVD can be expressed in terms of the curvature of the dispersion relation:

$$ \mathbf{GVD}(\omega_0) = -\frac{\partial^2 \omega(\kappa_0)}{c_0^2} \qquad (4.5) $$

where $c_0 = \partial_\kappa \omega(\kappa_0)$ denotes the group velocity at $\kappa_0$. This expression highlights the central role played by the second derivative of $\omega(\kappa)$ in controlling the spreading of wave packets.

As illustrated in Section 4.1, the interpolation method can be employed to prescribe not only the slope of the dispersion relation at selected points but also its local curvature. By selecting three interpolation points in the neighborhood of $(\kappa_0, \omega_0)$ and applying a finite-difference criterion, both the group velocity and the GVD can be tuned independently. It therefore becomes possible, for a given angular frequency $\omega_0$, to design lattices in which wave packets experience stronger or weaker spatial spreading over time. An example of such a parameterization is shown in Figure 4: after locally tuning the slope and curvature of the dispersion curves of three different models at the same frequency $\omega_0$, the corresponding wave packets centered around $\omega_0$ are observed to spread in space at different rates, depending on the imposed curvature.

<figure><img src="image_5.png" /><figcaption>Figure 4: (a) Dispersion curves of three parametrized models, each tailored in a frequency range around $\omega_0 = 12.0$. (b) Zoom on the customized region. Three interpolation points are selected to simultaneously control the slope and curvature of the dispersion relation near $\omega_0 = 12.0$ (see Table 3). For the three models, the group velocity is fixed at $c_0 = -2$, while the group velocity dispersion (GVD) takes values 0, 1.25, and 2.5, respectively. (c) Transient wave-packet analysis. The excitation imposed on node $x_0$ follows Eq. (4.1). The corresponding spectral content, shown in green in (b), is centered around $\omega_0$ and lies within the parametrized region of the dispersion curves.</figcaption></figure>

These results show that parametric interpolation not only alters the qualitative topology of dispersion branches, as in the case of rotons, but also provides a systematic means of adjusting quantitative properties such as the group velocity dispersion. This flexibility in shaping propagating-wave dynamics naturally extends to evanescent modes within band gaps, which can likewise be parametrically controlled, as discussed in the next subsection.

## 4.3 Designing evanescent wave decay via parametric interpolation

For a fixed angular frequency $\omega_0$, the interpolation method introduced earlier can be used to identify wave solutions associated with propagation constants $\lambda$ of unit modulus. Such solutions correspond to propagating waves that maintain constant amplitude in the absence of dissipative effects. However, even without material dissipation, the periodicity of the lattice gives rise to frequency intervals in which admissible propagation constants satisfy $|\lambda| \neq 1$. The corresponding modes, referred to as evanescent waves, exhibit exponential spatial decay along a preferred direction.
The propagation constant can be expressed in polar form as

$$ \lambda = |\lambda|e^{-i\kappa_R}, \qquad (4.6) $$

or equivalently,

$$ \lambda = e^{-i\kappa}, \qquad (4.7) $$

where $\kappa = \kappa_R + i\kappa_I$ and $\kappa_I = \ln(|\lambda|)$. When $|\lambda| < 1$, the amplitude decays to the right with factor $e^{n\kappa_I}$. Moreover, if $\lambda$ is a solution at frequency $\omega_0$, then $1/\lambda$ is also a solution, so the analysis may, without loss of generality, be restricted to $|\lambda| < 1$.

The proposed interpolation framework naturally extends to such evanescent solutions by solving the dispersion relation in the complex plane. For the Euler-Bernoulli beam lattice, the dispersion relation exhibits two propagating branches separated by a band gap. Within this interval, the only admissible solutions correspond to evanescent modes. The band-gap boundaries [$\omega_0, \omega_1$] can be determined by selecting two target points, ($\lambda_0 = -1, \omega_0$) and ($\lambda_1 = -1, \omega_1$). Figure 5 shows the dispersion curves of two models parametrized using these two target points so as to produce the same band gap over [$\omega_1, \omega_2$]. A third target point, which differs between the two models, ($\lambda_3, \omega_3$) with $|\lambda_3| < 1$ and $\omega_3 \in [\omega_1, \omega_2]$, is introduced to construct a parametrized evanescent solution within the band gap. The modulus of $\lambda_3$ directly controls the rate of exponential decay of the corresponding mode.

<figure><img src="image_6.png" /><figcaption>Figure 5: Two uniform nonlocal beam-lattice models exhibiting the same band gap [6.0, 12.0]. The colormap represents the modulus of $\lambda$, providing a complementary visualization of $\Im(\kappa)$. An evanescent wave solution at $\omega = 6.5$ is parametrized with different decay rates: $\kappa_I = -0.2$ in the first case and $\kappa_I = -0.7$ in the second. Transient analysis confirms that the spatial attenuation follows approximately $e^{\kappa_I}$. The excitation applied to the reference node $x_0$ follows Eq. (4.1). The parameters used are listed in Table 4.</figcaption></figure>

This example shows that evanescent modes can be actively designed and tailored within the proposed framework. Such control over spatial decay complements the previous results on rotons and group velocity dispersion, thereby further expanding the range of achievable wave-propagation behaviors in nonlocal lattices.
# 5 Conclusion

This work introduced a systematic interpolation-based method for tailoring dispersion relations in uniform nonlocal lattices. By prescribing selected frequency-wavenumber points as interpolation constraints, the method provides localized and flexible customization without requiring closed-form dispersion expressions or global curve-fitting procedures. The approach was developed in a general theoretical setting and illustrated through the Euler-Bernoulli beam lattice, a representative model capturing both bending and rotational dynamics.

Several applications were presented to highlight the range of achievable phenomena. The construction of roton-like extrema demonstrated how qualitative topological features of dispersion branches can be engineered. The adjustment of group velocity dispersion showed that quantitative wave-packet spreading can be controlled through local curvature tuning. Finally, the parametrization of evanescent modes within band gaps illustrated that exponential decay rates can likewise be designed, extending the scope of dispersion control beyond propagating waves.

Together, these examples establish parametric interpolation as a versatile and computationally tractable tool for dispersion engineering in nonlocal lattices. Future work may extend the framework to higher-dimensional systems, explore the inclusion of damping and dissipation mechanisms, or integrate the method with optimization schemes for application-driven design. Such developments would further enhance the applicability of nonlocal lattices in advanced metamaterial concepts for vibration isolation, acoustic control, and wave-based technologies.

# Acknowledgements

The authors acknowledge the French National Research Agency (ANR) for its financial support through the JCJC-Archi-Noise project (grant no. ANR-23-CE51-0045).
# A Appendix

## A.1 General formulation of the dispersion relation in uniform nonlocal models

The derivation of the dispersion relation for uniform nonlocal lattices, presented here, provides the theoretical foundation for the tuning procedure introduced in Section 2.2. It formalises the connection between the physical parameters $\beta$ and the spectral properties ($\lambda, \omega$) of the system, thus enabling the targeted adjustment of dispersion curves.

Consider a uniform nonlocal lattice of order $P$. Such a system is described by specifying the interactions between its constituent nodes, where each interaction depends solely on its nonlocal order $p \in \{1, \dots, P\}$.

Let $x_l$ and $x_r$ denote two nodes separated by a distance $pl$. The equations of motion for this pair can be expressed in block-matrix form as

$$ \begin{bmatrix} \mathbf{M}_p^+ & \mathbf{O}_p^+ \\ \mathbf{O}_p^- & \mathbf{M}_p^- \end{bmatrix} \partial_t^2 \begin{bmatrix} \mathbf{u}_l \\ \mathbf{u}_r \end{bmatrix} + \begin{bmatrix} \mathbf{A}_p^+ & \mathbf{B}_p^+ \\ \mathbf{B}_p^- & \mathbf{A}_p^- \end{bmatrix} \begin{bmatrix} \mathbf{u}_l \\ \mathbf{u}_r \end{bmatrix} = \mathbf{0} \qquad (A.1) $$

where the matrices $\mathbf{M}_p^\pm$, $\mathbf{O}_p^\pm$, $\mathbf{A}_p^\pm$, and $\mathbf{B}_p^\pm$ are $d \times d$ blocks determined by the physical parameter $\beta_p$. By assembling two such elements, the linear dependence of a given node $x_n$ on its $p$-th order neighbors $x_{n-p}$ and $x_{n+p}$ can be expressed as

$$ \begin{bmatrix} \mathbf{O}_p^- & \mathbf{M}_p^- + \mathbf{M}_p^+ & \mathbf{O}_p^+ \end{bmatrix} \partial_t^2 \begin{bmatrix} \mathbf{u}_{n-p} \\ \mathbf{u}_n \\ \mathbf{u}_{n+p} \end{bmatrix} + \begin{bmatrix} \mathbf{B}_p^- & \mathbf{A}_p^- + \mathbf{A}_p^+ & \mathbf{B}_p^+ \end{bmatrix} \begin{bmatrix} \mathbf{u}_{n-p} \\ \mathbf{u}_n \\ \mathbf{u}_{n+p} \end{bmatrix} = \mathbf{0} \qquad (A.2) $$

Continuing this assembly process for all $p \le P$ yields the stiffness matrix $\mathbf{K}$ from Eq. (2.1) in the compact form

$$ \mathbf{K} = \begin{bmatrix} \mathbf{B}_p^- & \mathbf{B}_{p-1}^- & \dots & \mathbf{A} & \dots & \mathbf{B}_{p-1}^+ & \mathbf{B}_p^+ \end{bmatrix} \qquad (A.3) $$

where

$$ \mathbf{A} = \sum_{p=1}^{P} (\mathbf{A}_p^- + \mathbf{A}_p^+) \qquad (A.4) $$

A similar construction applies to the mass matrix $\mathbf{M}$. In the present work, a lumped-mass assumption is adopted, namely $\mathbf{O}_p^\pm = \mathbf{0}$ and $\mathbf{M}_p^- = \mathbf{M}_p^+$ are diagonal.

Seeking a Floquet-Bloch wave solution of Eq. (2.1) in the form of Eq. (2.2) gives

$$ -\omega^2 \mathbf{M} \mathbf{\Phi} u_n + \mathbf{K} [\lambda^{-P} \mathbf{\Phi} \quad \dots \quad \mathbf{\Phi} \quad \dots \quad \lambda^{-P} \mathbf{\Phi}]^T u_n = 0 \qquad (A.5) $$

This directly leads to the dispersion relation given in Eq. (2.3).

## A.2 Monoatomic chain with uniform nonlocal interactions

This subsection considers a particular case of the uniform nonlocal lattice model. The system consists of a one-dimensional periodic chain of identical nodes, each having a single degree of freedom and mass $m$, separated by a constant lattice spacing $l$.

The mechanical interactions between nodes are purely axial and described by springs. Nearest-neighbor couplings ($p=1$) have stiffness $\beta_1$, while higher-order nonlocal interactions connect each node to its $p$-th neighbors ($p=2, \dots, P$) on both sides, with stiffness $\beta_p$.

Under these assumptions, the general dispersion relation

$$ m \omega^2 = 2 \sum_{p=1}^{P} \beta_p (1 - \cos(p\kappa)) \qquad (A.6) $$
Using the trigonometric identity $1 - \cos(p\kappa) = 2\sin^2(p\kappa/2)$, Eq. (A.6) becomes

$$m\omega^2(\kappa) = 4 \sum_{p=1}^{P} \beta_p \sin^2(p\kappa/2). \qquad (A.7)$$

Following the tuning approach presented in Section 2.2, a set of P wave numbers is selected:

$$\mathbf{k} = [\kappa_1, \dots, \kappa_P], \qquad (A.8)$$

together with the corresponding target squared angular frequencies:

$$\omega^2 = [\omega_1^2, \dots, \omega_P^2]. \qquad (A.9)$$

Substituting these values into Eq. (A.7) leads to a linear system in the unknown stiffnesses $\beta_p$:

$$\mathbf{A} \beta = m \omega^2, \qquad (A.10)$$

where $\beta = [\beta_1, \dots, \beta_P]^T$ and the matrix $\mathbf{A} \in \mathbb{R}^{P \times P}$ is defined by

$$\mathbf{A} = 4 \begin{bmatrix} \sin^2(\kappa_1/2) & \dots & \sin^2(P\kappa_1/2) \\ \vdots & \ddots & \vdots \\ \sin^2(\kappa_P/2) & \dots & \sin^2(P\kappa_P/2) \end{bmatrix}. \qquad (A.11)$$

The solution of Eq. (A.10) provides the stiffness values $\beta_p$ ensuring that the dispersion curve passes exactly through the prescribed $(\kappa_i, \omega_i)$ points. To guarantee physical admissibility ($\beta_p \ge 0$), the system can be solved using a non-negative least squares (NNLS) approach.

## A.3 Euler-Bernoulli beam lattices

### A.3.1 Dispersion relation for the Euler-Bernoulli beam lattice

This appendix derives the equations of motion for the Euler-Bernoulli beam lattice model introduced in Section 3. The starting point is the governing equation for a single Euler-Bernoulli beam, as given in [24]:

$$EI \frac{\partial^4 w}{\partial x^4} + \rho A \frac{\partial^2 w}{\partial t^2} = 0 \qquad (A.12)$$

where $E$ denotes the Young's modulus, $I$ the second moment of inertia, $\rho$ the density, and $A$ the cross-sectional area.

Consider an element of this beam of length $p$, discretized using standard finite element theory. Each end of the element possesses two degrees of freedom: the transverse displacement $w$ and the cross-section rotation $\theta$. Let $\mathbf{u}_l$ and $\mathbf{u}_r$ denote the nodal displacement-rotation vectors at the left and right ends, respectively. The discrete equation of motion reads:

$$\mathbf{M}_e \partial_t^2 \begin{bmatrix} \mathbf{u}_l \\ \mathbf{u}_r \end{bmatrix} + \mathbf{K}_e \begin{bmatrix} \mathbf{u}_l \\ \mathbf{u}_r \end{bmatrix} \qquad (A.13)$$

with the stiffness matrix

$$\mathbf{K}_e = \frac{EI}{p^3} \begin{bmatrix} 12 & 6p & -12 & 6p \\ 6p & 4p^2 & -6p & 2p^2 \\ -12 & -6p & 12 & -6p \\ 6p & 2p^2 & -6p & 4p^2 \end{bmatrix} \qquad (A.14)$$

and the consistent mass matrix

$$\mathbf{M}_e = \frac{\rho A p}{420} \begin{bmatrix} 156 & 22p & 54 & -13p \\ 22p & 4p^2 & 13p & -3p^2 \\ 54 & 13p & 156 & -22p \\ -13p & -3p^2 & -22p & 4p^2 \end{bmatrix} \qquad (A.15)$$
In the present work, masses are concentrated at the lattice interfaces. Following the approach of Madine and Colquitt [23], a lumped mass matrix is adopted:

$$ M_e = \text{diag}[1.0 \ \mu \ 1.0 \ \mu] \qquad (A.16) $$

with $0 < \mu \ll 1$ a small rotational-to-translational inertia ratio.

Equation (A.13) matches the general interaction form (A.1). By direct identification, the interaction matrices for the Euler-Bernoulli beam lattice are obtained as those given in Section 3. This establishes the explicit link between the beam element formulation and the dispersion relation used in the main text.

### A.3.2 Analytical dispersion relation for the Euler-Bernoulli beam lattice

In this appendix, we derive an explicit analytical form of the dispersion relation for the Euler-Bernoulli beam lattice. This allows us to write the interpolation system (2.5) in closed form for this specific case.

We start from the definitions of the interaction matrices $A_p^-, A_p^+, B_p^-$ and $B_p^+$ given in Section 3. The sum of the stiffness contributions reads:

$$ \sum_{p=1}^{P} \mathcal{K}_p(\lambda, \beta_p) = \sum_{p=1}^{P} \frac{\beta_p}{p^3} \left( \begin{bmatrix} 24 & 0 \\ 0 & 8 \end{bmatrix} + \begin{bmatrix} -12(\lambda^p + \lambda^{-p}) & 6p(\lambda^p - \lambda^{-p}) \\ -6p(\lambda^p - \lambda^{-p}) & 2p^2(\lambda^p + \lambda^{-p}) \end{bmatrix} \right) \qquad (A.17) $$

Using Euler's identity $\lambda = e^{i\kappa}$, the above expression becomes:

$$ \sum_{p=1}^{P} \mathcal{K}_p(\lambda, \beta_p) = \begin{bmatrix} \alpha & \gamma \\ -\gamma & \delta \end{bmatrix} \qquad (A.18) $$

where the scalar coefficients are given by:

$$ \alpha = 24 \sum_{p=1}^{P} \beta_p [1 - \cos(\kappa p)]/p^3 \qquad (A.19) $$

$$ \delta = 4 \sum_{p=1}^{P} \beta_p [2 + \cos(\kappa p)]/p \qquad (A.20) $$

$$ \gamma = 12i \sum_{p=1}^{P} \beta_p \sin(\kappa p)/p^2 \qquad (A.21) $$

The dispersion relation follows from the determinant condition:

$$ \sum_{p=1}^{P} \mathcal{K}_p(\lambda, \beta_p) - \omega^2 \mathbf{M} \qquad (A.22) $$

where the mass matrix $M$ is defined in Section 3. A straightforward computation yields:

$$ (\alpha\delta + \gamma^2) - (\mu\alpha + \delta)\omega^2 + \mu\omega^4 = 0 \qquad (A.23) $$

This quartic in $\omega$ reduces to a quadratic in $\omega^2$, whose two positive roots $\omega_\pm \le \omega_+$ are:

$$ \omega_{\pm}^2(\lambda, \beta) = \left[ \mu\alpha + \delta \pm \sqrt{(\mu\alpha - \delta)^2 - 4\mu\gamma^2} \right] / 2\mu \qquad (A.24) $$

In particular, in the case where $P = 1$, we can notice that the curves $\omega_\pm$ are injective on $[0, \pi]$ as a sum of injective curves. It is useful to note that $\omega_\pm^2$ can be decomposed into a linear and a nonlinear part with respect to the stiffness coefficients $\beta$:

$$ \omega_{\pm}^2(\lambda, \beta) = \mathcal{L}(\lambda, \beta) \mp G(\lambda, \beta) \qquad (A.25) $$
with

$$ \mathcal{L}(\lambda, \beta) := (\mu\alpha + \delta)/2\mu, \qquad (A.26) $$

$$ G(\lambda, \beta) := \sqrt{(\mu\alpha - \delta)^2 - 4\mu\gamma^2}/2\mu. \qquad (A.27) $$

If, for a fixed $\lambda_i$, one prescribes $\omega_-(\lambda_i, \beta) = a$ and $\omega_+(\lambda_i, \beta) = b$, the interpolation system (2.5) simplifies to:

$$ \left\{ \begin{array}{l} \mathcal{L}(\lambda_i, \beta) = (b+a)/2 \\ G(\lambda, \beta) = (b-a)/2 \end{array} \right. \qquad (A.28) $$

## A.4 Calculation of stiffnesses

### A.4.1 Positivity constraint on stiffness parameters

In the case of beam-type interactions, the interpolation system (2.5) is nonlinear. To solve it, iterative numerical schemes of the Newton type are employed. Throughout this work, an important design constraint has been systematically enforced: all presented models are obtained using strictly positive stiffness values. This choice ensures the physical realism of the constructed lattices and avoids configurations that would correspond to unstable or non-physical mechanical responses.

Accordingly, a positivity constraint is incorporated into the numerical solvers. In contrast, for the special case of nonlocal spring-interaction lattices discussed in Section A.2, the system (2.5) is linear. In such cases, a direct solution almost always exists for any set of prescribed target points, yielding stiffness values that reproduce the desired dispersion curve. However, negative stiffnesses may arise in these unconstrained solutions. To prevent this, the system is preferably solved using a constrained optimization approach, such as the Non-Negative Least Squares (NNLS) method developed by Lawson and Hanson [25]. This formulation guarantees physically admissible, strictly positive stiffness distributions while maintaining accurate interpolation of the prescribed dispersion points.

### A.4.2 Approximation of target dispersion profiles

In some cases, it may be desirable to tune a nonlocal lattice model so that its dispersion profile approximates a predefined target curve. To illustrate this, the interpolation method is applied here to nonlocal spring-interaction models. This scenario has been investigated in [13], three target dispersion profiles from that study are revisited using the present approach.

Without enforcing stiffness positivity, the proposed interpolation method yields models whose dispersion curves closely match the target profiles. However, when the positivity constraint on the stiffnesses is imposed, globally accurate reconstructions become infeasible. In such situations, it is necessary to restrict attention to a specific region of interest along the dispersion curve and focus the interpolation within that interval.

Initially, interpolation points are chosen empirically to provide a reasonable match to the overall target, denoted as case **WN** (without constraint). Although this approach effectively reproduces the qualitative shape of the curve, it may yield negative stiffness values. To address this issue, the same set of interpolation points is used within a Non-Negative Least Squares (NNLS) formulation to enforce stiffness positivity, yielding the case **OP 1**. This procedure guarantees physical admissibility but generally reduces the global accuracy of the reconstruction.

To improve the trade-off between fidelity and feasibility, a localized optimization strategy is introduced. A region of interest is selected on the target curve, and new interpolation points are chosen within that interval. The NNLS regression is then reapplied, yielding the case **OP 2**. This local parametrization improves the agreement between the model and the target curve in the selected frequency range, while maintaining positive stiffnesses throughout.

The different results obtained under the three interpolation strategies, unconstrained (WN), globally constrained (OP 1), and locally constrained (OP 2), are summarized in Table 3 and illustrated in Figure 6, highlighting the trade-off between accuracy and physical admissibility.
<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison between interpolation strategies. The target dispersion curve is shown in black. Red markers correspond to global interpolation points (WN), while blue markers indicate locally optimized points (OP 2) within the selected region of interest.</figcaption></figure>
### A.4.3 Numerical parameters used in the figures

The following tables summarize the numerical parameters employed in the simulations associated with each figure. They include the prescribed interpolation points ($\kappa_i, \omega_i$) and the corresponding optimized stiffness values $\beta_p$.

<table><caption>Table 1: Examples of parameterization.</caption><thead><tr><th colspan="9">Parameters of figure 2</th></tr><tr><th colspan="3">(a)</th><th colspan="3">(b)</th></tr></thead><tbody><tr><td>$\kappa_1 = 1.5$</td><td>$\omega_{1-} = 2.0$</td><td>$\beta_1 = 0.0668$</td><td>$\kappa_1 = 1.0$</td><td>$\omega_{1+} = 6.0$</td><td>$\beta_1 = 0.1078$</td></tr><tr><td>$\kappa_2 = 2.0$</td><td>$\omega_{2-} = 2.0$</td><td>$\beta_2 = 0.3490$</td><td>$\kappa_2 = 1.5$</td><td>$\omega_{2+} = 6.0$</td><td>$\beta_2 = 0.1740$</td></tr><tr><td>$\kappa_3 = 3.0$</td><td>$\omega_{3-} = 2.0$</td><td>$\beta_3 = 0.4526$</td><td>$\kappa_3 = 2.0$</td><td>$\omega_{3+} = 6.0$</td><td>$\beta_3 = 0.1780$</td></tr><tr><td></td><td></td><td></td><td>$\kappa_4 = 3.0$</td><td>$\omega_{4+} = 6.0$</td><td>$\beta_4 = 0.1533$</td></tr><tr><th colspan="3">(c)</th><th colspan="3"></th></tr><tr><td>$\kappa_1 = 1.5$</td><td>$\omega_{1-} = 1.0$</td><td>$\beta_1 = 0.0116$</td><td></td><td></td><td></td></tr><tr><td>$\kappa_2 = 2.0$</td><td>$\omega_{2-} = 1.0$</td><td>$\beta_2 = 0.0768$</td><td></td><td></td><td></td></tr><tr><td>$\kappa_3 = 3.0$</td><td>$\omega_{3-} = 1.0$</td><td>$\beta_3 = 0.1790$</td><td></td><td></td><td></td></tr><tr><td>$\kappa_4 = 1.5$</td><td>$\omega_{4+} = 6.0$</td><td>$\beta_4 = 0.2347$</td><td></td><td></td><td></td></tr><tr><td>$\kappa_5 = 2.0$</td><td>$\omega_{5+} = 6.0$</td><td>$\beta_5 = 0.3316$</td><td></td><td></td><td></td></tr><tr><td>$\kappa_6 = 3.0$</td><td>$\omega_{6+} = 6.0$</td><td>$\beta_6 = 0.2269$</td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 2: Parameterization of rotons.</caption><thead><tr><th colspan="9">Parameters of figure 3</th></tr><tr><th colspan="3">(a)</th><th colspan="3">(b)</th></tr></thead><tbody><tr><td>$\kappa_1 = 1.8$</td><td>$\omega_{1+} = 7.9$</td><td>$\beta_1 = 0.2337$</td><td>$\kappa_1 = 1.0$</td><td>$\omega_{1-} = 6.0$</td><td>$\beta_1 = 0.75$</td></tr><tr><td>$\kappa_2 = 2.0$</td><td>$\omega_{2+} = 8.0$</td><td>$\beta_2 = 0.0908$</td><td>$\kappa_2 = \pi$</td><td>$\omega_{2-} = 6.0$</td><td>$\beta_2 = 19.3979$</td></tr><tr><td>$\kappa_3 = 2.2$</td><td>$\omega_{3+} = 7.9$</td><td>$\beta_3 = 0.5148$</td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 3: Parameterization of the dispersion velocity group.</caption><thead><tr><th colspan="9">Parameters of figure 4</th></tr><tr><th colspan="3">(a)</th><th colspan="3">(b)</th></tr></thead><tbody><tr><td>$\kappa_1 = 1.3$</td><td>$\omega_{1+} = 12.4$</td><td>$\beta_1 = 0.8601$</td><td>$\kappa_1 = 1.3$</td><td>$\omega_{1+} = 12.6$</td><td>$\beta_1 = 0.4744$</td></tr><tr><td>$\kappa_2 = 1.5$</td><td>$\omega_{2+} = 12.0$</td><td>$\beta_2 = 0.2290$</td><td>$\kappa_2 = 1.5$</td><td>$\omega_{2+} = 12.0$</td><td>$\beta_2 = 1.8977$</td></tr><tr><td>$\kappa_3 = 1.7$</td><td>$\omega_{3+} = 11.6$</td><td>$\beta_3 = 0.0821$</td><td>$\kappa_3 = 1.7$</td><td>$\omega_{3+} = 11.8$</td><td>$\beta_3 = 0.2367$</td></tr><tr><th colspan="3">(c)</th><th colspan="3"></th></tr><tr><td>$\kappa_1 = 1.3$</td><td>$\omega_{1+} = 12.8$</td><td>$\beta_1 = 0.1937$</td><td></td><td></td><td></td></tr><tr><td>$\kappa_2 = 1.5$</td><td>$\omega_{2+} = 12.0$</td><td>$\beta_2 = 3.2726$</td><td></td><td></td><td></td></tr><tr><td>$\kappa_3 = 1.7$</td><td>$\omega_{3+} = 12.0$</td><td>$\beta_3 = 0.1344$</td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>Table 4: Parameterization of evanescent waves.</caption><caption>Parameters of figure 5</caption><thead><tr><th colspan="3">(a)</th><th colspan="3">(b)</th></tr></thead><tbody><tr><td>κ₁ = π</td><td>ω₁ = 6.0</td><td>β₁ = 0.5855</td><td>κ₁ = 1.0</td><td>ω₁ = 6.0</td><td>β₁ = 0.7330</td></tr><tr><td>κ₂ = π</td><td>ω₂ = 12.0</td><td>β₂ = 0.0625</td><td>κ₂ = π</td><td>ω₂ = 12.0</td><td>β₂ = 0.8492</td></tr><tr><td>κ₃ = π - 0.2i</td><td>ω₃ = 6.5</td><td>β₃ = 4.4422</td><td>κ₃ = π - 0.7i</td><td>ω₃ = 6.5</td><td>β₃ = 0.4595</td></tr></tbody></table>

<table><caption>Table 5: Stiffness coefficients $β_p$ obtained for three target dispersion curves (Figures 6), using three strategies: unconstrained interpolation (WN), non-negative least squares applied to the same points (OP 1), and localized constrained interpolation within a selected region of interest (OP 2).</caption><caption>Parameters of figure 6</caption><thead><tr><th rowspan="2">β<sub>p</sub></th><th colspan="3">(a)</th><th colspan="3">(b)</th><th colspan="3">(c)</th></tr><tr><th>WN</th><th>OP 1</th><th>OP 2</th><th>WN</th><th>OP 1</th><th>OP 2</th><th>WN</th><th>OP 1</th><th>OP 2</th></tr></thead><tbody><tr><td>β<sub>1</sub></td><td>31.9336</td><td>30.0132</td><td>4.0398</td><td>-0.4055</td><td>0</td><td>0</td><td>1.3113</td><td>1.2185</td><td>1.2185</td></tr><tr><td>β<sub>2</sub></td><td>15.1452</td><td>12.2734</td><td>0</td><td>0.8525</td><td>0.8019</td><td>0.5316</td><td>-0.2029</td><td>0.0551</td><td>0.0551</td></tr><tr><td>β<sub>3</sub></td><td>0.2799</td><td>0</td><td>0</td><td>0.6320</td><td>0.4303</td><td>0.7955</td><td>-0.1634</td><td>0</td><td>0</td></tr><tr><td>β<sub>4</sub></td><td>-4.3528</td><td>0</td><td>6.1786</td><td>0.0495</td><td>0</td><td>0</td><td>0.0900</td><td>0.0417</td><td>0.0417</td></tr><tr><td>β<sub>5</sub></td><td>-1.2462</td><td>0</td><td>13.4806</td><td>-0.0463</td><td>0</td><td>0</td><td>0.0843</td><td>0.0381</td><td>0.0381</td></tr><tr><td>β<sub>6</sub></td><td>3.0751</td><td>0.2038</td><td>16.2458</td><td>0.0446</td><td>0</td><td>0.0704</td><td>0.0067</td><td>0</td><td>0</td></tr><tr><td>β<sub>7</sub></td><td>1.7682</td><td>0</td><td>11.2344</td><td>0.0328</td><td>0</td><td>0</td><td>0.0190</td><td>0</td><td>0</td></tr><tr><td>β<sub>8</sub></td><td>-6.0159</td><td>0</td><td>0</td><td>-0.0148</td><td>0</td><td>0</td><td>0.0108</td><td>0</td><td>0</td></tr><tr><td>β<sub>9</sub></td><td>4.8657</td><td>2.9468</td><td>8.6173</td><td>0</td><td>0</td><td>0</td><td>0.0041</td><td>0</td><td>0</td></tr><tr><td>β<sub>10</sub></td><td>-1.2064</td><td>0</td><td>0.7174</td><td>0</td><td>0</td><td>0</td><td>0.0057</td><td>0</td><td>0</td></tr><tr><td>β<sub>11</sub></td><td>-0.5215</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β<sub>12</sub></td><td>0.0741</td><td>0</td><td>0</td><td>1.2250</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β<sub>13</sub></td><td>0.3604</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β<sub>14</sub></td><td>-0.0839</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β<sub>15</sub></td><td>-0.0908</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β<sub>16</sub></td><td>-0.1176</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β<sub>17</sub></td><td>0.2524</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β<sub>18</sub></td><td>-0.1473</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β<sub>19</sub></td><td>0.0266</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β<sub>20</sub></td><td>0.0026</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>
References

[1] Leon Brillouin. *Wave propagation in periodic structures*. MGH, 1946.

[2] Mahmoud I Hussein, Michael J Leamy, and Massimo Ruzzene. Dynamics of phononic materials and structures: Historical origins, recent progress, and future outlook. *Applied Mechanics Reviews*, 66(4): 040802, 2014.

[3] A Carcaterra, F Coppo, F Mezzani, and S Pensalfini. Long-range retarded elastic metamaterials: Wave-stopping, negative, and hypersonic or superluminal group velocity. *Physical Review Applied*, 11(1):014041, 2019.

[4] Amirsajjad Rezaei, Federica Mezzani, and Antonio Carcaterra. Wave propagation with long-range forces and mistuning effects. *Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science*, 235(14):2622-2628, 2021.

[5] Julio Andrés Iglesias Martínez, Michael Fidelis Groß, Yi Chen, Tobias Frenzel, Vincent Laude, Muamer Kadic, and Martin Wegener. Experimental observation of roton-like dispersion relations in metamaterials. *Science advances*, 7(49):eabm2189, 2021.

[6] Ke Wang, Yi Chen, Muamer Kadic, Changguo Wang, and Martin Wegener. Nonlocal interaction engineering of 2d roton-like dispersion relations in acoustic and mechanical metamaterials. *Communications Materials*, 3(1):35, 2022.

[7] Zhenxiao Zhu, Zhen Gao, Gui-Geng Liu, Yong Ge, Yin Wang, Xiang Xi, Bei Yan, Fujia Chen, Perry Ping Shum, Hong-xiang Sun, et al. Observation of multiple rotons and multidirectional roton-like dispersion relations in acoustic metamaterials. *New Journal of Physics*, 24(12):123019, 2022.

[8] Luca Iorio, Jacopo Maria De Ponti, Federico Maspero, and Raffaele Ardito. Roton-like dispersion via polarisation change for elastic wave energy control in graded delay-lines. *Journal of Sound and Vibration*, 572:118167, 2024.

[9] Sharat Paul, Md Nahid Hasan, and Pai Wang. Nonlocality-enabled inverse design of dirac-type and higher-order degeneracies for traveling and evanescent waves in phononic crystals. *arXiv preprint arXiv:2503.21020*, 2025.

[10] Yi Chen, Jonathan LG Schneider, Ke Wang, Philip Scott, Sebastian Kalt, Muamer Kadic, and Martin Wegener. Anomalous frozen evanescent phonons. *Nature Communications*, 15(1):8882, 2024.

[11] Ke Wang, Jonathan LG Schneider, Yi Chen, Julio Andrés Iglesias Martínez, Muamer Kadic, and Changguo Wang. Roton dispersion and displacement oscillatory for a nonlocal mechanical system. *European Journal of Mechanics-A/Solids*, page 105842, 2025.

[12] Arnab Banerjee, Kamal K Bera, and Sondipon Adhikari. Enhanced damping and bandwidth in roton-like dispersion of a beyond nearest neighbor periodic chain. *Physics Letters A*, 524:129821, 2024.

[13] Arash Kazemi, Kshiteej J. Deshmukh, Fei Chen, Yunya Liu, Bolei Deng, Henry Chien Fu, and Pai Wang. Drawing dispersion curves: Band structure customization via nonlocal phononic crystals. *Phys. Rev. Lett.*, 131:176101, 2023.

[14] Robyn Edge, Erica Paul, Katie H Madine, Daniel J Colquitt, Timonthy Starkey, and Gregory James Chaplain. Discrete euler-bernoulli beam lattices with beyond nearest connections. *New Journal of Physics*, 2025.

[15] Flavia Guarracino, Massimiliano Fraldi, and Nicola M Pugno. Local-to-non-local transition laws for stiffness-tuneable monoatomic chains preserving springs mass. *Philosophical Transactions A*, 382(2279): 20240037, 2024.

[16] Yi Chen, Ke Wang, Muamer Kadic, Sebastien Guenneau, Changguo Wang, and Martin Wegener. Phonon transmission through a nonlocal metamaterial slab. *Communications Physics*, 6(1):75, 2023.
[17] RG Edge, SAR Horsley, TA Starkey, and GJ Chaplain. Spring-loaded: Engineering dispersion relations using complex couplings. *arXiv preprint arXiv:2504.03539*, 2025.

[18] F Ongaro, PH Beoletto, F Bosia, Marco Miniaci, and NM Pugno. Closed-form solutions for wave propagation in hexagonal diatomic non-local lattices. *International Journal of Mechanical Sciences*, 291:110095, 2025.

[19] Esmaeal Ghavanloo and S Ahmad Fazelzadeh. Wave propagation in one-dimensional infinite acoustic metamaterials with long-range interactions. *Acta Mechanica*, 230(12):4453-4461, 2019.

[20] Yi Chen, Mahmoud AA Abouelatta, Ke Wang, Muamer Kadic, and Martin Wegener. Nonlocal cable-network metamaterials. *Advanced Materials*, 35(15):2209988, 2023.

[21] GJ Chaplain, IR Hooper, AP Hibbins, and TA Starkey. Reconfigurable elastic metamaterials: Engineering dispersion with beyond nearest neighbors. *Physical Review Applied*, 19(4):044061, 2023.

[22] Sharat Paul, Md Nahid Hasan, Henry Chien Fu, and Pai Wang. Complete inverse design to customize two-dimensional dispersion relation via nonlocal phononic crystals. *Physical Review B*, 110(14):144304, 2024.

[23] KH Madine and DJ Colquitt. Dynamic green's functions in discrete flexural systems. *The Quarterly Journal of Mechanics and Applied Mathematics*, 74(3):323-350, 2021.

[24] Karl F Graff. *Wave motion in elastic solids*. Courier Corporation, 2012.

[25] Charles L Lawson and Richard J Hanson. *Solving least squares problems*. SIAM, 1995.