Breakdown of Conventional Winding Number Calculation in One-Dimensional Lattices with Interactions Beyond Nearest Neighbors

Topological insulators hold promises to realize exotic quantum phenomena in electronic, photonic, and phononic systems. Conventionally, topological indices, such as winding numbers, have been used to predict the number of topologically protected domain-wall states (TPDWSs) in topological insulators, a signature of the topological phenomenon called bulk-edge correspondence. Here, we demonstrate theoretically and experimentally that the number of TPDWSs in a mechanical Su-Schrieffer-Heeger (SSH) model can be higher than the winding number depending on the strengths of beyond-nearest-neighbor interactions, revealing the breakdown of the winding number prediction. Alternatively, we resort to the Berry connection to accurately characterize the number and spatial features of TPDWSs in SSH systems, further confirmed by the Jackiw-Rebbi theory proving that the multiple TPDWSs correspond to the bulk Dirac cones. Our findings deepen the understanding of complex network dynamics and offer a generalized paradigm for precise TPDWS prediction in potential applications involving localized vibrations, such as drug delivery and quantum computing.

One illustrative example of employing a topological invariant to determine the nontrivial topologically protected domain-wall states (TPDWSs) can be seen in the one-dimensional (1D) Su-Schrieffer-Heeger (SSH) model [28,29] as shown in Fig. S1a, b in Supplementary Note 1. Initially introduced to study solitons in polyacetylene, the SSH model was later adapted in mechanical systems to identify TPDWSs via a winding number calculation [30][31][32].As discussed in the references above and the "Analysis of One-Dimensional Su-Schrieffer-Heeger Model" subsection in Methods, the two arrangements of isomers with different spring stiffness c 1 and c 2 in Fig. S1b represent two topologically distinct phases.When c 1 > c 2 , the origin is excluded in the contour plot in the complex plane of the off-diagonal term in the stiffness matrix, C(k), where k is the wave number in the reciprocal space, and thus, the winding number n = 0, signifying a trivial intra-cell-hopping phase.In contrast, when c 1 < c 2 , the contour winds about the origin once, i.e., n = 1, indicating a topologically nontrivial inter-cell-hopping phase, see Fig. S1c in Supplementary Note 1.These gaugedependent winding numbers can also be evaluated via the Zak phase [33] measuring the rotation of eigenvectors in the unit cell, see Fig. S2 in Supplementary Note 1.Alternatively, the winding number can be directly calculated from C(k): where σ 3 is the third Pauli matrix, and C ′ is similar to an effective Hamiltonian, which is a chiral matrix obtained from: where σ 0 is the identity matrix.Seaming two phases with different ns creates a domain wall where a localized mode emerges.Such a domain-wall state is topologically protected due to the intrinsic topological phase difference between the two domains, i.e., the bulk-edge correspondence.
The above discussion has been well understood and applied to systems of higher dimensions, such as the 2D quantum valley Hall effect in phononic crystals [12,18,32,[34][35][36][37].Most of these studies can be simplified using mass-spring systems considering only the nearest neighbor (NN) interactions.Recently, arising attention has been devoted to mechanical metamaterials with lattice interactions beyond nearest neighbors (BNNs), achieving roton-like acoustic dispersion relations under ambient conditions similar to those observed in correlated quantum systems at low temperatures [38][39][40][41][42][43].In addition, intriguing topological states also arise due to such BNN coupling, including the increased winding number corresponding to a higher number of edge states due to larger BNN differences, as reported in previous studies [44][45][46].
In this study, we report the impact of BNN couplings on bulk-edge correspondence in addition to the increased winding number.Despite being commonly believed that the number of TPDWSs is governed by the difference in winding numbers between the two domains, we prove that the number of TPDWSs in the SSH model is not dictated by the winding number, but by the Jackiw-Rebbi (JR) indices associated with the JR zero modes [47].In previously studied SSH models, the winding numbers and JR indices happened to predict the same number of TPDWSs.However, our investigation reveals that such a coincidence is not generic, i.e., in the presence of BNN couplings, these two indices can significantly deviate from each other.In such a generic setup, we prove analytically and verify numerically and experimentally, that the JR indices always correctly predict the number of TPDWSs, while the relationship between the winding number and TPDWSs generally fails.We note that the discrepancy between the two is not a special scenario in an SSH model but a rather generic phenomenon across all topological Maxwell lattices and chiral matters [48].We, thus, propose to use the Berry connection with distinguishable local winding numbers as an alternative topological index to identify TPDWs, which also applies to a broader range of lattices with BNN interactions beyond those presented in the main text.

Mass-Spring Model Analysis
To start with, we add third-nearest neighbors (TNNs) with spring stiffness c ′ to a 1D mass-spring chain of lattice spacing a with a NN spring stiffness c and restrict the motion of identical masses to the horizontal direction, as presented in Fig. 1a.The phonon dispersion of a pair of masses reveals that when c ′ < 1/3c, the acoustic and optical phonon bands cross at k = π/a, protected by the space inversion symmetry (SIS).When c ′ > 1/3c, two additional band crossings emerge in the irreducible Brillouin zone (IBZ), as presented in Fig. 1b, c.Derivation of the exact locations of the Dirac points due to the existence of the TNN is presented in Eqns.10-17 in the "Analysis of One-Dimensional Su-Schrieffer-Heeger Model with Third-Nearest Neighbors" subsection in Methods.The additional band folding due to strong c ′ s results in negative group velocities in the acoustic phonon branch, corresponding to the backward wave observed in the previous study [38].We then break the SIS by applying a small perturbation to the NN spring stiffness c, i.e., making c 1 = 0.8c and c 2 = 1.2c, while maintaining all TNNs identical, which opens a band gap between acoustic and optical bands, as shown in Fig. 1d.The winding number calculation from C ′ (k) in Eqn. 13 for the two isomers of such a system suggests that, regardless of the strength of c ′ , the difference between two phases is always one, indicating one TPDWS at the domain boundary of the two phases.Note that, with the existence of c ′ , the contour plots in the complex plane are no longer circular as those shown in Fig. S1c.With weak c ′ (for example, c ′ = 1/10c), they present oval shapes, as shown in Fig. 1e, while strong c ′ (such as c ′ = c) creates two additional loops along the path, Fig. 1f.In either case, the circuit winds around the origin exactly once when c 1 < c 2 , indicating the topological charge being n = 1, while excluding it when c 1 > c 2 , thus n = 0, yielding a consistent n difference.
To verify the number of TPDWSs predicted with winding numbers, we consider a supercell containing 301 masses with a domain-wall mass at the center connected to soft or stiff springs on both sides, as shown in Fig. 2a, b, about which are the symmetrically arranged two phases with different ns.Bloch conditions are applied at the two ends of the chain to mimic an infinitely long chain for phonon dispersion calculation.Details of the Bloch conditions for the supercell are presented in the "Supercell Analysis of the Su-Schrieffer-Heeger Model" subsection in Methods.As predicted, when connected by weak TNNs, i.e., c ′ = 1/10c, only one symmetric (asymmetric) TPDWS exists within the bulk bandgap when the domain-wall mass is connected by soft (stiff) NN springs, as shown in the supercell band structures in Fig. 2c with corresponding edge and bulk mode shapes presented in Fig. 2e-j.However, when c ′ = c, we identify two additional edge modes in the bulk bandgap, Fig. 2d, violating the aforementioned winding number prediction.Mode shapes of these emerging bands in the bulk bandgap shown in Fig. 2k-m, p-r, confirm the localization of displacements at the domain wall, distinguishable from the bulk modes, as shown in Fig. 2n,o,s,t.When the domain-wall mass is connected by two soft (stiff) springs c 1 = 0.8c (c 2 = 1.2c), we obtain one (two) symmetric and two (one) asymmetric displacement fields about the domain-wall mass, as can be seen in Fig. 2k-m (p-r).
Spatial Fourier transform (SFT) of the mode shapes presented in Fig. 3 show a significantly widened peak width in all these edge modes compared to their bulk counterparts, suggesting a faster spatial decay of the vibration from the domain wall, evident of an edge mode.It is worth noting that, due to strong TNN interactions, the additional band crossing at k = π/2 as shown in Fig. 1c when c 1 = c 2 (the location of which is expressed in Eqn. 17 in the "Analysis of One-Dimensional Su-Schrieffer-Heeger Model with Third-Nearest Neighbors" subsection in Methods) results in the global peak (valley) of the bulk acoustic (optical) mode occurring at k = π/2, instead of k = π where a local peak (valley) appears, as in Fig. 2d.Hence, sharper SFT peaks appear at k = π/2 in the bulk modes closest to the bandgap (i.e., SC1/2-B1/2), as shown in Fig. 3b.In-gap TPDWSs are presented as widened peaks located at k = π/2, π, and 3π/2, as shown in Figs.3b and S3 in Supplementary Note 2, indicating their rapid spatial decay away from the domain wall with hybridized wavelengths.

Jackiw Rebbi Zero Modes
These domain-wall states can also be characterized by a massless Dirac theory.For c 1 > c 2 (c 1 < c 2 ), the breaking of the SIS introduces a positive (negative) mass to each Dirac point.Due to the mass sign flipping at the domain wall, one TPDWS in the bandgap, known as the JR zero mode, is expected to arise at the domain wall for each Dirac cone [47].Thus, the number of crossings presented in Fig. 1c within the IBZ  equals that of TPDWSs when the bandgap is open.Such an agreement also strongly resembles those in the quantum valley Hall effect in 2D, where the number of in-gap TPDWSs also matches that of bulk Dirac cones [12].A comprehensive demonstration of the existence of TPDWSs due to the hybridization of JR modes corresponding to the three Dirac points within the IBZ and their analytical solutions characterizing the spatial decay are presented in the "Derivation of the Jackiw Rebbi Zero Modes" subsection in Methods.As presented in Figs. 3, and S3, SFT of the analytical JR zero modes match well with those obtained from the supercell analysis, indicating successful prediction of TPDW properties using the JR theory, as well as the breakdown of the winding number prediction.In principle, regardless of the strengths of c ′ , the analytical solutions of JR modes in the SSH model should all stay at ω/ω 0 = (c 1 + c 2 + 2c ′ )/m (or ω 2 /ω 2 0 = 0 if plotting the eigenvalues of C ′ (k) in Eqn.13).However, this is not the case if we simply flip the arrangement of springs about the domain-wall mass as presented in Fig. 2a, b.One reason is that the domain wall setup breaks the chiral symmetry of the stiffness matrix C(k) since the total spring constants about the domain-wall mass is 2c 1 + 2c ′ (2c 2 + 2c ′ ) for SC1 (SC2), different from c 1 + c 2 + 2c ′ around other masses in the supercell.Such a perturbation can be readily fixed by pinning the interface to the ground using an additional spring with a constant of c 2 − c 1 (c 1 − c 2 ) for SC1 (SC2), as shown in Fig. S4a (b) in Supplementary Note 2. The two TPDWSs due to different domain walls for weak TNNs (for example, c ′ = 1/10c) then become degenerate at ω/ω 0 = (c 1 + c 2 + 2c ′ )/m, or ω 2 /ω 2 0 = 0 if removing the diagonal elements of C(k), as shown in Fig. S4c.However, for the ones with strong TNNs (for example, c ′ = c) where multiple TPDWSs exist at one domain wall, such a fix can only bring one of the TPDWSs to ω/ω 0 = (c The other two TPDWSs are located symmetrically above and below it, as shown in Fig. S4d.The remaining shift in frequency is due to the hybridization of two JR modes with the same parity.Details about TPDWS symmetry and parity, as well as their hybridization conditions causing the frequency shift, are all discussed in "Jackiw Rebbi Mode Parity and Hybridization" subsection in Methods.The shifting of the zero-frequency(energy) domain-wall modes (after pinning the domain-wall mass) to finite frequencies is not unique to the SSH model, instead, it is a generic feature existing in most known 0D topological modes in various types of topological insulators (such as corner modes in 2D higher-order topological insulators), whose energies are also sensitive to local perturbations near the localized modes [49,50].

Berry Connection
The question then arises as to how to determine TPDWSs using a topological descriptor associating the spectral evolution of the eigenvector with these states.A closer examination of the contour plots in Fig. 1e, f and the winding number calculation in Eqn. 1 suggest that, although the difference in n is one regardless of the TNN strength, trajectories of the contour plots, or the integrand of Eqn. 1, i.e., the Berry connection, varies with c ′ , where C ′ is expanded to Eqn. 13 in the "Analysis of One-Dimensional Su-Schrieffer-Heeger Model with Third-Nearest Neighbors" in Methods to include the TNNs.Since the number of TPDWSs depends on the topological invariant difference due to different gauges, we plot ∆B(k , for unit cells with different c ′ in Fig. 4 to describe its topology.When c ′ = 0, only one peak exists at k = π/a in ∆B, corresponding to the Dirac point at k = π/a in the band structure in Fig. 1c.As c ′ increases while c ′ < 1/3c, this peak at k = π/a decreases and widens until it splits into two smaller peaks (such as when c ′ = 1/3c).As c ′ continues to increase, the valley at k = π/a dips below ∆B(k) = 0 while the two positive peaks drift apart with locations matching Dirac points as expressed in Eqn.17, until k = π/3a and 5π/3a, as discussed in the "Analysis of One-Dimensional Su-Schrieffer-Heeger Model with Third Nearest Neighbors" in Methods.Meanwhile, the two peaks and one valley are further sharpened as c ′ increases.The integral around each peak (valley) is ±1, i.e., yielding a local winding number.It is worth noting that the total integral over the IBZ does not change as c ′ varies, yielding a consistent winding number of n = 1.The transition from one peak in ∆B into two peaks and one valley agrees with the change of TPDWS counts with corresponding c ′ .Moreover, the locations of the peaks/valleys informing of the TPDWS wavelengths also agree with those calculated from JR zero modes demonstrated in "Derivation of the Jackiw Rebbio Zero Modes" and the supercell calculation presented in "Supercell Analysis of the Su-Schrieffer-Heeger Model" in Methods, the results of which are plotted in Figs. 3, S3, and S4u, v in Supplementary Note 2 if fixing the domain wall to make C(k) chiral.
To understand such a transition from one peak to two peaks and one valley with increased c ′ , one can draw an analogy between the evolution of ∆B as c ′ decreases and the inter-valley mixing of the Berry curvature in our previously studied valley Hall effect [12].In the current SSH model with TNNs, when c ′ ≫ c, the perturbation induced by the NNs, c, is relatively small, resulting in minimal inter-valley mixing between the two peaks and one valley in ∆B, which distinctively exist in the IBZ, matching the three TPDWSs within the bulk bandgap.As c ′ weakens, the difference in c becomes more prominent, introducing stronger SIS perturbation, and thus an enhanced peak-valley mixing closer to the valley at π/a, and eventually merging all into one single peak at k = π/a, leaving only one TPDWS in the bandgap.
It is worth noting that, compared to conventional winding number calculations, the Berry connection prediction alluded to is not limited to making correct TPDWS predictions in lattices with identical TNNs.As discussed in subsection "Topologically Protected Domain-Wall States beyond Equal Third-Nearest Neighbors" in Methods, one can also predict the number and the wavelengths of TPDWS when TNNs are nonidentical, i.e., c ′ 1 ̸ = c ′ 2 , as well as for systems with interactions beyond TNNs, whose results are shown in Figs.S5 and S6, respectively, in Supplementary Note 3. The additional TPDWs due to BNNs always arise in pairs with local integrals of Berry connection around the peaks and valleys being ±1, respectively.Information about wavelengths acquired from such a Berry connection analysis is also unattainable using conventional winding number calculations.Hence, the Berry connection provides a generalized methodology supplying enriched information about TPDWSs in lattices with complex networks.

Laser-Assisted Experimental Characterization
We proceed now to conduct experiments on 1D specimens adapted from an existing TNN model [38,39], as shown in Fig. 5a-d.Information regarding the experimental specimens is listed in subsection "Experimental Fabrication and Characterization" in Methods.As presented in Fig. 5d, each unit cell contains a pair of masses connected by alternating stiff and soft NN bars and identical TNN frames.The domain-wall mass behind the frame labeled as E1 in Fig. 5 f is connected by two stiff struts (blue bars), about which are placed with 8 unit cells with opposite stiff and soft NN arrangements.
The lattice specimens are hung by a string from the top and are excited in the ydirection using an electrodynamic shaker (PCB 2007E01, powered by a Krohn-Hite 7500 amplifier) placed off-centered near the bottom left end for torsional excitation.Velocities of the left and right ends of each frame in Fig. 5b, c in the y−direction are measured by a scanning laser Doppler vibrometer (SLDV, Polytec PSV-500), and their differences are recorded as the torsional velocities about the z-axis.Note that shear deformations in the y-direction will also be recorded simultaneously.However, these shear modes do not exist in our frequency range of interest.For comparison, unit cell and supercell analyses containing a domain wall with the same dimensions and material parameters are also performed using the finite element method with COMSOL Multiphysics with results presented in Fig. S7 in Supplementary Note 4.
We then prescribe a chirp excitation sweeping from 800 to 1400 Hz to the lattice with strong TNNs, Fig. 5b (500 to 700 Hz to the one with weak TNNs, Fig. 5c) and measure torsional velocities.For the specimen with strong TNNs and without a domain wall, we achieve an excellent agreement between the spatiotemporal spectral response obtained from a discrete Fourier transform of the torsional velocity sampled along the axial direction of the specimen (i.e., the z-axis in Fig. 5) and the acoustic and optical torsional phonon branches predicted from the unit cell analysis, as shown in Fig. 5e, presenting a roton-like dispersion relation.Frequency responses of the experimentally measured torsional velocities of the frames near the domain wall and in the bulk of the specimen with strong (weak) TNNs reveal three (one) distinct peaks within the bulk bandgap with amplified torsional velocities in proximity to the domain wall, as evident in Fig. 5g (h).These in-gap peaks agree with the TPDWSs predicted in supercell band structures containing matching lattice configurations with Bloch boundary conditions applied at two ends, presented as dark solid lines in Fig. 5g, h.Snapshots of the torsional velocity fields experimentally measured at bulk and TPDWS frequencies are shown in Fig. 5i and j, corresponding to strong and weak TNNs, respectively, agreeing with mode shapes calculated using the finite element method plotted in Fig. S7 in Supplementary Note 4. Symmetries of these measured and calculated TPDWSs concur with those in the toy model presented in Fig. 2p-r and f, i.e., the three torsional TPDWSs in the lattice with strong TNNs investigated the experimentally measured Fig. 5i, simulated Fig. S7a, and the mass-spring model in Fig. 2p-r, are asymmetric, symmetric, and symmetric about the domain wall (corresponding to the red, yellow, and green peaks in Fig. 5g), and the one TPDWS in the lattice with weak TNNs shown in Figs.5j, S7b, and 2f are all asymmetric (corresponding to the red peak in Fig. 5h), confirming that these domain-wall modes are indeed TPDWSs predicted in theory.

Conclusions
We theoretically and experimentally reveal the breakdown of the conventional winding number prediction of TPDWSs in SSH lattices with BNN interactions.We, instead, propose to count the local winding numbers by calculating the Berry connection characterizing the evolution of eigenvectors in the reciprocal space to obtain the correct number of TPDWSs.Moreover, Berry connection offers more insights into the TPDWSs, including their wavelengths and spatial decay rates.Further, we demonstrate that these TPDWSs are the phonon realization of JR zero modes, analytically validating the Berry connection prediction.Note that the discordance between the total winding number in IBZ and the counting of topological modes isn't an exception; rather, it is a common characteristic universally seen across topological states governed by winding numbers.For example, a similar discrepancy is also evident in topological Maxwell lattices and chiral matters [48], where the accurate enumeration of topological modes requires the winding number to be adjusted by the addition of integer numbers.Such amendments encapsulate the intricate aspects of the system's physics, such as lattice structures, gauge choice, and local counting.Furthermore, our study provides a more generalized paradigm in accurate topological state predictions in lattices beyond 1D with TNNs, and is applicable to a broader range of complex systems with multinodal interactions [21,51], especially at the nano- [52] and microscales [53], where BNN interactions commonly exist.Successfully identifying and achieving mechanical/vibration topological states in complex systems can also inspire solutions in other realms of science where predictable and precise vibration modes are critical, such as facilitating drug delivery [54,55] and advancing quantum information processing using phonons [56][57][58] in quantum technologies, such as quantum computing.

Analysis of One-Dimensional Su-Schrieffer-Heeger Model
The governing equations of a 1D SSH lattice unit cell shown in Fig. S1 in Supplementary Note 1 can be expressed as: where displacements of the two masses in the n-th cell are denoted as u n 1 and u n 2 , respectively, and can be expressed using a plane-wave solution in combination with Bloch-Floquet periodic boundary conditions: where ω is the vibration frequency, u n are the displacements of the n-th cell with , k is the wave number, which is inversely proportional to the wavelength λ, i.e., k = 2π/λ, a denotes the lattice constant, ũ(k) are displacements within the unit cell.Substituting this expression in Eqns. 4 and 5 gives: where C(k) is the stiffness matrix of the periodic system: Assume a = 1 and divide the stiffness matrix C(k) by c 2 , we can then plot the offdiagonal element of C(k), i.e., ρ(k) = c 1 /c 2 + e ika , and project ρ(k) to a complex plane, as shown in Fig. S1c in Supplementary Note 1.
The winding number difference between the two gauges can also be characterized by the Zak phase: where ũ− (k) = [ρ * (k)/|ρ(k)|, 1]/ √ 2 is the eigenvector corresponding to the smaller eigenvalue of the matrix C(k).Writing ρ(k) = |ρ(k)|e iϕ(k) , we find . This implies that the Zak phase measures the change in the phase of the first component (in this particular choice of gauge) of the eigenvector as wavenumber k changes from −π/a to π/a.When c 1 > c 2 , Z = 0, suggesting no changes in the phase difference of the eigenvectors across IBZ, see the blue curve in Figs.S1c and S2a.When c 1 < c 2 , Z = 1, indicating a phase change of 2π within the IBZ, see the red curve in Figs.S1c and S2b in Supplementary Note 1.

Analysis of One-Dimensional Su-Schrieffer-Heeger Model with Third-Nearest Neighbors
Adding TNNs, as in Fig. 1a, modifies the governing equations to: Plugging the same wave ansatz in Eqn.(6) in the above set of equations, we obtain an equation of the same form as Eqn.( 7) with the following stiffness matrix: The chiral matrix then becomes: The contour plots of the off-diagonal element of this matrix ρ(k) = c 1 + c 2 e ika + c ′ 1 e −ika + c ′ 2 e 2ika for a complete circuit of k from k = 0 and k = 2π/a for different values of TNN spring stiffness, c ′ 1 = c ′ 2 = c ′ , are plotted in Fig. 1e, f.We see that for c 1 < c 2 (c 1 > c 2 ) the winding number of ρ(k) around the origin is one (zero).
It is also instructive to find the values of k at which the Dirac points appear when To get this, we note that the band gap closes when the off-diagonal term in the matrix C(k) in Eqn. 12 is zero: where z = e ika .Clearly, there are three Dirac points when = 1 since z = e ika .Now, there are two cases: (i) is real, (ii) is complex.In the first case, we have: The solution of z in Eqn. 14 prevents η from being zero since it can blow up z.A negative η is also not possible because there is no negative spring constant.Therefore, the only valid solution is η = 1/3, resulting in all three Dirac points appearing at z = −1, i.e., k = π/a.
On the other hand, if is complex, the requirement of the existence of three Dirac points is: which is always true.Thus, if is complex, its absolute value is always 1, regardless of the value of η.Note that this condition is only valid when (1 + η)(1 − 3η) is complex, meaning η > 1/3 (since η ≥ 0).Hence, combining the complex and real solutions of z, we can identify the three Dirac points for η ≥ 1/3.The values of k at which these three Dirac points appear are: Note that in the limit of c ′ ≫ c or η → ∞, the three Dirac points appear at k = π/a and ±π/3a (the latter of which are equivalent to π/3a, and 5π/3a from 0 to 2π/a), corresponding to the band crossing locations presented in Fig. 1c.

Supercell Analysis of the Su-Schrieffer-Heeger Model
Let's consider supercell SC1 in Fig. 2a.To obtain the phonon dispersion and mode shapes of the supercell with a domain wall, Bloch boundary conditions are applied at The eigenvalues and eigenvectors of this matrix yield ω 2 and supercell mode shapes, respectively, as shown in Fig. 2.
Since matrix C(k) in Eqn.32 is not strictly chiral due to the sum of spring constants at the domain-wall mass is different from c 1 +c 2 +2c ′ , we add an additional spring with a constant of c 2 − c 1 (c 1 − c 2 at the domain-wall mass when it is connected by two soft (stiff) springs, c 1 (c 2 ), as shown in Fig. S4a (b) in Supplementary Note 2. We can then plug C(k) for such a setup into Eqn.13 to obtain C ′ (k) with strict chiral symmetry.The eigenvalues (ω 2 /ω 2 0 are henceforth symmetric about ω 2 /ω 2 0 = 0, as presented in Fig. S4c, d in Supplementary Note 2. Mode shapes of C ′ (k) are almost identical to those calculated from C(k) before adding an additional spring to the domain-wall mass, except that the sequences of the three modes in the case with strong TNNs are reordered, with the two above and below ω 2 /ω 2 0 = 0 sharing the same parity under inversion symmetry, opposite from the one located at ω 2 /ω 2 0 = 0, as shown in Fig. S4e-t in Supplementary Note 2.

Derivation of the Jackiw Rebbi Zero Modes
In this section, we demonstrate the existence of multiple TPDWSs as shown in Fig. 2d in using the JR theory.With an SIS, i.e., c 1 = c 2 , each band crossing point can be characterized by a massless Dirac theory.For c 1 > c 2 (c 1 < c 2 ), the breaking of SIS introduces a positive (negative) mass to each Dirac point.Due to the mass sign flipping at the domain boundary, one TPDWS in the bandgap, known as the JR zero mode, is expected to arise at the domain boundary for each Dirac cone [47], which explains the matching number of TPDWSs with c 1 ̸ = c 2 and band crossing points with c 1 = c 2 .The agreement of the two numbers also strongly resembles those in the quantum valley Hall effect in 2D, where the number of in-gap TPDWSs also matches that of bulk Dirac cones [12].Below, we provide a comprehensive demonstration of the existence of one TPDWS corresponding to each Dirac point and their analytical solutions characterizing the spatial decay observed in both the toy model analysis in Fig. 2 and the experimental observation in Fig. 5.

Jackiw Rebbi Mode Corresponding to Dirac Point at
ηc, and expanding matrix C ′ (k = π/a + δk) for small δk and m, we get: where σ i are Pauli matrices.In the above equation, m is the mass of the Dirac particle.Now, we create a domain wall at x = 0 with phase c 1 > c 2 in the region x < 0 and c 1 < c 2 in the region x > 0. Hence, m(x) > 0 for x < 0 and m(x) < 0 for x > 0, indicating the location dependence of mass m.Since the translation symmetry is broken due to the presence of the domain wall, we can replace δk with −i∂ x .We seek a zero frequency domain wall eigenmode ψ(x): Now, there are two cases: (i) (1 − 3η) > 0, and (ii) (1 − 3η) < 0. The eigenmodes corresponding to these two scenarios are discussed below.
, where f (x) is a scalar function, we obtain the following differential equation for f (x): where c 0 is a constant.Note that f (x) decays exponentially away from x = 0, since c(1−3η)a < 0 and m(x<0) c(1−3η)a > 0. Recall that the zero frequency domain wall mode is at k = π/a, the full expression of the mode is: (1 − 3η) < 0: Plugging the ansatz ψ(x) = f (x) 0 1 , where f (x) is a scalar function, we obtain the following differential equation for f (x): Similarly, f (x) decays exponentially away from x = 0 since m(x>0) c(1−3η)a > 0 and m(x<0) c(1−3η)a < 0. The zero frequency domain wall mode being at k = π/a leads to the full expression of the mode being: 4.4.2Jackiw Rebbi Mode Corresponding to Dirac Point at For simplicity, we will show here the existence of zero modes for η = 1 (this is what is considered in Fig. 2(e-g)), but the procedure applies to any η ≥ 1/3.For η = 1, the Dirac point is at k = ±π/2a.Away from the inversion symmetric point, when c 1 = c + m/2 and c 2 = c − m/2, expanding the matrix C ′ (k = ±π/2a + δk) for small δk and m: As in the case of k = π/a, we create a domain wall at x = 0 with phase c 1 > c 2 in the region x < 0 and the phase c 1 < c 2 in the region x > 0, implying the positiondependence of mass m(x), i.e., m(x) > 0 when x < 0 and m(x) < 0 when x > 0.
Since the translation symmetry is broken due to the presence of the domain wall, δk → −i∂ x .We seek a zero frequency domain wall eigenmode ψ(x): Plugging the ansatz ψ(x) = f (x) 1 0 , where f (x) is a scalar function, we obtain the following differential equation for f (x): where c 0 is a constant.Notice that f (x) decays exponentially away from x = 0 since m(x>0) 4ca < 0 and m(x<0) 4ca > 0. Recalling that the zero frequency domain wall mode is at k = ±π/2a, the full expression of the mode is: is only one asymmetric TPDW with a zero displacement amplitude at the domain wall mass (i.e., one parity-odd mode, SC2-E2) located at ω 2 /ω 2 0 = 0 and two symmetric ones with a non-zero-amplitude domain wall mass displacement (i.e., two parity-even modes, SC2-E1/3) with frequencies shifted up/down symmetrically about ω 2 /ω 2 0 = 0.An important point to note here is that only JR modes of the same parity can hybridize due to symmetry constraints.For example, in SC1 of Fig. S4a, the two odd modes hybridize, and one of the resulting hybridized modes shifts up in frequency, and the other one shifts down.Since the even mode cannot hybridize with the other two, it remains in the middle.Similar results also hold (albeit with parities flipped) for SC2.In either case, the JR mode at k = π always ends up mixing with one of the other two modes with the same parity with frequency shifts.Thus, these two shifted TPDWSs always present a peak at k = π in their SFT plots, as shown in Fig. S4v.On the other hand, the other mode at ω 2 /ω 2 0 = 0 not mixed with the one at k = π has SFT peaks only located at the other two Dirac points in the IBZ.
The shifting of the zero-frequency(energy) boundary/domain-wall modes (after removing the diagonal elements) to finite frequencies is not unique to the SSH model, instead, it is a generic feature that arises in most known 0D topological modes in various types of topological insulators (such as corner modes in 2D higher-order topological insulators), whose energies are also sensitive to local perturbations near the localized modes [49,50].

Topologically Protected Domain-Wall States beyond Equal
Third-Nearest Neighbors The Berry connection proposed in this work provides a generalized paradigm to predict the number of TPDWSs and their wave properties.The equal TNN scenario discussed in Results and Discussion is an example when the winding number fails, while the Berry connection succeeds in predicting the number of TPDWSs.In cases when the effect of TNN difference dominates, as in Fig. S5a in Supplementary Note 3, the Berry connection makes the same prediction of the TPDWSs as the winding number does, however with additional wave information.For example, with c ′ 1 = 3c − ∆c ′ and c ′ 2 = 3c + ∆c ′ , where ∆c ′ = 0.1c, the winding number difference between the two gauges presented in Fig. S5d is one, yet ∆B(k) shows two peaks and one valley from k = 0 to 2π, corresponding to three distinct edge modes similar to the scenario in Results and Discussion.With an enhanced c ′ difference, such as when ∆c ′ = 0.3c, we obtain three distinct peaks in ∆B(k) in IBZ at k = 1.23, π, and 5.05, as presented in Fig. S5b, c, corresponding to the three Dirac points predicted with Eqn. 17 with η = 3. Integration of these three local peaks all yields one, suggesting three TPDWSs existing in the bulk bandgap.These three TPDWSs also coincide with the winding number of three due to the all positive signs of the local peaks of ∆B(k).The contour plot of the off-diagonal element of C ′ (k) expressed in Eqn. 13 also shows winding numbers of -1 and 2 for the two gauges presented in Fig. S5e, yielding a difference of 3 between the two phases.Thus, these three TPDWSs are expected using either the Berry connection or the winding number calculation.Their existence can be confirmed by conducting a supercell analysis in the same manner described in Supercell Analysis of the Su-Schrieffer-Heeger Model.The band diagram, mode shapes, and the SFT of the TPDWS and bulk modes are presented in Fig. S5f-l.As can be seen from the SFT plots in Fig. S5l, peak locations of the TPDWSs all match those in the ∆B(k) plot presented in Fig. S5b, indicating consistent TPDWS wavelengths predicted with Berry connections.The evolution of ∆B(k) with ∆c ′ also suggests the breakdown of the winding number prediction fails when ∆c ′ < 0.2c due to the flip of the peak and valley in ∆B(k) at k = π due to a winding number difference of one (shown in Fig. S5d) as opposed to three from both the ∆B(k) calculation and the supercell analysis.Moreover, the Berry connection calculation is also applicable to lattices beyond TNN interactions.For example, with identical TNNs (c ′ ) and fifth nearest neighbors (FNNs), c ′′ , and nonidentical NNs (i.e., c 1 ̸ = c 2 ), the winding number difference is still one, as shown in Fig. S6a in Supplementary Note 3.However, the Berry connection ∆B(k) reveals three peaks and two valleys with local integrals ±1, Fig. S6b, suggesting five TPDWSs existing in the bulk bandgap.Indeed, from the supercell analysis with a similar setup as shown in Figs.2a and S5a, five edge modes emerge within the bulk bandgap, Fig. S6c, with their mode shapes presented in Fig. S6d-j.SFT of these five TPDWSs in Fig. S6k reveals that all of them are a hybridization of five wavelengths with wave numbers corresponding to the locations of peaks and valleys shown in Fig. S6b.One can prove in a similar fashion that as long as the differences in BNNs are sufficiently small, the winding numbers of two gauges will always yield zero and one, inconsistent with the actual number of TPDWSs, which can, nonetheless, be conveniently captured by the ∆B(k) calculation.

Experimental Fabrication and Characterization
The specimens are 3D-printed (Stratasys F170 FDM 3D Printer) using acrylonitrile butadiene styrene (ABS) with the following parameters: Young's modulus E=1.5 GPa, Poisson's ratio µ=0.35, and density ρ=1250 kg m −3 .As presented in Fig. 5d, each unit cell contains a pair of masses (green cubes) with side length W m =6 mm, connected by 5 mm-nearest-neighboring (NN) struts with alternating radii, r 1 =3.52 mm (blue) and r 2 =1.47 mm (yellow), to enable stiffer and softer NN interactions, respectively.Strong (weak) TNNs are established by a combination of red squared frames with side length W =16 mm, height H=4 mm (3.2 mm), and thickness t=1.33 mm (1.07 mm), and bars with radius r 3 =2.43mm (1.28 mm) connecting the masses and frames.Mode shapes of the three (one) edge modes and two bulk modes with strong (weak) TNNs modeled by COMSOL Multiphysics with these material properties and structural dimensions are presented in Fig. S7 in Supplementary Note 4. As we can see from Fig. S7a, the three TPDWSs in the lattice with strong TNNs, from high to low frequencies, are asymmetric, symmetric, and symmetric manners about the domain wall, while the one with weak TNNs, as shown in Fig. S7b is asymmetric about the domain wall.The Symmetries and locations of the deformed frames all match well with the ones obtained from experiments presented in Fig. 5i, j  To get the strict chiral stiffness matrix of the supercell, C(k), we add additional springs fixed to the ground to the domain-wall mass.Eigenvalues and eigen modes of the chiral matrix are presented in Fig. S4.   13 in the main text in the complex plane from k = 0 to 2π for d ∆c ′ = 0.1c and e ∆c ′ = 0.3c, respectively.f Band diagram (blue curves) of the supercell with spring constants listed above.Zoomed in are the three edge states (dashed blue curves) within the bulk bandgap (with bulk bands shown in red).Mode shapes of these domain-wall modes are presented in g-i with labels of SC-E(1-3).To distinguish the edge modes from the bulk ones, bulk mode shapes for bands below and above the bandgap are also plotted in j and k with labels of SC-B1/2.Red solid circles in these mode-shape plots denote the displacements of the green domain-wall mass.Although visualized in the vertical directions, all mass displacements are de f acto in the horizontal direction.Presented in l is the spatial Fourier transform of these mode shapes.

Supplementary Note 4: Finite Element Simulation of the Experimental Specimens
We conduct a finite element analysis of the experimental specimens using COMSOL Multiphysics to compare with experimental measurement.The simulated mode shapes are presented in Fig. S7. to compare with the ones obtained from experiments shown in Fig. 5i, j in the main text.

Fig. 1 c 1 Fig. 1 a
Fig. 1 Unit-cell analysis of the Su-Schrieffer-Heeger model with identical third-nearest neighbors.a

Fig. 2 Fig. 2 a
Fig. 2 Supercell analysis of the Su-Schrieffer-Heeger model with a domain wall with weak and strong identical third-nearest neighbors.SC1

Fig. 3 Fig. 3 a
Fig. 3 Spatial Fourier transform of supercell mode shapes in lattices with weak and strong identical third-nearest neighbors.

Fig. 4 Fig. 4 a
Fig. 4 Berry connection variation with third-nearest neighbor strength.a b c

Fig. 5 Fig. 5 a
Fig. 5 Experimental characterization of topological domain-wall states in lattices with strong and weak third-nearest neighbors.

1 c 1 >c 2 c 1
Fig. S1 a A polyacetylene chain created with Avogadro.b The simplified Su-Schrieffer-Heeger model consists of identical masses connected by alternating spring constants, c 1 and c 2 , with the unit cell circled in a dashed line.c Contour plots in the complex plane of ρ(k) for a complete circuit of k from k = 0 to 2π.
Fig. S4 a and b Supercells with the two arrangements of nearest neighbors with domain-wall mass connected to a fixed wall.c and d: ω 2 /ω 2 0 of with c c ′ = 1/10c and d c ′ = c.Blue and yellow bands correspond to supercell a and b, respectively.Red curves are from the unit cell analysis.Dashed blue and bold solid yellow bands are edge modes, denoted as SC1-E(1-3) and SC2-E(1-3), respectively.Corresponding mode shapes are presented in e-j and k-t.u and v: Spatial Fourier transform of mode shapes.Solid curves are the Jackiw Rebbi (JR) zero modes.

3
Fig. S5 a A supercell featuring unequal nearest neighbors with c 1 = 0.8c and c 2 = 1.2c with stiff springs (c 2 ) connected to the green domain-wall mass, and nonidentical third nearest neighbors, c ′ 1 = 3c − ∆c ′ and c ′ 2 = 3c + ∆c ′ , with strong neighbors (c ′ 2 ) connected to the green domain-wall mass and weak (c ′ 1 ) neighbors connected to the yellow domain-wall masses.b 2D and c 3D visualization of ∆B(k) from k = 0 to 2π with different ∆c ′ .d and e Contour plots of the off-diagonal element of the unit-cell stiffness matrix C ′ (k) in Eqn.13 in the main text in the complex plane from k = 0 to 2π for d ∆c ′ = 0.1c and e ∆c ′ = 0.3c, respectively.f Band diagram (blue curves) of the supercell with spring constants listed above.Zoomed in are the three edge states (dashed blue curves) within the bulk bandgap (with bulk bands shown in red).Mode shapes of these domain-wall modes are presented in g-i with labels of SC-E(1-3).To distinguish the edge modes from the bulk ones, bulk mode shapes for bands below and above the bandgap are also plotted in j and k with labels of SC-B1/2.Red solid circles in these mode-shape plots denote the displacements of the green domain-wall mass.Although visualized in the vertical directions, all mass displacements are de f acto in the horizontal direction.Presented in l is the spatial Fourier transform of these mode shapes.

b
Fig. S6 a Contour plots of the off-diagonal element of the stiffness matrix C ′ (k) in Eqn. 13 in the main text in the complex plane and b ∆B(k) from k = 0 to 2π for a lattice with identical third (c ′ = c) and fifth (c ′′ = c) nearest neighbors and nonidentical nearest neighbors c 1 = 1.2c and c 2 = 0.8c.c Band diagram (blue curves) of the supercell with spring constants listed in the figure.Zoomed in are the three edge states (dashed blue curves) within the bulk bandgap (with bulk bands shown in red).Mode shapes of these edge and bulk modes are presented in d-j with labels of SC-E(1-3) and SC-B1/2, respectively.Red solid circles in these mode-shape plots denote the displacements of the domain-wall mass.Although visualized in the vertical directions, all mass displacements are de f acto in the horizontal direction.Presented in k is the spatial Fourier transform of these mode shapes. .