To introduce our inverse design methodology for complex closed-loop locus geometries, we begin our discussion by revisiting the well-established simple 1D mechanical analog of the SSH model with identical masses, m , and alternating springs with stiffness parameters, k 1 and k 2 (Figure 1A). The governing dynamical equations describing the displacements (i.e., the first and second degrees of freedom u j and v j ) of the two masses in unit cell j , the stiffness matrix of the mechanical system, and the corresponding dispersion relation can be found in Appendix A. The topological features are primarily represented by the chiral matrix, C ( q ) , which is derived by eliminating the identical diagonal elements from the stiffness matrix (which is also the effective Hamiltonian) and is given C q = 0 − ϵ † − ϵ 0 , (1) where q is the wavenumber. The off-diagonal elements, ϵ and ϵ † , are complex conjugates; ϵ will be the focus of our discussion now. Here, ϵ is referred to as the “complex phase locus,” which is a complex function of the wavenumber q and stiffness parameters k 1 and k 2 ; it is expressed as, ϵ = k 1 + k 2 e i q . (2)

As detailed in Appendix A, the eigenvalues of C ( q ) , denoted as λ , can then be solved as: λ = ± | ϵ | / m . (3)

Results of the normalized eigenvalues Ω 2 = m λ / k (with k = k 1 + k 2 ) for q ∈ [ − π , π ] and different relationships between k 1 and k 2 are presented in Figure 1B.

We now focus on how ϵ dictates the topology of the diatomic lattice. If ϵ in Equation 2 is depicted in the complex plane for a complete cycle of the wavenumber q ∈ [ 0,2 π ] , the result is a perfect circle with a radius of k 2 and a center shifted horizontally by k 1 from the origin. This can be easily proven by transforming the expression ϵ into the form: ϵ = x + i y , (4)

with x = R ( ϵ ) ( y = I ( ϵ ) ) being the real (imaginary) part of ϵ . Using Equations 2, 4, the x and y values can be shown as: x = k 1 + k 2 ⁡ cos q , (5a) y = k 2 ⁡ sin q , (5b)

once again, representing a circle with the radius being k 2 as shown in Equations 5a, b and the parametric variable (or angle) being q . The shift of the circle’s center in the x -direction is also evident by the addition of k 1 in Equation 5a. As shown in Figure 1C, the center of the circle is at the coordinates ( k 1 , 0 ) , located on the positive real axis , since k 1 > 0 is required for dynamical stability ( k 2 > 0 is also enforced for the same reasoning). From the figure, we can count the number of times the locus ϵ winds about the origin of the complex plane—the winding number, ν —which characterizes the topological domain-boundary states of this simple SSH lattice. Such a winding number can be equivalently computed using the following integral: ν = ∫ − π π B q d q , (6) where B ( q ) is the Berry connection, defined as per Chiu et al. (2016): B q = 1 4 π i tr σ 3 C − 1 ∂ q C . (7) Here, σ 3 denotes the third Pauli matrix and tr [ ⋅ ] represents the matrix trace.

For the conventional diatomic lattice, where the complex phase locus is given by Equation 2, the Berry connection can be derived analytically. This involves evaluating the inverse and the q -derivative of the chiral matrix C ( q ) in Equation 1 and substituting them into Equation 7. The resulting expression for the Berry connection is as follows: B q = k 2 k 2 + k 1 ⁡ cos q 2 π k 1 2 + k 2 2 + 2 k 1 k 2 ⁡ cos q . (8)

By performing the integral in Equation 6 for the Berry connection (i.e., the integrand) in Equation 8 using MATLAB’s symbolic toolbox, we obtain the winding number: ν = 1 2 1 − s g n k 1 − k 2 k 1 + k 2 . (9)

The topological invariant—the winding number μ as expressed in Equation 9—can thus be summarized as follows: ν = 1 k 2 > k 1 0 k 2 < k 1 (10) As observed from Figures 1B, C, if k 1 < k 2 , the locus ϵ will wind the origin once, ν = 1 , suggesting the choice of the unit cell being topological and vice versa.

In the following sections, we will inversely design diatomic lattices with complex networks based on the shape of ϵ given by the relationship between x and y in Equation 4. Three examples of complex phase loci will be presented to demonstrate how to design a mechanical network system to achieve desired topological characteristics.

In this section, the inverse design methodology for a diatomic lattice based on an elliptical complex phase locus is detailed. Herein, the masses, m , of the inverse-designed lattices are assumed to be constant throughout this study, and emphasis is placed on the spring interconnections. Using Equation 4, the parametric form the x and y coordinates for an ellipse, with its center being on the positive real axis , are given by: x = k 1 + k 2 ⁡ cos q , (11a) y = k 3 ⁡ sin q . (11b)

Here, k 1 shifts the center of the ellipse from the origin of the complex plane, similar to its role in Equation 5a for a circular locus. The parameters k 2 and k 3 define the semi-major and semi-minor axes of the ellipse, respectively. Substituting Equations 11a, 11b into Equation 4, we obtain the following: ϵ = k 1 + k 2 + k 3 2 e i q + k 2 − k 3 2 e − i q . (12) Notice that Equation 12 reduces to Equation 2 when k 2 = k 3 , as a circle is a special case of an ellipse. The parameters k 1 , k 2 , and k 3 not only define the geometry of the ellipse but also determine the properties of the springs interconnecting the diatomic lattice.

To inversely design a lattice from the complex phase locus expressed in Equation 12, we need to understand the roles of each term constituting ϵ in conjunction with the Bloch theorem. Recall that the terms ϵ and ϵ † are the off-diagonal elements in the chiral matrix C in Equation 1, which describes the coupling between the unit cell’s degrees of freedom. For instance, ϵ represents the coupling between the second degree of freedom in the current unit cell v j and the first degree of freedom u j . The latter can be in the current cell or any other cell, expressed through the term e ± i q . Similarly, ϵ † couples the first degree of freedom to the second in any cell.

Thus, we can establish the interconnecting network between the two degrees of freedom based on the following three characteristics of Equation 12:

• Characteristic 1: The term k 1 couples the two degrees of freedom within the current unit cell j , i.e., u j and v j , as it is independent of the wavenumber q .

Based on these three characteristics, the degrees of freedom are connected as described, resulting in the diatomic lattice in Figure 2A, with the unit-cell equations of motion being: m u ̈ j + k 1 + k 2 u j − k 1 v j − k 2 + k 3 2 v j − 1 − k 2 − k 3 2 v j + 1 = 0 , (13a) m v ̈ j + k 1 + k 2 v j − k 1 u j − k 2 − k 3 2 u j − 1 − k 2 + k 3 2 u j + 1 = 0 . (13b)

We can observe that if the degrees of freedom in Equations 13a, b are condensed via the Bloch theorem, the resulting stiffness matrix will have constant diagonal elements of k = k 1 + k 2 despite the introduction of a new spring constant k 3 . Interestingly, the diagonal elements of the stiffness matrix for the diatomic lattice in Figure 2A are identical to those of a conventional diatomic lattice (Figure 1A) without next-nearest-neighbor connections, (see Appendix A). This similarity arises from the cancellation of the equal magnitude but opposite sign stiffnesses k 3 / 2 from the nearest and next-nearest couplings.

To reveal the topological properties of the diatomic lattice with an elliptical phase locus in Figure 2A, the winding number is calculated by performing the integral in Equation 6 for the chiral matrix C in Equation 1 with ϵ defined in Equations 11a, b. This process yields a closed-form expression for the winding number: ν = 1 2 1 − s g n k 1 − k 2 k 3 . (14)

As with the conventional diatomic lattice with a circular locus in Figure 1C, the winding number is defined via Equation 14, indicating that it only depends on k 1 and k 2 , provided that k 3 > 0 to ensure dynamical stability. Figures 2B, C illustrate the dispersion relation and the corresponding complex phase loci, respectively, for three cases: (i) topologically trivial k 2 < k 1 , (ii) topological transition k 2 = k 1 with the sole bandgap closed, and (iii) topologically nontrivial k 2 > k 1 . The last case can be inferred by the winding of the ellipse around the origin of the complex plane (i.e., ν = 1 as indicated by Equation 10), which is clearly seen in the right-most panel of Figure 2C.

Next, we shift our focus to a more complex family of loci known as epitrochoids, whose general form is depicted in Figure 3A. Epitrochoids are closed-loop roulettes generated by tracing a point fixed at a distance k 4 (initially pointing to the negative real axis) from the center of a circle with radius k 3 . Termed the “epicycle,” this circle rolls counterclockwise without slipping around the outside surface of a stationary circle of radius k 2 , completing a closed loop. The stationary circle is assumed to have its center on the positive real axis, with a horizontal distance of k 1 from the origin of the complex plane. Here, we define an integer n that takes a value greater than 2; its relationship to the inversely-designed lattice will be explained shortly. Additionally, depending on the epitrochoid parameters, this integer n determines the number of internal loops the epitrochoid may exhibit—specifically n − 1 internal loops.

From the general epitrochoid locus, various special cases can be derived by selecting different sets of parameters. For instance, when k 3 = k 4 and n ≥ 2 , we obtain epicycloids (Figure 3Bi). When k 3 = k 2 and n = 2 , we obtain a limaçon (also known as “Pascal’s snail”), similar to the system briefly discussed by Chen et al. (2018). For the latter, a looped limaçon (Figure 3Bii) occurs when k 4 / k 3 > 1 , transitioning to a non-looped limaçon when k 4 / k 3 = 1 , at which point the curves take the shape of an epicycloid with a single cusp (often called a “cardioid”). When 1 / 2 < k 4 / k 3 < 1 , the limaçon curve loses its inner loop and features a dimple instead (Figure 3Biii). Further reduction of the ratio below one-half smooths out the curve to a convex limaçon until k 4 approaches 0, at which point the limaçon converges to a circle (Figure 3Biv). Note that n = 1 also results in a circle of radius k 2 + k 3 − k 4 , which is not of interest in this work. As discussed next, the parameters k 1 through k 4 will be related to the spring constants of a diatomic lattice via the inverse design methodology.

Again, the inverse design of a diatomic lattice with an epitrochoid complex phase locus starts by explicitly stating x and y parametric expressions in Equation 4 for an epitrochoid, which are given by: x = k 1 + k 2 + k 3 cos q − k 4 ⁡ cos k 2 k 3 + 1 q , (15a) y = k 2 + k 3 sin q − k 4 ⁡ sin k 2 k 3 + 1 q . (15b)

To construct the complex phase locus ϵ with the epitrochoid roulette, the factor multiplied by the wavenumber in the cosine/sine arguments in Equations 15a, b must be an integer n . Therefore, k 2 / k 3 + 1 = n , mandating that k 2 = ( n − 1 ) k 3 . Thus, the parametric equations in (Equations 15a, b) reduce to: x = k 1 + n k 3 ⁡ cos q − k 4 ⁡ cos n q , (16a) y = n k 3 ⁡ sin q − k 4 ⁡ sin n q . (16b)

Substituting Equations 16a, b back into Equation 4 provides the expression of ϵ as: ϵ = k 1 + n k 3 e i q − k 4 e i n q . (17)

Following the methodology presented in Section 2.2, the diatomic lattice corresponding to the complex phase locus in Equation 17 is shown in Figure 4. Note that due to the negative stiffness k 4 coupling (as implied by Equation 17), an elastic foundation must be implemented to prevent dynamical instability. Consequently, adding an elastic foundation stiffness equal to twice k 4 is recommended. This adjustment results in the following equations of motion for the unit cell: m u ̈ j + k 1 + k 2 + k 3 + k 4 u j − k 1 v j − k 2 + k 3 v j − 1 + k 4 v j − n = 0 , (18a) m v ̈ j + k 1 + k 2 + k 3 + k 4 v j − k 1 u j − k 2 + k 3 u j + 1 + k 4 u j + n = 0 . (18b)

After utilizing k 2 + k 3 = n k 3 , Equations 18a, b become m u ̈ j + k 1 + n k 3 + k 4 u j − k 1 v j − n k 3 v j − 1 + k 4 v j − n = 0 , (19a) m v ̈ j + k 1 + n k 3 + k 4 v j − k 1 u j − n k 3 u j + 1 + k 4 u j + n = 0 . (19b)

The topology of the diatomic lattice with an epitrochoid locus depends on the choice of the three stiffness parameters k 1 , k 3 , and k 4 , as k 2 and k 3 are dependent (as discussed above). Since k 1 + n k 3 + k 4 in Equations 19a, b represents the sum of the spring constants connected to each mass in the lattice and remains constant, we define this spring constant as k = k 1 + n k 3 + k 4 . Recall that the integer n denotes the connection between one of the degrees of freedom in the current cell j and its counterpart in the ( j − n ) th cell. Therefore, the choice of k 1 , k 3 , and k 4 is constrained by a given set of k and n .

As illustrated in Figure 3, the parameter k 1 determines the center position of the complex locus, while k 2 = ( n − 1 ) k 3 governs the size of the main loop. The stiffnesses k 3 and k 4 control the size of the epitrochoid’s internal loops, if present. When the internal loops are much smaller than the main loop (i.e., when k 4 ≪ k 3 ), the system exhibits only two possible winding numbers: 0 and + 1 , depending on whether k 2 + k 3 = n k 3 is smaller or greater than k 1 . This behavior aligns with that of a simple SSH model with weak or negligible non-local interactions. As the non-local spring constant k 4 increases, the internal loops grow in size and begin to wind around the origin counter-clockwise. The winding number then assumes positive odd or even integers when n is even or odd, respectively. This phenomenon can be attributed to the symmetry of the locus about the real axis. When k 4 ≫ n k 3 and k 1 is relatively small, the size of the internal loops approaches that of the main loop, resulting in the maximum possible winding number, equal to n . A summary of possible winding numbers for an epitrochoid complex locus is shown in Table 1, and examples of winding numbers for various cases of n are provided in Appendix B.

Additionally, it is important to note that the maximum number of Dirac points in the first Brillouin zone, corresponding to the number of band crossings when the bandgap is closed, is n . This corresponds to the difference between the maximum and minimum values of ν . This observation aligns with the phonon realization of the Jackiw–Rebbi theory (Jackiw and Rebbi, 1976) in lattices with non-local interactions, as discussed by Rajabpoor Alisepahi et al. (2023).

For example, consider the case of n = 5 . Here, five distinct topological states are observed, corresponding to winding numbers of ν = 0,1,2,4 , and 5 , as depicted in the barycentric triangle in Figure 5. When the combination of these stiffness parameters varies along the transition boundaries between regions of different winding numbers, the bandgap closes when | ϵ | = 0 , indicating topological transitions. Note that the numerical procedures for calculating μ follows Pal et al. (2018).

It is worth noting that the bandgap’s width increases significantly near the corners of the Barycentric triangle, which correspond to the cases k 1 → k , n k 3 → k , or k 4 → k (Figure 5). This behavior arises because the bandgap width is strongly influenced by the degree of space inversion symmetry (SIS) breaking within a unit cell. Near the corners of the Barycentric triangle, SIS is significantly disrupted due to large disparities between the nearest neighbors (near the top two corners) and the beyond-nearest neighbors (near the bottom corner). At the center of the triangle, where the strengths of k 1 , n k 3 , and k 4 are comparable, the symmetry-breaking effects of nearest and beyond-nearest neighbors compete. This competition leads to intermediate values of the winding number ν and weakens the overall symmetry-breaking effect, resulting in a narrower bandgap. Therefore, the area of each ν in the barycentric triangle does not directly correlate with the bandgap width. Instead, the bandgap width is determined by the specific combination of k 1 , n k 3 , and k 4 , which governs the extent of SIS breaking.

To highlight the difference between the topologically distinct regions, five representative combinations of stiffness parameters are selected and listed in Table 2. Their corresponding dispersion relations, calculated using Equation 3 in conjunction with Equation 17, as well as their complex phase loci, are presented in Figure 6.

This section explores complex phase loci in the shape of hypotrochoids. These closed-loop roulettes are generated like the epitrochoids discussed in Section 2.3, with the definitions of parameters k 1 through k 4 remaining unchanged (Figure 7A). However, there are two key differences in terms of the rotation of the epicycle: (1) it rolls without slipping on the inside surface of the stationary circle in the counterclockwise direction, and (2) the distance k 4 fixed to the epicycle is initially oriented to the positive real axis. Like the epitrochoid case, an integer n ≥ 1 is defined, and a hypotrochoid may exhibit n + 1 internal loops, depending on its parameters. This integer n plays a crucial role in the inverse design of a diatomic lattice with a hypotrochoid complex locus, as will be elaborated shortly.

Under special conditions, a generic hypotrochoid complex phase locus can transform into more simplified geometries (Figure 7B). For instance, it becomes a hypocycloid when k 3 = k 4 and n ≥ 2 (Figure 7Bi) and an ellipse when n = 1 , k 2 = 2 k 3 and k 4 ≠ k 3 (Figure 7Bii). In the special case of an ellipse, the minor and major axes are equal to k 3 − k 4 and k 3 + k 4 , respectively, reinforcing that a zero k 4 results in a circle (with a radius of k 3 ) (Figure 7Biii). Additionally, the ellipse collapses to a line when k 3 = k 4 , resulting in a looped line, often referred to as the “Tusi mechanism” (Figure 7Biv).

Following Equation 4, a hypotrochoid roulette in the complex plane is defined via the following x and y parametric functions: x = k 1 + k 2 − k 3 cos q + k 4 ⁡ cos k 2 k 3 − 1 q , (20a) y = k 2 − k 3 sin q − k 4 ⁡ sin k 2 k 3 − 1 q . (20b)

As in the case of the epitrochoid, the factor multiplied by the wavenumber in the cosine/sine arguments in Equations 20a, b must be an integer. Therefore, k 2 / k 3 − 1 = n (or k 2 = ( n + 1 ) k 3 ), which, in combination with Equations 4, and 20a, b yields the following complex phase locus: ϵ = k 1 + n k 3 e i q + k 4 e − i n q . (21)

After constructing the diatomic lattice corresponding to ϵ in Equation 21, as shown in Figure 8, the following equations of motion for the unit cell are derived: m u ̈ j + k 1 + k 2 − k 3 + k 4 u j − k 1 v i − k 2 − k 3 v j − 1 − k 4 v j + n = 0 , (22a) m v ̈ j + k 1 + k 2 − k 3 + k 4 v j − k 1 u j − k 2 − k 3 u j + 1 − k 4 u j − n = 0 , (22b)

which, after using k 2 − k 3 = n k 3 , simplify to: m u ̈ j + k 1 + n k 3 + k 4 u j − k 1 v i − n k 3 v j − 1 − k 4 v j + n = 0 , (23a) m v ̈ j + k 1 + n k 3 + k 4 v j − k 1 u j − n k 3 u j + 1 − k 4 u j − n = 0 . (23b)

Topological transitions in lattices featuring hypotrochoid loci are achieved by tuning the stiffness parameters k 1 , k 3 , and Equations 23a, 23b k 4 . As in the case of epitrochoid loci, k 1 determines the center position of the complex locus, k 2 = ( n + 1 ) k 3 governs the size of the main loop, and k 3 and k 4 control the size of the internal loops, if present. Mirroring the epitrochoid loci, when the internal loops are significantly smaller than the main loop (i.e., when k 4 ≪ k 3 ), the winding number is either 0 or +1, depending on whether k 2 − k 3 = n k 3 is smaller or greater than k 1 . As the non-local spring constant k 4 increases, the internal loops may begin to wind around the origin in a clockwise direction, and the winding number assumes negative odd or even integers, matching the parity of the integer n . This progression continues until k 4 ≫ n k 3 , at which point the size of the internal loops approaches that of the main loop. Assuming that k 1 remains relatively small, the system reaches its maximum possible winding number, ν = − n . A summary of possible winding numbers for hypotrochoid complex loci is provided in Table 3, and examples of winding numbers for various cases of n are detailed in Appendix B.

Similar to the epitrochoid complex phase loci, the maximum number of Dirac points in the first Brillouin zone is n , matching the difference between the maximum and minimum values of ν .

Using n = 5 as an example, we numerically calculate the winding numbers for all combinations of stiffness parameters within the constraint k = k 1 + n k 3 + k 4 . As with the epitrochoid complex phase loci, we identify five distinct winding numbers of ν = 0 , + 1 , − 1 , − 3 , and , − 5 , separated by boundaries where the bandgap closes (precisely at | ϵ | = 0 ), as shown in the barycentric triangle in Figure 9. Like the epitrochoid case, when the combination of stiffness parameters relatively small along the transition boundaries between regions of different winding numbers, the bandgap closes, signaling topological transitions. Furthermore, the bandgap reaches its maximum size when k 1 , n k 3 or k 4 approaches k (Figure 9). Five representative combinations of stiffness parameters (listed in Table 2) are selected to showcase their corresponding dispersion relations, winding number, and complex phase loci (Figure 10).