## 2 Modal Analysis and Bandgap Formation in Mechanical Metastructures
This technique, widely used in modal analysis, produces convergent solutions to the boundary value problem as formulated. Using modal decomposition, the beam's deflection in the domain \(D\) is expressed as a sum of modal shapes in one direction. This assumes that the behavior of the beam can be accurately represented by a finite number of modes, as expressed in Equation (12):  
\[w(x,t)=\sum_{m=1}^{N_{m}}\phi_{m}(x)z_{m}(t), \tag{12}\]  
Here, \(\phi_{m}(x)\) denotes the spatial mode shape, and \(z_{m}(t)\) is the time-dependent modal coordinate for the \(m\)-th mode. These modal representations are crucial in modeling the dynamics of a flexible beam with integrated discrete resonators. Incorporating the modal expansion from Equation (12) into the system's governing differential equation, given by Equation (1), yields Equation (13). This resultant equation effectively combines the modal decomposition with the system's differential operators, capturing the influence of the resonators. It provides a complete representation of the beam's dynamic response, encompassing both the modal characteristics and the interactive effects of the resonators. \[ \begin{split}\mathcal{L}\sum_{m=1}^{N_{m}}\phi_{m}(x)z_{m}(t)+ \mathcal{C}\frac{\partial}{\partial t}\sum_{m=1}^{N_{m}}\phi_{m}(x)z_{m}(t)+ \mathcal{M}\frac{\partial^{2}}{\partial t^{2}}\sum_{m=1}^{N_{m}}\phi_{m}(x)z_{ m}(t)-\\ \sum_{r=1}^{N_{r}}\biggl{(}k_{r}z_{r}(t)+c_{r}\frac{dz_{r}(t)}{ dt}\biggr{)}\delta(x-x_{r})=\mathcal{F}_{b_{m}}(x,t)\end{split}  \tag{13}\]  
Multiplying Equation (13) by \(\phi_{n}(x)\) and integrating over the domain \(D\), and applying the orthogonality conditions Equations (8)-(11) of the mode shapes, gives  
\[\bar{z}_{m}(t)+2\zeta_{m}\omega_{m}\dot{z}_{m}(t)+\omega_{m}^{2}z_{m}(t)-\sum _{r=1}^{N_{r}}m_{r}\omega_{r}^{2}z_{r}(t)\phi_{m}(x_{r})=\mathcal{Q}_{b_{m}}( x,t),\quad m=1,2,\ldots,N_{m} \tag{14}\]  
Similarly for resonators, substituting the modal expansion Equations (12) into (2) yields  
\[\bar{z}_{r}(t)+2\zeta_{r}\omega_{r}\dot{z}_{r}(t)+\omega_{r}^{2}z_{r}(t)+\sum _{m=1}^{N_{m}}\bar{z}_{m}(t)\phi_{m}(x_{r})=\mathcal{Q}_{b_{r}}(t),\quad r=1, 2,\ldots,N_{r} \tag{15}\]  
Here, \(N_{m}\) and \(N_{r}\) denote the number of modes and resonators, respectively. Each mode has a specific modal frequency \(\omega_{m}\), and each resonator has a mass \(m_{r}\) and its own natural frequency \(\omega_{r}\). The damping ratios \(\zeta_{m}\) for the modes and \(\zeta_{r}\) for the resonators quantify energy dissipation. To simplify, the superscript \"dot\" indicates time derivatives, and \"prime\" indicates spatial derivatives. Each equation in the model represents the dynamics of modal coordinates or resonator displacement as a second-order ordinary differential equation. The dynamics are influenced by modal and resonator parameters (natural frequencies \(\omega_{m}\) and \(\omega_{r}\), damping ratios \(\zeta_{m}\) and \(\zeta_{r}\)), their interactions, and base excitation forces (\(\mathcal{Q}_{b_{m}}\) and \(\mathcal{Q}_{b_{r}}\)). Decoupling of these equations is achieved through an orthogonal transformation, involving pre- and post-multiplication by the mode shape matrix. This leads to diagonalization of the mass and stiffness matrices, thanks to the orthogonality of eigenvectors to both matrices. The result is a set of decoupled ordinary differential equations. This normal mode method applies in the absence of damping or with proportional damping, where the damping matrix is a linear combination of the mass and stiffness matrices. The transformation becomes orthonormal when the mode shape is normalized to the mass matrix. Combining the structural dynamics represented in Equation (14) with the dynamics of resonators from Equation (15) enables the coupling of inertial terms and decoupling of stiffness in the system, facilitating analysis in the frequency domain.