Multiband Pure Topological States in Elastic Structures

Inspired by notions of topological physics, recent years have witnessed the rapid development of mechanical metamaterials with novel properties of topological states. However, most of the current investigations have either focused on discrete mass-spring lattices, with topological states limited to a single operating band, or on various elaborate continuous elastic systems, enduring the drawbacks of modal couplings. It remains largely unexplored how to design topological elastic systems that naturally possess multiple operating bands and are free from modal couplings. In this study, we design an elastic system based on fundamental mechanical elements (beams, rods and nuts), which is capable of supporting multiband pure topological states. Through an equivalent beam-spring model with lumped masses together with finite element analysis, we demonstrate that our proposed structure exhibits multiple Dirac points (DPs) at different frequencies. We show that simply adjusting the heights of nuts fastened on beams can lift the degeneracies, giving rise to two kinds of valley Hall phases characterized by opposite valley Chern numbers. The dispersion diagram of the supercell formed by unit cells with different topological indices shows that there simultaneously exist perfectly pure interface modes (i.e., no other modes coexist) within two frequency ranges. Furthermore, numerical simulations demonstrate that the domain wall formed by structures with distinct topological properties supports topologically protected interface waves over dual frequency ranges. Our results have potential for the design of mechanical systems that need to work under changeable working frequencies and may have significant impact on many diverse fields such as vibration control, energy harvesting and seismic isolation.


INTRODUCTION
In the last decade, the discovery and extensive study of topological insulators in condensed matters [1] have brought innovative ideas to the design of classical wave systems with novel properties [2]. In electronic systems, non-trivial topological phases, manifested by robust and backscattering-free boundary states [3], are generally characterized by topological invariants, e.g., Chern numbers, which are defined by the eigenstates on energy bands [4]. The concept of band structure has also been introduced and largely studied in various artificial periodical structures such as photonic crystals [5] and phononic crystals [6]. Consequently, various concepts in topological physics are naturally introduced into classical wave systems [2], which provides many unprecedented new avenues for wave manipulation [7], energy harvesting [8] and information transfer [9]. In return, taking advantage of scalability and controllability, classical wave systems, e.g., acoustic crystals and mechanical structures, provide powerful and versatile platforms for probing into topological physics and related unusual phenomena [10].
Robust topological states exist at interface between phases with distinct topological invariants [11]. Despite the origin of the topological states varies dramatically in the existing proposed structures, the common requirement in their realization is breaking a certain kind of symmetry, namely, temporal or spatial symmetry [12]. For example, analogous to quantum Hall effect [13], active components [14][15][16], the application of an external field [17][18][19] and synthetic gauge fields [20,21], have been widely exploited to break time-reversal symmetry, opening topological bandgaps that support one-way edge modes. These approaches make the system complicated and obviously difficult to implement in practice. Recently, spatial symmetry have been proven to be one of the most accessible way to achieve non-trivial topological phases, since it relies on simple manipulating the symmetry of the underlying lattice. For instance, quantum spin Hall effect has been emulated by solely changing the geometric parameters of the structure [22,23], which breaks the point group symmetry and results in topologically non-trivial bandgap sustaining two dispersive edge states with opposite propagation directions [19,24,25]. Particularly, simple rotation or mirror operations have succeeded in realizing topological interface states based on valley Hall effect [26][27][28][29][30].
Here, we propose a way to trigger the valley Hall phase transition by adjusting the nuts fastened on Euler-Bernoulli (EB) beams, which is amenable to practical implementation.
In spite of the fact that great achievements have been made in designing mechanical systems possessing topological states, the majority of those investigations have focused on discrete massspring systems [14,19,[31][32][33]. Thanks to its structural simplicity, these discrete systems are excellent candidates for illustrating abundant topological essentials associated with various exotic phenomena [34,35]. In particular, the resulting topological interface modes hosted by mass-spring lattices are ideally pure, meaning that there are no additional trivial bands coexisting within the same frequency range [19,31]. This monomodal feature is beneficial to witness highly identifiable interface guided wave propagation without any interference. Nevertheless, the difficulty of transfer to practical configurations and the nature of limited degree of freedom induced frequency band singleness [19,27,31] make the discrete systems hopeless of practical applications. In recent years, tremendous efforts have been devoted to achieving topological elastic wave transmission in various continuous elastic systems [7,[36][37][38][39][40]. Unlike discrete systems, continuous structures often possess a wide range of frequency bands due to their infinite degree of freedom [41]. In terms of potential application, this multimodal feature offers a tremendous degree of flexibility in controlling wave dynamic response, providing opportunities for designing structures possessing broadband or multiband topologically protected waves [42][43][44][45][46]. However, the nature of multimode inevitably leads to a hybrid of the bands inside the topological bandgaps, which can severely impedes the performance of interface states [38,40,47]. In other words, ubiquitous modal couplings greatly increases the complexity and difficulty of topological band engineering for elastic waves [45,46]. From these points of view, designing elastic structures with topological states featuring a combination of multiple bands (i.e., with multiple working frequency ranges) and purity (i.e., free from modal hybridization) remains elusive and is rarely studied in the present literature.
In this paper, we design an elastic system composed of fundamental mechanical elements (beams, rods and nuts), which is capable of supporting pure topological states over multi-frequency ranges. Through an equivalent beam-spring model and finite element analysis, we demonstrate that the honeycomb structure composed of EB beams exhibits multiple Dirac degeneracy points at different frequencies, and adjusting the heights of nuts fastened on beams can lift the degeneracies, giving rise to two kinds of valley Hall phases characterized by opposite valley Chern numbers. Dispersion analysis of the supercell constructed by unit cells with distinct topological indices together with numerical simulations for a finite-size structure demonstrate that the domain wall formed by structures belonging to distinct topological phases supports interface waves at dual frequency ranges, which are topologically protected manifested by the feature of being immune to backscattering of sharp corner. The present study can improve the shortcomings of the most existing phononic topological insulators with only a single operating band, and is expected to open up the possibility for the realization of frequency switchable topological states.

METHODS
This section presents the free vibration analysis for the proposed elastic honeycomb structure based on theoretical model and finite element method (FEM).

Beam-Spring Model With Lumped Masses
Considering the preconditions for the emergence of DPs [48], we construct a honeycomb lattice structure using congruent circular columns coupled with each other via slender rods, as shown in Figure 1A, where, each column is equipped with two nuts whose spacing is adjustable. Note that all the nuts are identical, and colors red and green are only used to distinguish nonequivalent sites within one unit cell. In order to facilitate theoretical analysis, we propose an equivalent beam-spring model to simplify the honeycomb lattice structure, where, the circular columns are considered as EB-beams, the slender rods are regarded as the massless springs and the nuts are replaced by lumped masses, as shown in Figure 1B. A theoretical framework based on the beamspring model for the transverse free vibration analysis will be elaborated, which will be applied to calculation of the band structure of the infinite honeycomb lattice.

Governing Equations
Considering an arbitrary unit cell located at the site n 1 , n 2 { }of the infinite periodic lattice, with n 1 , n 2 { }∈ Z 2 , the governing equations for the motion of transverse free vibration of the two EB beams within the unit cell are as follows: where w n1,n2 i and v n1,n2 i denote the transverse displacement of the beam labeled i (i = 1, 2) along x and y-axes, respectively. E, I, ρ, A and L are Young's modulus of elasticity, the cross-sectional moment of inertia, density, the cross-sectional area and length of the beam, respectively. m is the equivalent lumped mass, attached at heights of h 1 and L − h 1 for the beam 1, and at h 2 and L − h 2 for the beam 2 (see Figure 1B). h c as well as L − h c is the position at which the horizontal springs are connected with the beams. Moreover, p n 1 ,n 2 i and q n 1 ,n 2 i , which represent the forces acting on the beam i (i = 1, 2) within the unit cell n 1 , n 2 { }along x and y-axes, respectively, due to the coupling to adjacent beams via connecting springs, are expressed as: where T are unit vectors characterizing the directions of connecting springs of the same stiffness k s , and u n 1 ,n 2 i denotes the displacement vector formed by the two components of the transverse displacement of the beam element i (i = 1, 2) within the unit cell n 1 , n 2 { }, that is, Using the time-harmonic assumption, i.e., w i = W i exp (iωt), v i = V i exp (iωt), Eq. 1 can be rewritten as: where ω is the angular frequency, the capital symbols (W and V) represent the corresponding time-harmonic solutions of Eq. 1, P n 1 ,n 2 i and Q n 1 ,n 2 i ,accordingly, are the time-harmonic form of Eq. 2, which are not repeated here.

Bloch's Theorem
For free vibration analyses in an infinite periodic lattice, Bloch's theorem states that the time-harmonic solution of the equation of motion can be written as a periodic function with a spatial modulation [49], that is, where k is the Bloch wave vector and U(x) is periodic with respect to the basis vectors of the lattice, hence satisfies In this case, Eq. 5 leads to the well-known Floquet-Bloch condition: hereinafter, notations p = exp (ik ·a 1 ) and s = exp (ik ·a 2 ) are adopted for the sake of brevity. Note that Bloch's theorem allows the vibration response of an infinite and periodic structure to be completely described through the governing equations for a unit cell, with suitable boundary conditions. The properties at any other location can be related to those of this very unit cell according to the formula expressed by Eq. 6.
By employing the Bloch's theorem, our governing equation Eq. 4 will turn into: with have been used, the superscript (p) denotes complex conjugation, and T i t i •t i (notation • denotes the outer product of two vectors) with i = 1, 2, 3. Equations 7,8 can be succinctly written as: with where 1 n×n is a n-by-n identity matrix, and notation ⊗ denotes the Kronecker product of two matrices. It is noteworthy that Eq. 9 is a fourth order differential equation with respect to z, while the unknown U itself is composed of the transverse displacements of the beam in the Oxy plane.

Eigenvalue Problem
In order to obtain the eigenvalue problem associated with the free vibration problem presented above, we use Galerkin method and the superposition of modes [50] to deal with Eq. 9, in which the time-harmonic Bloch solutions can be expressed as: Frontiers in Physics | www.frontiersin.org June 2022 | Volume 10 | Article 909820 where trial functionsW (n) (n 1, 2, . . . , N) are taken to be the normalized modes shapes of an individual free-free EB beam, expressed as Eqs. S.10 and S.16 in Supplementary Material. ξ 1(n) , η 1(n) , ξ 2(n) , η 2(n) are the corresponding superposition coefficients, constituting the eigenvector to be solved in the following. N is the number of truncation order of the trial function been adopted. Substituting Eq. 11 into Eq. 9, followed by taking inner product with every term of the set of the trial functions, and using the orthogonality conditions [detailed as Eq. S.18 in Supplementary Material], one can obtain: where, ψ denotes the eigenvector, formed by the assembly of unknown coefficients in Eq. 11, i.e., , and, accordingly, in which, with ω 0(n) (n = 1, 2, . . . , N) being the nth order natural angular frequency of the free-free EB beam, corresponding to the normalized natural modesW (n) . U h is a vector formed by the value of each order normalized eigenmode at z = h, that is, and the matrix T , defined in Eq. 10, is given by: Eq. 12 is equivalent to the well-known generalized eigenvalue problem: with K = Ω − S, M 1 4N×4N + mH and Λ = ω 2 . Note that the Hermitian property of the matrix T (and so that K) ensures that Eq. 14 has real eigenvalues. Consequently, solutions of Eq. 14 unequivocally yield the eigenfrequencies and eigenmodes of our infinite honeycomb lattice based on the theoretical model, which will be detailed in the next section.

Irreducible Brillouin Zone
In general, to obtain the band structures, we have to sweep every wave vector k in the first Brillouin zone (BZ) of the reciprocal space, spanned by the basis vectors expressed as: where a 1 [3a/2, − 3 √ a/2, 0] and a 2 [3a/2, 3 √ a/2, 0] are two basis vectors defining the periodicity of the honeycomb lattice, a 3 = [0, 0, 1] is assumed to be the unit vector e z . Thanks to the symmetry of our hexagonal lattice, the band structure can be calculated by considering only the wave vectors along the boundaries of the first irreducible Brillouin zone (IBZ) (see the inset of Figure 2), whose edges can be expressed as: Frontiers in Physics | www.frontiersin.org June 2022 | Volume 10 | Article 909820 Following the ideas of the theoretical method proposed in Sec. 2.1, solutions of displacement fields of the two beams within one unit cell are all expanded by the trial function truncated to order N [i.e., Eq. 11]. Then for each point along Γ-K-M-Γ, the eigenvalue problem Eq. 14 need to be solved to obtain the corresponding natural angular frequency ω and eigenvector ψ.

Finite Element Analysis
In this study, we perform finite-element simulations for the elastic waves in the honeycomb lattice structures using COMSOL Multiphysics with the solid mechanics module. The solid honeycomb lattice model is established, shown as in Figure 1A, where the material of two circular columns and nuts is selected as aluminum (E al = 70 GPa, ρ al = 2,700 kg/m 3 and ] al = 0.33) and that of slender connecting rods is nylon (E ny = 2 GPa, ρ ny = 1,150 kg/m 3 and ] ny = 0.4). Throughout this paper, the diameter of the slender rods is set to be d = 0.4 cm, and that for the circular columns is D = 2 cm. The nuts denoted by oblate columns have a diameter of 2D and a height of h n = 1.5 cm. All other parameters needed in modeling can be found in Table 1. We calculate the band structure of the infinite and periodic honeycomb structure using a single unit cell, wherein the Floquet-Bloch periodic boundary condition is applied in both the a 1 and a 2 directions. The eigenfrequency study is performed for every wave vector k scanned along the boundaries of the IBZ whose edges are expressed as Eq. 16.

Multiple Dirac Points
We first consider a periodic honeycomb structure in which the nuts on the left and right beams are fastened at the same height, i.e., h 1 = h 2 = h 0 , with h 0 /L = 0.25, where L is the length of the beam. Other unknown parameters are listed in Table 1. Note that the cross-sectional moment of inertia of the beam I = πD 4 /64, the equivalent stiffness of the connecting spring k s = E ny πd 2 /[4 (a − D)] and the equivalent lumped mass m = ρ al π[(2D) 2 − D 2 ]h n /4 are all obtained corresponding to the parameters adopted in the finite element model. Figure 2 demonstrates the band structure of the initial unit cell with h 1 = h 2 = h 0 , obtained by theoretical method based on the equivalent beam-spring model and FEM. Despite some slight differences, which was caused by neglecting the bending action of slender rods in the equivalent beam-spring model, the theoretical results and FEM both show two Dirac degeneracy points at frequencies around 400 and 1,100 Hz at K point. It is worth noting that the emergency of DPs is ensured by C 3v symmetry of the underlying lattice [48], while the multiplicity of DPs is attributed to the multimode characteristic of transverse free vibration of beams (referring to the mode shapes displayed in Figure 5). Since that the emergency of DP is the prerequisite for engineering valley Hall phases [51], one can expect to achieve topological states within the dual working frequency ranges by simultaneously gaping the two DPs aforementioned. This will be detailed below.

Mirror Symmetry Breaking and Topological Phase Transition
It has been widely recognized that the Dirac degeneracy can be lifted by breaking spatial symmetries [48], which provide simple ways to achieve topological order. In this study, we break the mirror symmetry of the primitive lattice by adjusting the fastening heights of nuts such that h 1 is no longer equal to h 2 . We define the strength of perturbation as: where h 0 is the initial height of the nuts fastened on both the left and right beams. Note that h 1 and h 2 deviate from the unperturbed value h 0 with the same perturbation strength but opposite sign. In this case, two different symmetry-invariant configurations of the perturbed unit cell, shown as Cell A (β > 0) and Cell B (β < 0) in Figure 3, are created by breaking the mirror symmetry in two opposite ways. Figure 4A demonstrates the band structures of the honeycomb lattice arranged with unit cells of type-A (with β = 0.076), obtained by both the theoretical method and FEM. We can notice that the splittings of the two DPs occur and two complete bandgaps (indicated by cyan and gray region) emerge at frequencies around 400 Hz and 1,050 Hz. The consistency of   bandgap widths shows again that the theoretical model is in good agreement with FEM.
Interestingly, we find that the two bandgap widths show dramatically different trends with the adjustment height of the nuts (i.e., with different values of β). The inset of Figure 4A displays the evolution of bandgap widths as functions of β, wherein Δf quantifies the frequency difference between K + 1 (K + 2 ) and K − 1 (K − 2 ). It clearly indicates that the width of the first bandgap changes linearly with the change of β, while that of the second bandgap seemingly changes quadratically. Nevertheless, there exists an intersection between the two curves (with β of about 0.076 and Δf of about 41.5 Hz), which corresponds to the optimal point at which the two bandgaps have nearly the same size. The band structures presented in Figure 4A are just obtained upon this optimal point. We emphasize that this provides an accessible way to tune the bandgaps (namely, the operating frequency ranges of the topological states shown in the following text). Specifically, if we adjust the fastening heights of the nuts such that the perturbation strength β is larger than 0.076, then the lower frequency bandgap is dominant; while for β smaller than 0.076, the higher frequency bandgap dominates.
By analogy with the valley in quantum systems [26], two pairs of states of two band extreme at K point called as valley states appear near the two lifted DPs respectively, denoted by K + 1 (K + 2 ) and K − 1 (K − 2 ) marked in Figure 4A. In order to investigate the valley states in the honeycomb structures and to gain the full knowledge of the topological properties related to them, we demonstrate the evolution of bounding frequencies of the two gaps with β varying from −0.12 to 0.12 along with the modal shapes associated with the valley states for β = ±0.076, as shown in Figure 5. Both for the lower and higher bandgap, it is clearly shown that the two bands will cross each other when β transits from positive to negative, with degeneracies occurring when β = 0, which corresponds to the initial unperturbed unit cell. When β changes its sign (i.e., from unit cell of type A to B or vice versa), both the two bandgaps experience an opening-closing-reopening process, accompanied by the well-known band inversion phenomena [28], manifested as up and down exchange of eigenmodes corresponding to the valley states K + 1 (K + 2 ) and K − 1 (K − 2 ). Specifically, one can notice that, for β > 0 (Cell A), the eigenmodes for the upper (lower) band feature left-hand (right-hand) circularly polarized mode of the left (right) circular column entity; while for β > 0 (Cell B), the eigenmodes for the upper (lower) band feature right-hand (left-hand) circularly polarized mode of the right (left) circular column entity. This symmetry-breaking induced band-inversion phenomena evidently suggests that the topological phase transition occurs when β varies from positive to negative.
In addition, topological phase transition with the sign of β flipping can be quantitatively inferred by the valley Chern number, which is an index to characterize the topological nature of the band [27,33] and obtained by integrating the Berry curvature over an individual valley, expressed as: where, v = K(K′) is the BZ corner, the integration domain BZ(v) is a local area near the point v, and F denotes the Berry curvature associated with the bands forming the valley v. It has been shown that by means of k ·p perturbation method [19], the nonvanishing valley Chern numbers are related to the sign of geometrical perturbation β, that is, where sgn(β) means the sign of β, herein depends on the relative fastening height of the nuts on the left (h 1 ) and right (h 2 ) beam. Due to the fact that the two types of unit cell possess the opposite sign of β (see Figure 3), one can immediately realize that the valley Chern numbers of them also have opposite signs, representing two distinct valley Hall phases. For the sake of verification, we calculate the Berry curvature around valleys by the numerical method [14], using the eigenvectors stemming from the theoretical model and Bloch eigenmodes from FEM, respectively. For β = 0.076, the combined Berry curvatures associated with both the two bands above the two bandgaps are shown in Figure 4B (calculation details are presented in Supplementary Material). For both the results from our theoretical model and FEM, the integral of the Berry curvature over the full BZ is precisely zero with the total Chern number C = 0, which is ensured by time-reversal symmetry [3]. For both lower and higher frequency bands, the associated Berry curvatures are prominently localized near valleys K and K′, and feature exactly the same amplitude but opposite signs with respect to K and K′. By integrating these Berry curvatures (taking the results from FEM as an example), we obtain the valley Chern numbers C K(K′) = +(−)0.37 (lower bands) and C K(K′) = +(−)0.29 (higher bands).
Recalling that, endowed with ± β, the two types of unit cell are symmetry-invariant, implying that K(K′) for type-A unit cell is just equivalent to K′(K) for type-B. Therefore, if C K = +0.37 (+0.29) characterizes the topological properties associated to the first (second) bandgap for type-A unit cell, then C K′ = −0.37 (−0.29) corresponds to type-B unit cell. In spite of apparent discrepancy between numerical valley Chern numbers and the theoretical predictions (±1/2), the different signs, + and −, clearly distinguish different valley Hall phases. It is worth emphasizing that, although we only care about type-A unit cell above, the type-B unit cell with opposite value of perturbation strength (β = −0.076) has exactly the same dispersion diagram as that of A, though band inversion is implied with the eigenmodes associated with the corresponding valley states flipped [27], which echoes the opposite valley Chern numbers aforementioned. Therefore, topological phase transition occurs when we cross the domain wall between lattices formed by type-A and type-B unit cells, and according to the bulk-edge correspondence [1], topological states are bound to appear in the two bandgaps shared in common by these two lattices.

Interface Modes
To verify the existence of elastic valley Hall interface states, we consider a strip supercell composed of 14 unit cells with different topological valley Hall phases at either side of the interface, as shown in Figure 3. Note that the supercell is arranged periodically along the a 1 direction, with periodicity 3 √ a, and has finite size along the a 2 direction. Therefore, the wave vector k p along which the projected band structure of the supercell shall be calculated can be restricted to −π/( 3 √ a) ≤ k p ≤ π/( 3 √ a), in the direction of a 1 [3a/2, − 3 √ a/2, 0]. During the numerical calculation process, the Floquet-Bloch boundary conditions are in the a 1 direction, while the upper-right and lower-left edges are set free.  Figure 4A) as a function of β. The insets showing the displacement fields and modal polarization associated with eigenmodes of K − 1 and K + 1 states for β =±0.076 reveal the band inversion and thus manifest the topological phase transition as β varies from positive to negative. Note that the bandgap is filled with different colors to identify different topological phases (B) The same as (A), but the phase transition at frequency around the second DP is revealed.  Figure 3, with β =±0.076). Colour indicates the degree of interface confinement of the eigenmodes, defined by the ratio α in Eq. 20. The enlarged band diagram clearly shows perfectly pure interface modes (i.e., no other modes coexist) spanning the two bandgaps, highlighted by red dots (B) Displacement fields of the interface modes at two frequency ranges indicated by black arrows, together with those of the bulk modes.
The corresponding band structure of the strip supercell is presented in Figure 6A, wherein the color bar indicates the degree of localization of the mode at the interface by making a ratio α between the average displacement nearby the interface and in the whole supercell, defined as: where I 2 represents two unit cells on each side of the interface, SC means the whole supercell, and u x , u y , u z are three displacement components in x, y, and z axis, respectively. Notably, α would be close to 1 (red) for modes localized near the interface and much less than 1 or near 0 (blue) for bulk modes. From the enlarged view of the band structure, we can find that the interface modes (denoted by the red dots) appear simultaneously in the overlapping bulk bandgaps discussed above, with working frequencies around 400 Hz and 1,050 Hz. It is worth noting that the interface modes spanning the two bandgaps are perfectly pure, namely, no other modes coexist at the same frequency ranges. This extraordinary feature, as one of the essential points of this paper, has rarely been reported and can hardly been realized in the existing topological elastic structures wherein modal hybridization is ubiquitous and unavoidable [2].
In order to demonstrate the localized properties of the interface modes, we plot the displacement fields for some typical eigenmodes evaluated at four frequencies indicated by the black arrows, as shown in Figure 6B. It can be seen clearly that, for the interface mode, either in low-frequency or highfrequency range, only beams close to the interface have a strong vibration, while all the beams of the supercell have considerable displacement amplitude for the bulk mode. What's more, it is worth observing that both the low-frequency and high-frequency interface modes are gapless, signifying that two bandgaps resulted from the splittings of DPs (see Sec. 3.2) are just the expected operating frequency ranges of the topological elastic waves.

Topologically Protected Interface States at Multi-Working Frequencies
To gain further insight into the topological interface states observed in the previous section, we construct a typical configuration with zigzag waveguide (shown as the top-left panel of Figure 7) to carry the numerical simulation of topological elastic wave propagation property. Simulations are carried by the frequency domain study of solid mechanics in the COMSOL. During the simulation, tips of the two beams located near the black star are subjected to a horizontal harmonic displacement excitation with some specified frequency, while other beams and nuts, including those on the boundaries, are all set free. The displacement fields, indicating the dynamic response, are drawn to reveal the characteristics of the wave propagation at this very frequency.
The contours of displacement fields of dynamic responses correspond to the excitation frequencies 415 Hz and 1,088 Hz are displayed in the bottom-left and bottom-right panels of Figure 7, respectively. As we can see, for both lower and higher frequencies, elastic waves are localized and can pass the zigzag interface without any hindrance. This is exactly what we expected since the two excitation frequencies are both located in the bandgaps in which the valley Hall interface modes occupy, as we demonstrated in Sec. 3.3. For comparison, we map the dynamic response displacement field for excitation frequency 350 Hz, corresponding to the bulk band, as displayed in the topright panel of Figure 7. Unlike the former cases in which displacement fields are concentrated near the interface, for the bulk band frequency, the structure experiences almost a full-field dynamic response. These results are in complete agreement with the previous dispersion analysis, and provide strong evidence that our proposed structure supports topological elastic wave propagation with multi-working frequencies.

CONCLUSION
This study present a strategy to design elastic structures exhibiting pure topological states over multi-frequency ranges. We proposed an equivalent beam-spring model with lumped masses for investigating dispersion relations of the designed structures, together with numerical simulations using FEM for verification. The band structures for the initial unit cell obtained by these two ways, agreeing well with each other, document two Dirac points at different frequencies. The emergence of these deterministic Dirac points is ensured by the C 3v symmetry of the underlying lattice and the multiplicity is attributed to the nature of multi-order natural modes of the EB beams consisting our structures. We demonstrated that simply adjusting the heights of nuts fastened on beams can lift the degeneracies, resulting in distinct valley Hall phases characterized by opposite valley Chern numbers. The dispersion analysis of the supercell formed by unit cells featuring different topological phases indicated that there simultaneously exist perfectly pure interface modes (i.e., no other modes coexist) within two frequency ranges. These multiband pure topological states are strongly confirmed by numerical simulations. It is worth noting that the purity and the multiplicity of topological states hosted by our elastic structure just right correspond to the merit of discrete lattices and that of continuous systems, respectively. In addition, we pointed out that adjusting the heights of fastening nuts can not only trigger the topological phase transition but also provide an accessible way to tune the operating bandgaps for topological states, although we only focused on the scenario that the two frequency ranges of interest are roughly equal in this paper.
Our proposed strategy open new possibilities of managing multiple elastic wave modes and is undoubtedly promising for multiband applications, e.g., multiband filters and multiband waveguide. Besides, composed of ubiquitous mechanical elements, our proposed structure would be amenable to physical implementation, hence we envisage that this architecture can be exploited in practice as an excellent and versatile platform to explore other intriguing wave propagation phenomena brought about by richer topological notions.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.