Topological synchronization of coupled nonlinear oscillators

Synchronization of coupled oscillators is a ubiquitous phenomenon found throughout nature. Its robust realization is crucial to our understanding of various nonlinear systems, ranging from biological functions to electrical engineering. On another front, in condensed matter physics, topology is utilized to realize robust properties like topological edge modes, as demonstrated by celebrated topological insulators. Here, we integrate these two research avenues and propose a nonlinear topological phenomenon, namely topological synchronization, where only the edge oscillators synchronize while the bulk ones exhibit chaotic dynamics. We analyze concrete prototypical models to demonstrate the presence of positive Lyapunov exponents and Lyapunov vectors localized along the edge. As a unique characteristic of topology in nonlinear systems, we find that unconventional extra topological boundary modes appear at emerging effective boundaries. Furthermore, our proposal shows promise for spatially controlling synchronization, such as on-demand pattern designing and defect detection. The topological synchronization can ubiquitously appear in topological nonlinear oscillators and thus can provide a guiding principle to realize synchronization in a robust, geometrical, and flexible way.


I. INTRODUCTION
Nonlinear oscillators ubiquitously appear in a variety of fields from biology [1,2] to engineering [3][4][5][6]. They often exhibit frequency-synchronization [7], where (even inhomogeneous) interacting oscillators vibrate at the same frequency. Synchronization plays a crucial role in classical nonlinear systems, e.g., circadian rhythm [8], and even in quantum systems [9,10]. Controlling synchronization and desynchronization [11][12][13] of nonlinear oscillators is essential to fulfill their functions. However, nonlinear oscillators are often irregularly affected by their uncontrollable circumstances, and therefore it is desirable to investigate universal principles to robustly control and design nonlinear oscillators against such disorder.
On another front, topology represents the property of matter unchanged under continuous deformations and thus can provide a guiding principle to realize robust systems. Numerous studies in condensed matter physics have focused on such utility of topology, where a remarkable example is a topological insulator [14][15][16] exhibiting a metallic surface and an insulating bulk. Topological insulators exhibit topologically nontrivial bulk associated with edge-localized modes robust against disorder, as a consequence of a principle called bulkedge correspondence [15,16]. Such edge modes exhibit gapless dispersion relations and backscattering-free edge current. Because bulk topology is intrinsically related to Hamiltonians that govern linear dynamics, topology in physics has been studied mostly in linear systems. In recent studies, it has been revealed that topological edge modes can affect the nonlinear dynamics [17][18][19] in, e.g., photonics [20] and mechanical lattices [21], where topological edge soliton [22] and bulklocalized topological modes [23] can emerge. Despite these recent advancements, topology in nonlinear systems is still largely unexplored even at the conceptual level. * sone@noneq.t.u-tokyo.ac.jp Given these fundamental problems, we here propose a topological mechanism of robust edge-localized frequencysynchronization, where only nonlinear oscillators at the edge of the system synchronize, whereas dynamics of the bulk oscillators is chaotic. We demonstrate the emergence of such a synchronized state by numerical calculations of nonlinear oscillators with linear couplings corresponding to a topological Hamiltonian. We term the proposed state as a topological synchronized state (TSS).
By calculating the Lyapunov exponents [24] and vectors [25,26], we confirm the chaotic behavior of bulk oscillators and find that topological edge modes are represented by the edge-localized Lyapunov vectors. While such coexistence of synchronization and desynchronization is apparently reminiscent of chimera states [27][28][29][30][31][32][33][34], we emphasize that the mechanism of the proposed TSS is distinct from conventional chimera states because of its topological origin. Moreover, nonlinearity leads to another unique topological phenomenon: unconventional extra boundary modes appear at the emergent effective boundary. Such extra boundary modes can increase the number of synchronized oscillators in the TSS.
As applications of the TSS, we propose to arrange the synchronized oscillators in an on-demand pattern, and to detect defective structures by observing the lasing patterns around the defects. We also propose a concrete electrical-circuit realization of the TSS. While we focus on several types of non-Hermitian Hamiltonians [35][36][37][38][39][40] of lasing edge modes [41][42][43] to realize the TSS, we expect that the TSS can ubiquitously emerge in nonlinear oscillators with topological linear couplings. In addition, the robustness of the topological edge modes guarantees the stability of the TSS against detuning parameters and thus provides a robust and universal mechanism to control synchronization.

II. CONCEPT OF TOPOLOGICAL SYNCHRONIZATION
Before moving to detailed results of numerical analyses, we illustrate the general concept of TSS proposed in this article. Here, we consider coupled nonlinear oscillators that self-vibrate even if they are decoupled. Such oscillators maintain their self-oscillations at the balance of the injection and dissipation of energy and exhibit nonlinear features that cannot arise in linear harmonic oscillators, such as the robustness of the amplitude. Dynamics of nonlinear oscillators can be described as a limit cycle; one of the simplest models for this is the (complex) Stuart-Landau oscillators [44] described aṡ Z = (iω + α − β|Z| 2 )Z, where Z is the complex-valued state variable. In the following, we focus on the case that β is real. This model exhibits self-excited oscillations with frequency ω and amplitude α/β and is the normal form of a nonlinear oscillator around the parameters where it begins self-oscillations.
TSS is defined as a coexistence of synchronized edge oscillators and desynchronized bulk oscillators due to the nontrivial topology of the system. Such synchronization of the edge os-cillators can appear in coupled nonlinear oscillators where the linear coupling is described by a Hamiltonian that is topological in the wavenumber space. We can construct concrete models of the TSS by utilizing linearly coupled Stuart-Landau oscillators arranged on a square lattice (see Fig. 1(a)). In the present work, we assume that four oscillators exist at each site and they interact with ones at the same site as well as the nearest neighbor sites (four or more oscillators are necessary to realize the desirable properties of the linear coupling, such as its lasing edge modes discussed in the following). The dynamics of our model is described by the following set of equations, where x = (x, y), x = (x , y ) represent the location of the site in the square lattice, and j, k = 1, · · · , 4 distinguish four oscillators at each site. H jk (x, x ) represents the linear coupling between the oscillators at the same site or the nearest neighbor sites. We assume Z k (x ) = 0 for x corresponding to the outside of the system, which realizes open boundaries. We note that this choice of the boundary condition is just for the simplicity of the models, and TSS is independent of the detail of boundary conditions on the condition that both edges are disconnected. We adopt a Hamiltonian of topological lasing modes [41][42][43] as the matrix H jk (x, x ) describing the linear coupling, and obtain TSS under the open boundary condition.
To investigate the robustness of the TSS against disorders, we introduce the inhomogeneity into the system such that natural frequencies of each oscillator, ω j (x), are uniformly distributed around the mean value ω 0 with the distribution width being ∆ω.
In the following sections, we analyze three models of TSS. These models are constructed in the same spirits as described above, while the linear couplings are different. In our first model (Sec. III), we introduce the Hamiltonian featuring the exceptional edge modes [43] to describe the linear coupling. Exceptional edge modes are protected in an unconventional mechanism unique to non-Hermitian systems and enable us to realize TSS in a simple system. To show that we can also realize TSS by using conventional topological edge modes, we construct another model in Sec. IV. In Sec. V, we propose another model utilizing Hermitian linear couplings by modifying our first model. The properties of our three models are summarized in Table I. In Sec. VI, we propose applications of TSS to on-demand pattern designing and defect detection with their numerical demonstration. We also present a schematic of the possible realization of TSS using an electrical circuit. Section VII summarizes the main results and discusses several open problems.

A. Model and numerical calculation of its dynamics
To realize the TSS in a simple model, we first adopt the Hamiltonian of exceptional edge modes [43] as H jk (x, x ), which utilizes the robustness of singularities called exceptional points, topological gapless structures unique to non-Hermitian systems. The exceptional edge modes are protected by the topology of the dispersion relations around their branchpoint structures [35] and exhibit lasing behavior [41,42] amplifying edge oscillations. In nonlinear oscillators, such edge modes can lead to qualitatively different synchronization behavior between edge and bulk oscillators without judicious designing of the system.
In more detail, we construct the linear coupling starting from Phase diagram of the model of the TSS. The ratio of frequency fluctuations of the edge oscillators to those of the bulk ones is depicted at each parameter. We calculate the dynamics of the model of the TSS (Eqs. (1), (2)) under different strengths of the linear coupling and the inhomogeneity of the natural frequencies.
The red region (right bottom) shows parameter regimes where the TSS emerges. The circles and squares represent the classification of data of frequency fluctuations of the edge and bulk oscillators obtained by a numerical clustering method. The circles correspond to the parameters for which the presence of the TSS is identified. The phase boundary obtained here corresponds to the parameter points at which the bandgap of the effective Hamiltonian closes and thus the bulk band can change its topology. the Hamiltonian of a Chern insulator with two internal degrees of freedom, H QWZ = (u + cos k x + cos k y )σ z + sin k x σ x + sin k y σ y , called the Qi-Wu-Zhang (QWZ) model [45]. This corresponds to the lattice version of the Dirac Hamiltonian, which is a prototype of topological insulators and exhibits edge-localized modes associated with its topological invariants (Chern numbers). Then, we combine the Chern insulator and its time-reversal counterpart H * QWZ by the non-Hermitian coupling ibσ x ⊗ σ x and obtain the linear coupling that exhibits amplification of edge modes protected by exceptional points. Such Hamiltonian used in our prototypical model is expressed in the wavenumber space as where I 2 is the 2 × 2 identity matrix and σ x,y,z are the Pauli matrices We can obtain the following real-space description by the inverse Fourier transformation: where δ x,x represents the Kronecker delta and e x,y are the lattice vectors in the x or y direction. ibσ x ⊗ σ x represents the non-Hermitian coupling between the QWZ model and its time-reversal counterpart. We note that there are coupling terms with imaginary coefficients, and these can be realized by doubling the internal degrees of freedom or time-delayed coupling as seen in electrical circuits [46]. We numerically calculate the dynamics of our model under the open boundary condition and confirm the emergence of the TSS, i.e., the synchronization (desynchronization) of the oscillators at the edge (in the bulk). Figures 1(c),(d) show the frequency distribution obtained from the numerical simulations. One can see that the edge oscillators exhibit constant and homogeneous frequencies, indicating their frequency-synchronization, while the bulk ones oscillate at time-and space-varying frequencies. It is noteworthy that the fluctuations of natural frequencies are irrelevant to the TSS (see Appendix C for the related numerical calculation). The existence of TSS under the inhomogeneous natural frequencies indicates the robustness of TSS guaranteed by its topological nature.
We carefully examine the robustness of TSS by calculating the dynamics of our model at different strengths of the linear coupling and inhomogeneity of the natural frequencies. Figure 2 shows the phase diagram determined on the basis of frequency fluctuations of the edge and bulk oscillators, which signifies the presence/absence of the TSS. If the inhomogeneity of the natural frequencies is large and/or the linear coupling is weak, we obtain a fully desynchronized state. We note that the parameter where the TSS disappears corresponds to the point at which the dispersion relation of the Hamiltonian becomes gapless due to strong disorders and thus can change its topology. Therefore, the result in Fig. 2 supports the argument that TSS is topologically protected. We can also analyze the robustness of TSS by constructing a one-dimensional model imitating the dynamics of the synchronized edge oscillators in our model (see Appendix D and Eqs. (D1), (D2)).
Here, we emphasize that the TSS is defined as the coexistence of synchronized edge oscillators and desynchronized bulk oscillators with a topological origin. If we set α negative in our model (1), we can damp the bulk oscillations (see Appendix B), while the edge ones still exhibit synchronized oscillations (note that a similar synchronized state is observed in previous research [47]). However, such synchronization is not a TSS defined here because the bulk oscillators do not exhibit (chaotic) self-oscillations. Using another Hamiltonian as linear couplings, we can also realize cluster synchronization where the edge oscillators and bulk ones oscillate at different frequencies (see Appendix B), which is also out of the range of the TSS.

B. Chaos in bulk oscillators and the edge-localized Lyapunov vectors
We next show that the bulk oscillators exhibit chaotic dynamics. Chaos is characterized by the extreme sensitivity to initial conditions, i.e., exponential growth of the initial difference between trajectories. The infinitesimal rate of such exponential divergence is called the Lyapunov exponent [24], and its positivity gives a defining feature of chaos. We numerically calculate the Lyapunov exponents of the first model. Figure  3(a) shows the obtained ones and the existence of positive Lyapunov exponents indicating the chaos. We note that Stuart-Landau oscillators can exhibit two types of chaos, namely, amplitude chaos and phase chaos [48]. In this model, we can see amplitude-chaotic behaviors, e.g., phase slips via zero amplitude (see Appendix F).
Lyapunov exponents also tell us the effective degree of freedom of chaotic attractor, Lyapunov dimension defined as D L = i≤M λ i /|λ M +1 |+M , where λ i represents Lyapunov exponents arranged in descending order, and M is the smallest integer that satisfies i≤M +1 λ i < 0. Previous studies [30,31] have revealed that the Lyapunov dimension corresponds to the number of the desynchronized oscillators in the chimera state. In our first model, we obtain the Lyapunov dimension D L = 254.568, which is almost 90 percent of the total degree of freedom of oscillators except for the first and second oscillators from the edge, D bulk = 288. Our Lyapunov analysis thus clearly indicates that the bulk oscillators are desynchronized and exhibit chaotic dynamics.
Lyapunov vectors can be used to know the local geometric information of chaotic attractors [26]. Specifically, per- turbation parallel to a Lyapunov vector is amplified or attenuated in a forward and backward process expressed as |δ Z(±t)| ∼ |δ Z(0)|e ±λit , where λ i is the associated Lyapunov exponent. In our first model, we reveal that the Lyapunov vectors corresponding to the small Lyapunov exponents are localized at the edge of the system, thus realizing topological edge modes in nonlinear systems. To clearly show such localization to the edge, we define the following index of the proportion of amplitude of the edge oscillators, where v i is the ith component of the Lyapunov vector and we sum the squares of the components corresponding to the edge oscillators. This index takes a large value when a Lyapunov vector is localized furthermost to the edge. Figure 3(b) shows the proportions of the edge oscillators of the Lyapunov vectors in our first model. One can see the steep increase in the index, which indicates that the first about 80 percent of Lyapunov vectors are extended to the bulk, while the others are localized to the edge. We demonstrate such localization and delocalization of the Lyapunov vectors by plotting them in the real space (see Appendix G). It is noteworthy that the Lyapunov exponents decrease around the index where the edge proportions increase. Thus, we expect that perturbation to the edge oscillators is attenuated, which implies stability of their synchronization.
We also evaluate the degree of localization in the entire Lyapunov vectors by calculating inverse participation ratios (IPRs) defined as The IPR becomes large when a small number of components of a Lyapunov vector exhibit large amplitudes, that is, it is localized at only a few sites. Figure 3(c) shows the IPR of the Lyapunov vectors obtained in our first model. There is a slight increase in the IPR around the relative index of the Lyapunov vectors 0.8 as a result of the localization to the edge. The smallness of the increase in the IPR indicates that most of the edge-localized Lyapunov vectors spread at many edge sites, which is reminiscent of conventional topological edge modes in linear systems. We can find another region where the IPR steeply increases. Those large IPRs imply the strong localization of the Lyapunov vectors to a few edge sites. We note that the Lyapunov exponents also decrease in this region, indicating the strong damping of the perturbation corresponding to these localized Lyapunov vectors. While the proportions of the edge oscillators in Eq. (5) reveal the existence of the edge-localized Lyapunov vectors, IPRs allow us to classify such edge-localized Lyapunov vectors into two groups, the extended ones corresponding to conventional edge modes and the strongly localized ones induced by nonlinearity.
The edge-localized Lyapunov vectors and their negative Lyapunov exponents are related to the edge modes of the effective Hamiltonian obtained via the linearization of the equation. We note that if there are no nonlinear terms (i.e., β = 0 in our model), Lyapunov exponents and vectors are identical to the imaginary parts of the eigenvalues and eigenvectors of the effective Hamiltonian, respectively. It is noteworthy that non-Hermiticity is essential to realize nonzero imaginary parts of eigenvalues, and nonlinear oscillators ubiquitously exhibit such dissipative and injective effects. In nonlinear systems, we can obtain the effective Hamiltonian by linearizing the equation around the state at each moment. The linear equation governed by the sequence of such effective Hamiltonians describes the time evolution of the difference between the perturbed trajectory and the original one. We numerically When the edge oscillators are not fully amplified (corresponding to t = 10), the edge-localized loss is not large enough to create an effective boundary. After the edge oscillators are fully amplified, the inside oscillators begin to largely oscillate, which implies the appearance of the extra topological boundary modes.
check that the effective Hamiltonian of our first model exhibits edge modes with negative imaginary parts of the eigenvalues (see Appendix E), which indicates the (short-term) stability of the synchronization of the edge oscillators. We note that nonlinear terms lead to a random on-site loss in the effective Hamiltonian. However, topological modes are robust against such disorders and thus still appear from the present effective Hamiltonian.
We also note that some previous studies [49,50] discuss the role of the localization of eigenvectors of the linearized equation in the localized patterns. However, the origins of the localization are different between such localized patterns and the TSS. The localized patterns in the previous research rely on the real-space configuration of the oscillator networks and are explained from the perturbation analysis, while the TSS utilizes its nontrivial topology in the wavenumber space. It is also noteworthy that the TSS cannot be observed in topologically trivial systems (cf. Appendix I), which indicates the crucial role of the wavenumber topology in the TSS.

C. Emergence of extra topological modes by nonlinearity-induced boundary
We find that nonlinearity can induce unconventional extra edge modes, which should be the origin of the synchronization of the second oscillators from the edge as discussed in the previous section. Namely, nonlinear systems exhibiting selfexcited vibration, including our model, can generate effective boundaries and extra topological modes regardless of the initial conditions. Previous studies [23,51] have discussed that one can create effective boundaries by externally preparing a localized initial state, whereas effective boundaries in our first model spontaneously appear by utilizing the lasing edge modes.
To show the existence of such extra boundary modes, we rewrite the equation of our model (1)  amplify the edge oscillators as shown in Fig. 1(e). We note that recent research [39] has shown that on-site loss creates an effective boundary, and topological boundary modes can appear at the boundary between the regions with gain and loss. While the previous study has focused on linear dynamics with an externally introduced boundary, in our model the on-site loss is an emergent feature in the sense that it is spontaneously induced by the nonlinearity.
To explicitly demonstrate the presence of such extra boundary modes, we numerically diagonalize the state-dependent HamiltonianH with edge-localized loss. Here, we consider the on-site loss that should be induced by lasing edge modes of the original Hamiltonian H used in Fig. 1. Figure 4(a) shows the dispersion relation of this state-dependent Hamiltonian. We obtain eight gapless modes, which are the twice number of gapless modes compared to the original Hamiltonian without on-site loss (see Appendix H for the dispersion relation of the original Hamiltonian). We find positive imaginary parts of eigenvalues and the localization to the second site of the corresponding eigenvectors, which leads to the amplification of the second oscillators from the edge.
We also directly confirm the emergence of the extra boundary modes from the numerical simulation. We set the coefficient of nonlinear terms to be small compared to those of linear coupling H. Figure 4(c) presents the snapshots of the relative amplitude of each site obtained from the simulation. In the beginning, only the outermost oscillators have large amplitudes, while after a sufficiently long time the inner ones also begin to oscillate with large amplitudes. This behavior represents the emergence of the extra boundary modes after the nonlinearity-induced loss grows enough to balance with the linear coupling and create an effective boundary.
The presence of the extra boundary modes above can alter the number of the synchronized oscillators in the TSS. Edgelocalized on-site loss also exists in the effective Hamiltonian obtained from the linearization of the equation around the state at each moment. Therefore, the number of topological modes can increase in such an effective Hamiltonian via the same mechanism as in the state-dependent Hamiltonian. We confirm the existence of the extra boundary modes and negativity of the imaginary parts of their eigenvalues (see Appendix E), which leads to the increase in the number of dissipative Lyapunov vectors.

A. Model and its dynamics
While we utilize the Hamiltonian featuring exceptional edge modes in our first model (Eq. (2) and Fig. 1), we can also realize TSS by using a Hamiltonian of topological edge modes protected by conventional bulk topology. To demonstrate this, we consider another non-Hermitian Hamiltonian exhibiting lasing edge modes, which is described in the real space as where a, b, u, and u are real parameters, and σ i , δ x,x represent the ith component of the Pauli matrices and the Kronecker delta. As in our first model, this Hamiltonian is constructed from two layers of the QWZ model (a Chern insulator model), H QWZ and aH QWZ combined by the non-Hermitian coupling ibσ x ⊗ I 2 . Thus, we can rewrite the Hamiltonian as where ±iu I 2 represent the additional gain and loss. The modification of the hopping amplitude (parametrized by a) and non-Hermitian couplings enable the amplification of only the edge modes derived from the nontrivial topology of H QWZ . Figures 5(a),(b) show the frequency distributions of nonlinear oscillators with such topological linear couplings. We can confirm the emergence of the TSS, i.e., the edge oscillators are frequency-synchronized while the bulk ones are chaotic. We can also confirm the amplification of the edge oscillators in Fig. 5(c) as in Fig. 1(e).

B. Lyapunov analysis
We also calculate the Lyapunov exponents and vectors to confirm the chaos of the bulk oscillators and the existence of the edge-localized Lyapunov vectors (see Appendix A for the detailed numerical method). Figure 6(a) shows the Lyapunov exponents of the model considered here. There are positive Lyapunov exponents which indicate the chaos of the bulk oscillators. Figure 6(b) shows the proportions of the edge amplitudes (5) of the Lyapunov exponents. We can confirm the existence of the edge-localized Lyapunov vectors as in the model using exceptional edge modes [43] (cf. Eq. (2) and Fig. 3(b)). However, the behavior of the inverse participation ratios (IPRs) (6) shown in Fig. 6(c) is different from that in Fig. 3(c), because there are no steep increases in the IPRs. Thus, the strongly localized Lyapunov vectors disappear from the model considered here. This is because the nonlinear effects are not large enough to generate such strongly localized modes. We also note that nonlinearity-induced boundary modes are absent in this model as can be seen in Fig. 5. This may also be due to the smallness of nonlinearity.

A. Model and its dynamics
We also find that the TSS can emerge in nonidentical nonlinear oscillators with linear couplings described by a Hermitian Hamiltonian. To construct a model of such TSS, we consider the unitary transformation of the Hamiltonian considered in our first model (2) and obtain where I 2 is the 2 × 2 identity matrix and σ x,y,z are the Pauli matrices. This Hamiltonian is described in the real space as where δ x,x is the Kronecker delta and e x,y are the lattice vectors in the x or y direction. Here, we perform the unitary transformation of the Hamiltonian in our first model (2) to convert the non-Hermitian term ibδ x,x (σ z ⊗σ z ) jk into a diagonal one. Such diagonal non-Hermitian terms can be introduced as the modulation of the parameter α in the Stuart-Landau equation (cf. Eq. (1)). Thus, we can consider Stuart-Landau oscillators with this linear coupling as nonidentical oscillators whose linear coupling is Hermitian.
We numerically calculate the dynamics of such coupled nonidentical oscillators. Figure 7(a) shows the frequency distribution of oscillators at the first components of lattice points. We confirm the TSS, i.e., the edge oscillators synchronize while the bulk ones exhibit chaos. Unlike the other models of TSS, however, the second components of the oscillators exhibit chaotic motion even at the edge of the system (Fig. 7(c)). Such synchronization and desynchronization of the edge oscillators can be confirmed from the indicator of the local order parameter restricted to the first and fourth (second and third) components as shown in Fig. 7(b),(d). The first components of edge oscillators exhibit large indicators of the local order parameter indicating their synchronization, while those of the second components have small values.

B. Lyapunov analysis
We also calculate the Lyapunov exponents and vectors to confirm the chaos of the bulk oscillators as in Figs. 3, 6. Figure 8(a) shows the existence of positive Lyapunov exponents, which implies the chaos in this system. We compare the proportion of amplitude of the edge oscillators of the Lyapunov vectors (5) between the first and fourth components and the second and third ones, as shown in Fig. 8(b),(c). The Lyapunov vectors corresponding to small Lyapunov exponents are localized at the first and fourth components of the edge sites. On the other hand, the amplitudes of Lyapunov exponents at the second and third components are small at the large indices. Thus, the first and fourth components of edge oscillators are synchronized, while the other oscillators exhibit chaos.
Desynchronization of the second and third components of the edge oscillators is explained from the amplitude distribution of the lasing edge modes shown in Fig. 9. The lasing edge modes with the maximum imaginary part of the eigenvalue exhibit zero amplitude at the second and third components of the lattices and edge-localized amplitudes at the first and fourth components. Therefore, the lasing edge modes obtained from the linear coupling (9) only amplify the first and fourth components of the edge oscillators. Since such amplification leads to the different behavior between the edge and bulk oscillators (see also Appendix H), only the first and fourth components of the edge oscillators can synchronize, while the second and third components exhibit chaotic oscillations.

VI. APPLICATIONS OF TSS
Damped bulk oscillators in the TSS will create an effective edge, along which the oscillators can be newly synchronized. We here propose several applications of such emerging synchronized oscillators, e.g., on-demand synchronization pattern designing and defect detection. First, by placing damped oscillators into the bulk, one can arrange synchronized oscillators in an arbitrary pattern. We demonstrate such flexible pattern designing in Fig. 10 (see Appendix A for the numerical method to introduce damped oscillators). In Fig. 10(b), we define the indicator of the local order parameter as which becomes large when the phase of the oscillator matches those of the nearest neighbors (we note that this includes information of all the components unlike that defined in Eq. (11)). The oscillators around damped ones exhibit the homogeneous frequencies and the large values of the indicator M , thus being synchronized. We note that the synchronization pattern disappears in a topologically trivial system (see Appendix I), demonstrating the crucial role played by the topology of the linear coupling.
If some of the oscillators suddenly stop their selfoscillations, the oscillators around them begin to synchronize. Thus, by observing the appearance of a synchronized cluster, we can continuously monitor the formation of the broken oscillators, i.e., defects. To demonstrate such possibility of defect detection using the TSS, we numerically calculate the dynamics of our first model and damp some oscillators in the middle of the simulation. Figure 11(c) shows the spatial distribution of the indicator M (x, y) (12) in the steady-state regime. We find the synchronized oscillators that precisely lie along the edges of the defects. In contrast, when the sensors that detect the state of oscillators get broken and are unable to send signals to the central computing system (see Fig. 11(b)), our method obtains no enhanced signals indicating synchronization. Therefore, we can affirm the existence of broken oscillators from the large value of the indicator M (x, y) (12) even under the nonnegligible possibility of this kind of breakdown of sensors. We numerically check the absence of the signal of synchronization in the case of sensors' disorder (see Fig. 11(d)). Meanwhile, a system similar to our model can be realized by using an electrical circuit, which is shown in Fig. 12 as a schematic. Using nonlinear resistors, we propose to use van der Pol circuits [52,53] that simulate nonlinear oscillators. We also utilize negative impedance converters with current inversion [54] to omit the frequency dependence. The dynamics of voltages in this electrical circuit imitates our first model, while the dynamics of van der Pol circuits is different from the Stuart-Landau oscillators (see Appendix J for the derivation of circuit equations). We note that one can realize the TSS without fine-tuning of the circuit constant of each element due to the topological protection of the TSS (see Appendix D). The electrical-circuit implementation of topological synchronization can have a potential application in information processing, as is discussed in the previous studies of topological electrical circuits [19,55]. Our model should also be relevant to broad ranges of other physical systems. Especially, we expect that it is possible to realize the TSS in photonic systems, where topology [22,23,41], non-Hermiticity [39,56,57], and nonlinearity [3,4] can be intertwined.

VII. SUMMARY AND DISCUSSIONS
We proposed a robust mechanism to control synchronization, namely, topological synchronization where edge oscillators synchronize while bulk ones exhibit chaotic dynamics. We termed such coexistence state of chaotic bulk motions and synchronized edge oscillators as the topological synchronized state (TSS) and analyzed the localized behaviors of its models. The edge modes of the topological linear couplings between oscillators are represented by the Lyapunov vectors localized at the edge of the system. The mechanism relies on topology in the wavenumber space and thus is different from the conventional mechanism of chimera states; the latter can be explained from the self-consistency analysis and spatial variation of local order parameters [27,28]. We also demonstrated that nonlinearity induces emerging effective boundaries, which leads to extra topological modes unique to nonlinear systems. As examples of applications, we showed that one can realize arbitrary patterns of synchronized oscillators by using topological synchronization. We also revealed that chaotic nonlinear systems can be used to detect disordered oscillators. Our model can be realized by utilizing an electrical circuit. Our proposal provides a general method to robustly and geometrically control the nonlinear oscillators by combining nonlinearity and topology.
While in our numerical calculations, we focused on 10 × 10 and 20 × 20 lattices, TSS is independent of the system size, as long as the system is large enough to apply the band theory of bulk systems. This is because we here consider topology in the wavenumber space, and such analyses in the wavenumber space describe the bulk properties of infinite systems. Therefore, nontrivial topology and associated phenomena do not alter if we consider the larger lattice than analyzed in our calculations. We also used the Dirichlet-type open boundaries, while TSS appears under the other open boundary conditions (e.g. the Neumann boundary conditions), which may play an important role in some typical nonlinear systems such as fluids. We can understand this from the conventional argument in topological physics [15]. Such open boundaries connect topological systems and the vacuum, which have different topological invariants. However, we can change topological invariants only by closing the band gaps. Since the system and the vacuum has gapped bulk, such gap closing is possible only at the boundary. Therefore, boundary-localized gapless modes must appear in this situation. TSS shares the same origin with conventional topological edge modes and thus can appear under any type of open boundary conditions. It is noteworthy that the TSS should ubiquitously appear in nonlinear systems with a topological effective Hamiltonian. In this paper, we confirmed the emergence of TSS in ho- mogeneous Stuart-Landau oscillators whose linear couplings are described by Hamiltonians of two types of lasing edge modes protected by different topological mechanisms. We also discussed the TSS utilizing a Hermitian topological linear coupling and inhomogeneity of oscillators (cf . Table I). Furthermore, there is no need to judiciously adjust the parameters as discussed in Appendix D, which should be of practical advantage to realize the TSS in a variety of systems. Despite the broadness of models considered in this work, the TSS in identical oscillators induced by a Hermitian topological Hamiltonian still remains an intriguing problem. This may be possible by utilizing nonlinearity-induced non-Hermitian effects and realizing positive imaginary parts only in bulk modes of the effective Hamiltonian obtained by linearization of the equation (cf. Appendix E). Furthermore, we anticipate that other nonlinear oscillators than the Stuart-Landau model (e.g., the Fitzhugh-Nagumo model [58,59]) can exhibit the TSS, because topological edge modes should be stable against the change of the on-site loss induced by the nonlinear terms. the time evolution of the complex-valued state variables, we deal with real and imaginary parts of those variables separately as Z i (x, y) = X i (x, y) + iY i (x, y). Real and imaginary coefficients are substituted as aZ → a(X, Y ) T and ibZ → −ibσ y (X, Y ) T , where a, b are real, and σ y is the Pauli matrix. We arrange the 20 × 20 sites which have four complex Stuart-Landau oscillators for each. We impose the open boundary condition both in the x and y direction. We set the natural frequency of each oscillator ω 0 + Ω with Ω being a random value from uniform distributions ranging [−∆ω, ∆ω]. The numerical calculation starts from the random initial condition, where real and imaginary parts of the initial state variables are the random value from uniform distributions ranging [−0.1, 0.1]. In Fig. 1, we set the time step dt = 0.005 and use the parameters u = −1, b = 0.5, α = 0.5, β = 1, ω 0 = 1, and ∆ω = 0.2. In Figs. 1(c),(d), we demonstrate the frequency distribution. We define the phase of each oscillator as Im log Z i (x, y). We calculate and plot the time-averaged variation of the phase, setting the time window as T o = 10. In Figs. 1(c),(d), we only plot the frequency of the first component of oscillators at each site. We also plot the amplitude of the first component of oscillators |Z 1 (x, y)| in Fig. 1(e).
In Fig. 5, we show the dynamics of nonlinear oscillators in the same way as in Fig. 1. The parameters used are u = −1, u = 0.02, a = 2, b = 0.5, α = 0.5, β = 1, ω 0 = 1, and ∆ω = 0.2. We set the time step dt = 0.005. We plot the time-averaged variation of the phase, i.e., the frequency of each oscillator at the first component of the lattice point in Figs. 5(a),(b). We also plot the amplitude of the first component of oscillators |Z 1 (x, y)| in Fig. 5(c).
In Fig. 7, we also calculate the dynamics of nonlinear oscillators in the same way as in Fig. 1. The parameters used are u = −1, b = 0.5, α = 1, β = 1, ω 0 = 1, and ∆ω = 0.2. We set the time step dt = 0.005. We plot the time-averaged variation of the phase, i.e., the frequency of each oscillator at the first (second) component of the lattice point in Fig. 7(a) (Fig. 7(c)). We also plot the indicator of the local order parameter (12) of the first and fourth (second and third) components in Fig. 7(b) (Fig. 7(d)).

Calculations on the phase diagram
To obtain Fig. 2, we numerically calculate the dynamics of our first model (Eqs. (1), (4)) at different strengths of the linear coupling and inhomogeneity of the natural frequencies. We introduce the parameter of the strength of the linear coupling c by rewriting the equation of our model (1) as We change the parameter c from 0.2 to 2 at intervals of 0.2. We also change the inhomogeneity of the natural frequencies ∆ω from 0 to 2 at intervals of 0.1. We fix the other parameters as u = −1, b = 0.5, α = 0.5, β = 1, and ω 0 = 1. At each parameter, we calculate the dynamics and frequencies of oscillators in a 20 × 20 square lattice by using the fourthorder Runge-Kutta method as in Fig. 1. Then, at time t = 100, we calculate the standard deviations of the frequencies of the outermost oscillators and the bulk ones whose x and y coordinates satisfy 6 ≤ x ≤ 15 and 6 ≤ y ≤ 15. We plot the standard deviation of the frequencies of the outermost oscillators divided by that of the bulk ones in Fig. 2.
We also perform the numerical clustering of twodimensional data composed of the pair of the standard deviations of the frequencies of the edge and bulk oscillators. We utilize the k-means method to conduct such clustering. We set the number of clusters as k = 2. We determine the symbol used at each point in Fig. 2 based on which cluster the corresponding two-dimensional data point belongs to.

Numerical methods of Lyapunov exponents and Lyapunov vectors
We calculate the Lyapunov exponents and vectors of our first model (Eqs. (1), (4)) by utilizing the Shimada-Nagashima algorithm [60] and Ginelli's algorithm [61]. In those calculations, we arrange 10 × 10 sites and impose the open boundary condition. We deal with real and imaginary parts of the state variables separately as in the numerical calculation of our model (see Appendix A 1). We first calculate the dynamics of our first model from the random initial condition as in Fig. 1. Here, we set the time step as ∆t = 0.1. The obtained trajectory is used to determine the infinitesimal rate of change of the difference between the perturbed and original trajectories. We allow an initialization period of T ini = 1000 to stabilize the dynamics into a chaotic attractor. We repeatedly perform the QR decomposition based on the matrices derived from the linearization of the equation and the obtained trajectory from t = T ini to t = T ini + 20000. The averages of the logarithmic diagonal entries of the R matrices converge to Lyapunov exponents. To omit the effect of the initial condition, we use the R matrices obtained from the time t = T ini + 10000 to t = T ini + 20000 to calculate the Lyapunov exponents. We also utilize those R matrices to calculate the Lyapunov vectors. We iteratively multiply the inverse of obtained R matrices to a randomly-obtained upper-triangle matrix. Finally, the product of the obtained upper-triangle matrix and the Q matrix obtained in the calculation of the Lyapunov exponents represents the set of the Lyapunov vectors.
We also utilize the Shimada-Nagashima algorithm [60] and Ginelli's algorithm [61] to calculate the Lyapunov exponents and vectors of the other models (Eqs. (7), (10)). The parameters of models are set to be equal to those utilized in the calculation of dynamics (Figs. 6, 8). Furthermore, in the model utilizing the Hermitian linear coupling, we set the initialization period as T ini = 100000. The other parameters used in the calculations of Lyapunov exponents and vectors are the same as those in the calculations on our first model.

Numerical calculations on the emergent nonlinearity-induced boundary
To show the emergence of the extra topological boundary modes by the emergent nonlinearity-induced boundary, we diagonalize the state-dependent HamiltonianH(Z) under the open (periodic) boundary condition in the x (y) direction. We perform the inverse Fourier transformation of the Hamiltonian in Eq. (2) used in our first model only in the x direction. We divide the real and imaginary parts of the state variables as in the simulation of the dynamics of our model and consider the twice size of the matrix compared to the Hamiltonian in Eq. (2). We consider the 1 × 20 supercell structure. To determine the strength of the nonlinear loss terms (i.e., −β|z| 2 z), we first diagonalize the Hamiltonian without nonlinear loss terms. We obtain four edge modes v i j (x) (i = 1, · · · , 4 is the index of the edge modes and j = 1, · · · , 4 represents the component of oscillators at each site) with maximum imaginary parts of eigenvalues Im E 0 at the wavenumber k = 0. Then, we introduce the on-site loss −ia(|v 1 where the coefficient is set to be a = Im E 0 to balance gain and loss applied to the edge modes v i j (x). Thus, this on-site loss corresponds to the case that the edge modes in the Hamiltonian without the nonlinearity-induced on-site loss are fully amplified. Finally, we diagonalize the Hamiltonian with this nonlinear on-site loss term at each wavenumber and obtain the dispersion relation in Fig. 4(a) and the eigenvector of the edge mode in Fig. 4(b).
We also calculate the dynamics of the extra boundary modes. As in the calculation of the dynamics in Fig. 1, we use the fourth-order Runge-Kutta method and consider the 20 × 20 square lattice with the open boundary to simulate the dynamics of nonlinear oscillators with linear couplings described by the Hamiltonian in Eq. (4). By using the parameters α = 5 × 10 −5 and β = 1 × 10 −4 , we realize the situation that the nonlinearity has a negligible effect on the dynamics at the initial stage of the simulation. We consider the random initial condition, where real and imaginary parts of the initial state variables are the random value from uniform distributions ranging [−0.001, 0.001]. We set the time step dt = 0.005. The other parameters used are u = −1, b = 0.8, ω 0 = 0.1, and ∆ω = 0.2.

Numerical demonstrations of the on-demand pattern designing
To demonstrate the possibility of arranging synchronized oscillators in a desired pattern, we calculate the dynamics of our first model (Eqs. (1), (4)) under the existence of damped oscillators (see Fig. 10). We implement damped oscillators by setting α = −100 at the corresponding sites. We arrange 15 × 30 sites and impose the periodic boundary condition both in the x and y direction. As in the numerical simulation of the dynamics of our model, we utilize the fourth-order Runge-Kutta method. We start the calculation from the random initial condition, where real and imaginary parts of the initial state variables are the random value from uniform distributions ranging [−0.1, 0.1]. We set the time step dt = 0.005 and use the parameters u = −1, b = 0.3, α = 0.5, β = 1, ω 0 = 0.1, and ∆ω = 0.2. We calculate the frequency of each oscillator from the time-averaged variation of the phase, setting the time window as T o = 20. We also set the time window to calculate the indicators of local order parameters as T = 100.

Numerical demonstrations of the defect detection
In the demonstration of the defect detection by using the topological synchronized state (TSS) (cf. Fig. 11), we calculate the dynamics of our first model (Eqs. (1), (4)) under the periodic boundary condition by the fourth-order Runge-Kutta method. To implement the broken oscillators in Fig. 11(c), we change the parameters α at the corresponding sites from modes, there is a possibility that bulk oscillators stably synchronize in some models, which indicates the absence of the TSS in such models. Here we consider the following Hamiltonian of topological lasing modes: where σ x,y,z are the Pauli matrices, and I 2 is the 2 × 2 identity matrix. This Hamiltonian exhibits lasing modes localized only at the right side of the system. Figure 13 shows the band structure calculated under the open (periodic) boundary condition in the x (y) direction, which shows the existence of lasing edge modes. We arrange two Stuart-Landau oscillators (cf. Eq. (1)) at each site of a 20 × 20 square lattice and combine them by the linear coupling described by the Hamiltonian (B1). We assume the open (periodic) boundary condition in the x (y) direction. By calculating the dynamics of this model as in Fig. 1, we find that the bulk and left-side oscillators also synchronize. Figure 14 shows the frequency of the oscillator at each site. One can confirm that the bulk and left-side oscillators exhibit a homogeneous frequency that is different from that of the rightside oscillators. Thus, the bulk and left-side oscillators form another synchronized cluster, which indicates the absence of the TSS. We note that similar cluster synchronization is studied in one-and two-dimensional nonlinear topological systems [19].

Coexistence of synchronized edge oscillators and damped bulk ones
By tuning the parameters, we show that our model (Eqs. (1), (2)) can exhibit a different coexistence state of synchronization where the bulk oscillators are damped while the edge ones are synchronized. Here, we set the parameter α in the Stuart-Landau equation (cf. Eq. (1)) to be negative and larger than −|b|. The negativity of α implies the damping of selfoscillations of the bulk oscillators, while the lasing edge modes of linear couplings still lead to synchronized oscillations at the edge of the sample in this parameter region. Figure 15 shows the frequency and amplitude distributions obtained from the numerical calculations of our model at α = −0.2. We can confirm that the bulk oscillators exhibit almost zero amplitudes, while the edge ones synchronize. It is noteworthy that the bulk oscillators slightly oscillate at frequencies close to that of the edge ones due to their interaction. We define the TSS as the coexistence of the synchronized edge oscillators and the desynchronized bulk oscillators, while the synchronized state observed in this calculation exhibits bulk oscillators without self-oscillations and thus is not categorized in the TSS.

Appendix C: TSS without fluctuations of the natural frequencies
Here we show that the fluctuations of the natural frequencies are unnecessary to realize the TSS as in conventional chimera states [28,29]. We calculate the dynamics of the first model of the TSS (Eqs. (1), (2)) without fluctuations of natural frequencies, i.e., ∆ω = 0. Figure 16 presents the frequency distribution obtained under the homogeneous natural frequency. We can confirm that the edge oscillators vibrate at a homogeneous and constant frequency, while the bulk oscillators exhibit inhomogeneous and unstable frequencies. Therefore, the TSS appears even under the homogeneous natural frequency. We set the parameters as u = −1, b = 0.5, α = 0.5, β = 1, and ω 0 = 1 and consider the random initial condition as in Fig. 1.  (Fig. 1(c),(d)), the edge oscillators exhibit homogeneous and constant frequencies, while the bulk ones vibrate at space-and time-varying frequencies. Therefore, the edge oscillators synchronize, while the bulk ones desynchronize, which indicates the emergence of the TSS. The parameters used are u = −1, b = 0.5, α = 0.5, β = 1, ω0 = 1, and ∆ω = 0.

Stability of the TSS via a one-dimensional model
Here we discuss the stability of the TSS via the onedimensional chain of Stuart-Landau oscillators which effectively describes the edge oscillators in our first model (Eqs. (1), (2)). Since the edge oscillators in our model exhibit larger amplitudes than those of the bulk ones (cf. Fig. 1(e)), the interaction between the edge and bulk oscillators is smaller than that between the edge ones. Therefore, we can assume that the linear couplings between the edge and bulk oscillators only have a perturbative effect on the dynamics of the synchronized edge oscillators, and thus the model without the hopping terms in the x direction, H 1d jk (y, y ) = uδ y,y + δ y+1,y + δ y−1,y 2 ( can describe the edge oscillators in our first model, where δ y,y is the Kronecker delta and σ i and I 2 are the ith component of the Pauli matrices and the 2 × 2 identity matrix respectively. We introduce ξ j (y) to describe the perturbative effect from the bulk oscillators, while we show that the TSS is robust against such perturbation and inhomogeneity of parameters. We note that the chaotic dynamics in the bulk vanishes the long-term correlation of interaction between the edge and bulk oscillators, and thus ξ j (y) can be considered as noise.

Re E wavenumber ky
Im E wavenumber ky First, we obtain the stationary synchronized solution of Eq. (D1) under the absence of the interaction from bulk oscillators and the homogeneous natural frequencies. We conduct the Fourier transformation of the Hamiltonian used in the onedimensional model (D2) and obtain H 1d (k) = (u + cos k y ) I 2 ⊗ σ z + sin k y I 2 ⊗ σ y + ibσ x ⊗ σ x .
(D3) This Hamiltonian exhibits the maximum imaginary part of the eigenenergies Im E = √ b 2 − u 2 at k y = 0 and the corresponding eigenvector (ib, By utilizing this eigenvector, we construct the following homogeneous stationary (periodic) solution: where ω 0 is the homogeneous natural frequency. Since the gapped bulk bands require |u | < |b|, the absolute values of the components of this stationary solution are the same, which leads to the same intensity of the nonlinear terms. The stationary solution (D4) represents the frequency-synchronization of the edge oscillators, where they oscillate at the same frequency ω 0 .
To confirm the stability of the frequency-synchronization of the edge oscillators, we conduct linear stability analysis around the stationary solution (D4). We linearize the Stuart-Landau equationŻ = (iω + α − βZ 2 )Z and obtain d dt where Z = X + iY (X, Y are real) and Z 0 = X 0 + iY 0 is the state around which the equation is linearized. The eigenvalues of the coefficient matrix in this equation are iω +α−X 2 0 −Y 2 0 , iω + α − 3X 2 0 − 3Y 2 0 and have negative real parts when we linearize the equation around the stationary solution (D4). Therefore, these terms attenuate the fluctuations from the stationary solution, which leads to the linear stability of the onedimensional chain of nonlinear oscillators. We numerically diagonalize the linearized equation of Eq. (D1) and obtain the dispersion relation in Fig. 17. We confirm that all the eigenvalues have nonpositive imaginary parts, which indicates the linear stability of the stationary solution. Therefore, the synchronized state in the one-dimensional model (D1) is robust against the perturbations, such as the inhomogeneity of the natural frequencies and the perturbative interactions from bulk ones. Since the one-dimensional chain of oscillators and its stationary solution describes the synchronized edge oscillators in our first model (Eqs. (1), (2)), the TSS is also stable against the existence of disorders.

Robustness of the electrical circuit realization of the TSS
TSS is robust against disorder due to its topological nature and thus is realizable without a fine-tuning of parameters. Especially, while there should be inhomogeneity in circuit constants of elements constructing the TSS circuit (cf. Fig. 12), TSS remains stable against such disorders. To confirm the stability of TSS against the inhomogeneity of circuit constants, we numerically calculate the dynamics of our first model (Eqs. (1), (2)) under the fluctuating coupling strengths. We introduce the fluctuations of parameters corresponding to the inhomogeneous circuit elements by setting the proportion of deviation in the coupling strength from a site (x; a) to a site (x ; b) to be the same as that from a site (x ; b) to a site (x; a). We also consider the fluctuation of the parameters α and β in Eq. (1), which is determined from the characteristics of the nonlinear current resources in the circuit of TSS (cf. Fig. 12). We set the maximum width of these fluctuations as 10 percent of the mean values, and determine them as the distributions follow the uniform distribution. Figure 18 shows the frequency distribution obtained from the simulation. One can confirm the emergence of the TSS, that is, the edge oscillators synchronize while the bulk ones exhibit chaotic behavior. Therefore, the TSS is robust against the fluctuations of parameters in Eq. (1), which should be of practical advantage to experimentally realize the TSS.

Appendix E: Dispersion relation of the effective Hamiltonian
Here we calculate the dispersion relation of the effective Hamiltonian derived from the linearization of the equation of our first model (Eqs. (1), (2)). Such effective Hamiltonian includes the on-site loss induced by the nonlinear terms, which depends on the state around which the linearization is performed. Specifically, such on-site loss terms are obtained from the linearization of the Stuart-Landau oscillatorṡ Z = (iω + α − βZ 2 )Z as d dt where Z = X + iY (X, Y are real) and Z 0 = X 0 + iY 0 is the state around which the equation is linearized. We note that the coefficient matrix in this equation has eigenvalues Thus, if the square of the amplitude is X 2 0 + Y 2 0 = α/β (equal to the amplitude of the isolated Stuart-Landau oscillator), the obtained on-site term leads to attenuation of the fluctuation.
To determine the intensities of those on-site loss terms, we first simulate the dynamics of our first model (Eqs. (1), (2)) in a 20 × 20 square lattice under the periodic (open) boundary condition in the x (y) direction. Here we consider the random initial condition as in Fig. 1. We obtain the state variables of the oscillators in the 10th row at t = 200 and linearize the equation around them. In the calculation of the dispersion relation, we consider a 1 × 20 supercell structure. Finally, we numerically diagonalize the obtained effective Hamiltonian and plot the dispersion relation as shown in Fig. 19. We find bulk modes exhibiting positive imaginary parts of eigenvalues, which lead to the chaotic dynamics of the bulk oscillators. There are also gapless modes with negative imaginary parts of eigenvalues that are localized at the edge of the sample. These edge modes should correspond to the edge-localized Lyapunov vectors in our first model (cf. Figs. 3, 23). We note that the non-Hermiticity of the effective Hamiltonian plays an important role, because nonzero imaginary parts of eigenvalues correspond to nonzero Lyapunov exponents and are unique to non-Hermitian Hamiltonians. It is also noteworthy that the bulk bands with the positive imaginary parts of the eigenvalues are derived from the repulsion of the bulk bands of the original linear couplings with zero imaginary parts (cf. Figs. 24, 25). To realize such degenerated imaginary bulk bands and lasing edge modes, (at least) four oscillators seem to be necessary at each site.
While nonlinearity generates random on-site loss terms in the effective Hamiltonian, topological modes are robust against such disorder. Especially, previous research [43] shows that the edge modes of the Hamiltonian describing the linear coupling of our first model (2) are robust against perturbative on-site gain and loss. In contrast to conventional topological insulators, the bulk bands are also gapless in Fig. 21. However, the edge and bulk modes are not degenerate, and thus edge modes remain in the effective Hamiltonian. We note that the state-dependent Hamiltonian analyzed in Fig. 4(a),(b) also includes random on-site loss terms, while they have no effects on the robust existence of the topological modes in the state-dependent Hamiltonian.
We note that the imaginary parts of the eigenvalues of the effective Hamiltonian represent the short-term stability or instability of the system, while the Lyapunov exponents represent the long-term stability or instability. The Lyapunov exponents and vectors correspond to the time-averaged behavior of the effective Hamiltonian and thus are different from the eigenvalues and eigenvectors of the effective Hamiltonian obtained at each period. However, since the edge modes are robust against the unsteady disordered effect of nonlinear terms, the edge-localized modes and their dissipative behaviors remain even after the time averaging.   2)). The relative indices correspond to those in Fig. 3. We confirm the large IPRs at small relative indices (smaller than about 0.05) and large relative indices (larger than about 0.95). The parameters used are u = −1, b = 0.5, α = 0.5, β = 1, ω0 = 0.2, and ∆ω = 0.2.

Amplitude chaos in the TSS
Amplitude oscillators such as Stuart-Landau oscillators can exhibit two types of chaos, amplitude chaos and phase chaos. Amplitude chaos [48] accompanies phase slips where the amplitudes of oscillators go to zero, and thus the phases jump. Such amplitude-dependent behavior is unique to amplitude oscillators and thus cannot be described by their phase equations.
We confirm the amplitude chaos in our first model (cf. Eqs. (1), (2)). We numerically calculate the dynamics of our model as in Fig. 1. Figure 20 plots the amplitudes of the first component at each site and emphasizes the sites with small amplitudes. A part of bulk oscillators exhibit almost zero amplitudes, and thus phase slips can occur. Therefore, amplitude chaos appears in the bulk of our model. We note that our model can also exhibit the phase chaotic behavior (i.e. chaos irrelevant to the phase slips observed in the numerical calculation).

Inverse participation ratios of the Lyapunov vectors in the wavenumber space
While in the main text we analyze the inverse participation ratios (IPRs) of the Lyapunov vectors in the real space, here we calculate and discuss the IPRs of our first model (cf. Eqs. (1), (2)) in the wavenumber space. Especially, we find that the IPRs in the wavenumber space classify the bulk-extended Lyapunov vectors into two groups. To obtain the IPRs in the wavenumber space, we conduct the Fourier transformation of the Lyapunov vectors of our model, and calculate the following value: whereṽ i (k) is the Fourier component of the ith component of the Lyapunov vectors corresponding to the wavenumber k. Figure 21 presents the IPRs of the Lyapunov vectors in the wavenumber space. We find that the first and last 10 percents of the Lyapunov vectors exhibit large IPRs in the wavenumber space, while the others exhibit small IPRs. We note that the Lyapunov vectors corresponding to the small Lyapunov exponents also exhibit large IPRs in the real space. The Lyapunov vectors exhibiting large Lyapunov exponents and large IPRs in the wavenumber space should belong to a different class from the other bulk-extended ones. The large IPRs in the wavenumber space are related to the band structure of the effective Hamiltonian obtained from the linearization of the equations of motion (see Appendix E and Fig. 19). The bulk modes with large positive imaginary parts of eigenenergies are concentrated around the wavenumber k y = 0. Therefore, we can describe the Lyapunov vectors corresponding to the large Lyapunov exponents as superpositions of the bulk modes only around k y = 0, which leads to the large IPRs in the wavenumber space. We can discuss the large IPRs of some edge-localized Lyapunov vectors and the small IPRs of the other ones in a similar way.
Previous research of Hamiltonian chaos [62] has studied the participation ratio of coefficients of the eigenstates in a quantum chaotic system when they are described as the superposition of the eigenstates of decoupled harmonic oscillators. Such participation ratios classify the eigenstates in chaotic systems into regular and ergodic ones. In our model, the Lyapunov vectors exhibiting large (small) IPRs in the wavenumber space corresponds to regular (ergodic) ones. However, we leave the full understanding of the topological synchronized state in terms of Hamiltonian chaos to a future work.

Appendix G: Real-space distribution of frequencies and
Lyapunov vectors in a 10 × 10 square lattice Related to the calculation in Fig. 3, we calculate the dynamics of the first model (cf. Eqs. (1), (2)) of the TSS in a 10 × 10 square lattice. Figure 22 shows snapshots of the frequencies of oscillators. We can confirm the frequency-synchronization of the edge oscillators and the desynchronization of the bulk oscillators. We thus expect that the TSS can appear independently of the system size.
We also calculate the Lyapunov vectors [26,61] in the 10 × 10 lattice model and find that some of them are localized to the edge of the system (cf. Fig. 3(b),(c)). Here we plot the Lyapunov vectors in the real space and directly demonstrate such localization of the Lyapunov vectors. Figure 23 shows some examples of Lyapunov vectors. One can see that the Lyapunov vector associated with a positive Lyapunov exponent (corresponding to Fig. 23(a)) has large amplitudes at a large number of bulk sites. Meanwhile, most of the bulk sites have small amplitudes in the Lyapunov vector of the relative index around 0.875 (corresponding to Fig. 23(b)), which indicates its localization to the edge. Furthermore, the Lyapunov vector of the largest index is strongly localized to corner sites as shown in Fig. 23(c), which leads to the large value of its inverse participation ratio (see Eq. (6) and Fig. 3(c)). Thus, these observations of the Lyapunov vectors in the real space are consistent with the results in Fig. 3. linear couplings is H(k) = aH QWZ + iu I 2 ibI 2 ibI 2 H QWZ − iu I 2 (H1) = (u + cos k x + cos k y ) × a + 1 2 (I 2 ⊗ σ z ) + a − 1 2 (σ z ⊗ σ z ) + sin k y a + 1 2 (I 2 ⊗ σ y ) + a − 1 2 (σ z ⊗ σ y ) in the wavenumber space, where H QWZ is the Hamiltonian of the Qi-Wu-Zhang model [45] and I 2 , σ x,y,z are the 2 × 2 identity matrix and the Pauli matrices respectively. This Hamiltonian exhibits nonzero Chern numbers and lasing edge modes, i.e., gapless localized modes with the positive imaginary parts of the eigenvalues. We check the existence of such edge modes by calculating the dispersion relation under the open (periodic) boundary condition in the x (y) direction. Figure 25 shows the dispersion relation of the lasing edge modes. One can find the gapless modes with positive imaginary parts of eigenvalues larger than those of the bulk modes. As discussed in the main text, these edge modes lead to the TSS in the linearly coupled Stuart-Landau oscillators (cf. Fig. 5).

Appendix I: Absence of synchronized pattern in a topologically trivial system
The synchronization pattern discussed in Fig. 10 requires the topology of linear coupling. To clearly demonstrate this, we simulate the dynamics of our first model (Eqs. (1), (2)) in a topologically trivial parameter regime where the linear couplings exhibit no edge modes. We set the parameter u = −3 (cf. Eq. (2)) and the other parameters to be equal to those used in Fig. 10. Figure 26 shows the frequencies and the indicators of order parameters in the nontopological parameter region. We can confirm that the frequencies of the oscillators around damped ones are close to those away from damped regions and the indicators of order parameters are small, which implies the absence of TSS. Therefore, the topology and edge modes of the linear couplings play a vital role in the synchronization pattern designing.

Appendix J: Circuit equation for realizing the TSS
Previous studies [46,54,[63][64][65][66] have proposed and realized topological edge modes in electrical circuits. The relation between the input current I a and the electrical potential V a at each node is described by Kirchhoff's law where C ab (C a ) is the inverse of the impedance between site a and b (site a and ground). iJ ab plays the role of the effective Hamiltonian in electrical circuits. Resistors play the role of imaginary couplings, while capacitors and inductors lead to real terms in the effective Hamiltonian. By tuning the impedances of circuit elements, one can realize the effective Hamiltonian imitating topological materials. If we introduce capacitors with capacitance C as current resources, we can which governs the time evolution of the voltage at each node. To realize the electrical circuit of topological lasing modes, we utilize negative impedance converters with current inversion (ICINs), which is discussed in a previous study [54]. The ICINs realize nonreciprocal couplings (C ab = C ba ) that cannot be realized by using resistors, inductors, and capacitors. Input current I in and output current I out are described as where V in (V out ) is the electric potentials at the input (output) side. We assume R = R in the electrical circuit of the TSS (cf. Fig. 12). Then by tuning the parameters of circuit elements, one can construct the effective Hamiltonian iJ that is equal to the Hamiltonian of the topological insulator laser [42,43] used in our first model (2). It is noteworthy that impedances of capacitors and inductors and thus the effective Hamiltonian iJ depends on the driving voltage frequency.
To directly realize an electrical circuit whose voltage dynamics is described by nonlinear equations similar to our model, we next omit capacitors and inductors and construct a frequencyindependent circuit of topological lasing modes. To effectively realize imaginary terms without capacitors and inductors, we prepare twice the number of nodes compared to the circuit in Fig. 12 and express one complex-valued state variable by using the voltages at two nodes as Z = V 1 + iV 2 . Then, real and imaginary couplings are substituted as where C is real. These real couplings are realized by the substitution of circuit elements shown in Fig. 12(b). Each node in the topological circuit of the TSS consists of van der Pol circuits [52,53], which realize nonlinear oscillators. To construct the van der Pol circuits, we utilize nonlinear resistors whose conductance is described asC(V i ). The dynamics of the van der Pol circuits is derived by Kirchhoff's law as where C is the capacitance of the capacitor and I in is the input current. When we can expandC(V i ) asC(V i ) = α − βV 2 i + O(V 3 i ), the dynamics of the van der Pol circuits can be approximated by the Stuart-Landau oscillators in the leading order. Therefore, the entire dynamics of the circuit imitates the first model of the TSS (Eqs. (1), (2)).