Topological synchronization of quantum van der Pol oscillators

To observe synchronization in a large network of classical or quantum systems demands both excellent control of the interactions between the nodes and very accurate preparation of the initial conditions due to the involved nonlinearities and dissipation. This limits the applicability of this phenomenon for future devices. Here, we demonstrate a route towards significantly enhancing the robustness of synchronized behavior in open nonlinear systems that utilizes the power of topology. In a lattice of quantum van der Pol oscillators with topologically motivated couplings, boundary synchronization emerges in the classical mean field as well as the quantum model. In addition to its robustness against disorder and initial state perturbations, the observed dynamics is independent of the underlying topological insulator model provided the existence of zero-energy modes. Our work extends the notion of topology to the general nonlinear dynamics and open quantum system realm with applications to networks where specific nodes need special protection like power grids or quantum networks.

For many quantum mechanical applications dissipation is often regarded as an undesirable yet unavoidable consequence because it potentially degrades quantum coherences and renders the system classical. However, interactions with the environment can also be considered a fundamental resource for striking collective effects typically impossible in Hamiltonian systems [1][2][3]. A hallmark of such collective behavior in nonequilibrium systems is the phenomenon of synchronization: in the complete absence of any time-dependent forcing from the outside, a group of oscillators adjusts their frequencies such that they spontaneously oscillate in unison [4,5]. Synchronization is intimately related to the phenomenon of self-sustained oscillations, where a system maintains a periodic motion in an autonomous fashion [6]. In this sense, interactions with the environment and nonlinearities in the equations of motion represent prerequisites for synchronization. With the recent developments in quantum technology which allow one to exquisitely tailor both the system and environmental properties, synchronization has emerged in the quantum domain with various different examples ranging from nonlinear oscillators to spin-1 systems, superconducting qubits, optomechanics and ensembles of atoms that demonstrate synchronization in the quantum regime [7][8][9][10][11][12][13][14][15][16][17][18][19].
In both quantum and classical systems, unavoidable imperfections -local deformations caused by ambient conditions as well as long-term degradation -have a significant impact on the collective behavior and can even destroy synchronicity altogether. Moreover, collective synchronized dynamics in large networks often sensitively depends on fine-tuned initial conditions. It is therefore desirable for a reliable performance of future devices to identify universal principles to enhance the robustness of synchronization. In this work, we demonstrate that the * cwwaechtler@berkeley.edu power of topology can be exploited for this task. Topological insulators describe a special class of solids exhibiting an insulating bulk but symmetry protected conducting surface states, known as topological edge states [20][21][22][23][24]. These edge states display a surprising immunity to a wide range of local deformations, inherently avoiding backscattering over broad energy ranges and even circumventing localization in the presence of disorder [25].
While the robustness of topological systems is well established for isolated systems, it is a nontrivial problem to determine whether this feature is still available in open setups. Recently, the study of out-of-equilibrium systems has moved to the forefront of research in topological band structures, motivated largely by the desire to create interesting states of matter with properties beyond those achievable in equilibrium [26][27][28][29][30][31][32]. While topological band and bandgap structures are inherently tied to linear systems, we currently observe rapidly growing interest to generalize topological concepts to nonlinear systems, which in turn has sparked the field of nonlinear topological photonics and the phenomenon of topological lasing [33][34][35][36][37].
Despite these recent advancements, applying topological concepts to synchronization phenomena is still largely unexplored, with only a few studies on classical nonlinear oscillators [38][39][40], yet with already some remarkable results. For example, Sone et al. [40] showed that linearly coupled Stuart-Landau oscillators arranged on a square lattice may exhibit chaotic bulk motion and robust synchronized edge oscillations when the coupling Hamiltonian promotes topological lasing modes. The authors termed this coexistence state 'topological synchronized state', which may be observed for Hermitian as well as non-Hermitian coupling Hamiltonians. Furthermore, the involved nonlinearity may also induce emerging effective boundaries, which lead to extra topological modes unique to nonlinear systems. However, to the best of our knowledge, similar studies for quantum systems are entirely missing so far. In this work, we systematically investigate the potential of topological insulator models to promote topological synchronization in classical and quantum systems. In order to provide such a comprehensive study, we choose the van der Pol (vdP) oscillator [4,5] as our building block to construct topological lattices. Since the quantum version of the vdP oscillator [7,8,10] reduces to its classical analog at the mean field level (specifically to the Stuart-Landau oscillator, which is the weakly nonlinear limit of the classical vdP oscillator), we are able to study the influence of topological band structures in both regimes. We show not only that boundary synchronization emerges in the classical as well as the quantum regime if the underlying lattice possesses nontrivial topology, but also that it is protected over a wide range of local disorder. Our results pave the way towards utilizing topology as a powerful tool for enhancing the robustness of collective dynamics in nonlinear open quantum systems.

A. Topological van der Pol oscillator network
We study a lattice of N sites, each consisting of a harmonic oscillator labeled by the index j = 1, ..., N and we assume identical frequencies ω 0 of each oscillator. The Hamiltonian of the system can be written as a tightbinding Hamiltonian (1) where a † j (a j ) denote creation (annihilation) operators of bosonic particles at lattice site j. This general form of the system Hamiltonian allows us to realize different topological lattices simply by modifying the couplings λ jj ′ . Recently, it has been realized that a quantum harmonic oscillator subject to one-phonon gain with rate κ 1 and two-phonon loss with rate κ 2 represents the quantum analogue of the classical vdP oscillator [7][8][9][10][11]. In an open quantum system approach and using the nota- , the dynamics of the system density matrix (t) describing a network of vdP oscillators is then given by the master equatioṅ Throughout this work we focus on the weak dissipation regime, that is we assume that κ 1 , κ 2 ≪ ω 0 . For realistic systems, this regime corresponds to long coherence times and weak coupling to the reservoirs such that a description of the system dynamics in terms of a Lindblad master equation [cf. Eq. (2)] is valid. Furthermore, we are interested in the regime where κ 2 > κ 1 as in this case each oscillator remains close to its ground state [7].

A. Mean field approximation and linear stability analysis
From the full quantum model defined in Eq.
(2), one may derive classical equations of motion for the expectation values α j = ⟨a j ⟩ ∈ C by performing a mean field approximation. The governing equation of the complexvalued mean field amplitudes α = (α 1 , . . . , α N ) is then given byα where H = H 0 + H top denotes the matrix corresponding to the Hamiltonian H [cf. Eq. (1)] with diagonal matrix H 0 = ω 0 1 and the matrix describing the coupling between oscillators H top . Furthermore, ⊙ denotes the Hadamard product defined as (α ⊙ α) n = α n ⋅ α n .
In the absence of any coupling between lattice sites (λ jj ′ ≡ 0) each oscillator will approach its respective periodic steady stateᾱ j (t) =Ā exp(−iω 0 t + ϕ j ) with A = κ 1 2κ 2 and arbitrary phases ϕ j . However, the interactions enable the emergence of collective synchronized motion, whose existence is intimately related to the stability of the fixed point α FP defined viaα FP = 0. From Eq. (3), we find that the fixed point is given by α FP = 0. The stability of α FP may be analyzed in terms of a corresponding model linearized around α FP : where we have omitted second and higher order terms. The Jacobian matrix J with entries J jj ′ = ∂α j ∂α j ′ is evaluated at α FP resulting in the matrix In particular, Eq. (4) only contains the linear terms of Eq. (3). In general, the fixed point α FP is stable as long as the real part of all eigenvalues of J are negative. In contrast, if at least one eigenvalue has positive real part, the fixed point is dynamically unstable and the corresponding eigenstate will exponentially grow. However, the nonlinear damping in Eq. (3) counteracts this exponential amplitude increase, leading eventually to stable oscillations.
Since H top is a Hermitian matrix, it has real eigenvalues µ (l) , where l = 1, . . . , N is an index for the different eigenvalues. Moreover, the diagonal entries of J are given by −iω 0 + κ 1 2, such that ν (l) = −i(ω 0 + µ (l) ) + κ 1 2 are the eigenvalues of J. Hence, for κ 1 > 0 (which we assume throughout this work), all eigenstates are linearly unstable. Moreover, for κ 2 > κ 1 we expect small oscillation amplitudes which remain close to the fixed point even if the coupling λ jj ′ is comparable to the intrinsic frequency ω 0 . Consequently, the dynamics of the oscillator lattice (after relaxation) may be well approximated by a superposition of eigenvectors, where µ (l) and α (l) are the eigenvalues and their respective eigenvectors of H top , and c (l) denotes scalar coefficients. This may suggest that the system dynamics is fully described in terms of the linear model and thus by Eq. (4). However, our results -discussed in Sec. III -indicate that such a simple picture is not sufficient to capture the full dynamical evolution of the system in the presence of dissipation and nonlinearities.

B. Quantum fluctuations
Solving the full open quantum system given by Eqs. (1) and (2) is a non-trivial task due to the large number of interacting oscillators and the involved nonlinearities in the dissipators. In order to gain further insight, we follow a similar derivation to the ones found in Refs. [11,41]: We move to a displaced frame via the displacement operator and define the density matrix in the displaced frame as . Interestingly, the terms linear in the operators a j and a † j vanish as long as Eq. (3) is fulfilled, such that we are left with terms quadratic and cubic in the operators; see App. A for details. By neglecting terms of order O(a 3 j ), we obtain an effective master equation of Lindblad form, Note, that the full dynamics of the system is now determined by Eq. (8) in combination with Eq. (3) because the time-dependent mean field amplitudes α j (t) appear in the effective Hamiltonian (9) and the time-dependent dissipation rates. Hence, in general both equations have to be solved simultaneously. However, as Eq. (3) is independent of the dynamics of the density matrix α (t), the mean field dynamics may be solved first and treated as time-dependent inputs for Eq. (8). In App. C, we compare this effective model to the full quantum model described by Eq. (2) for a single and two coupled quantum vdP oscillators in terms of Wigner functions.
As the Hamiltonian (9) is quadratic, the dynamics of the effective quantum model is fully described by the covariance matrix C mn (t) = Tr[ α (t){X m , X n } 2] with the quadratures X 2j−1 = (a j + a † j ) The equation of motion of the covariance matrix C(t) is given by [42] where B(t) and D(t) are determined from Eq. (8); see App. B. Throughout this work, we consider a pure coherent state as the initial state (0) = ⊗ j α j (0)⟩ ⟨α j (0) . In the comoving frame, such an initial condition corresponds to α (0) = ⊗ j 0 j ⟩ ⟨0 j which results in an initial covariance matrix of diagonal form, i.e., C(0) = 1 2 ⋅ 1, which reflects the Heisenberg uncertainty principle. Note, that the time-dependency of B(t) and D(t) arises from the time-dependent mean field amplitudes α j (t) (see also see App. B), such that Eqs. (3) and (10) are also solved simultaneously as mentioned previously.

C. Synchronization in classical and quantum systems
Synchronization of classical oscillatory systems can take many different forms. A primary example often discussed is the case of two strongly nonlinear systems which oscillate with different natural frequencies, and that adjust their rhythm due to a mutual weak coupling referred to as frequency synchronization. Additionally, the two oscillators may also synchronize their phases such they oscillate in-phase (or anti-phase). In the case of weak coupling, one can show that only phases rather than amplitudes of oscillators are relevant (see e.g. [43]). Conversely, a relatively strong coupling between two oscillators can affect not only their phases but also their amplitudes, and the features of synchronization in the presence of strong interactions are nonuniversal [4]. For example, on a finite lattice even if all oscillators are identical and they only interact locally with their nearest neighbors, synchronicity across the lattice is not guaranteed; examples will be discussed later in Sec. III.
In this work, we are concerned to what extent topological lattices and, in particular, edge states affect the synchronicity of self-oscillatory systems. As the lattices studied here host edge states that are zero dimensional, the corresponding oscillators are located at the boundaries and only interact via many bulk oscillators. If these bulk oscillators are neither frequency nor phase synchronized, one may not expect that the oscillators are able to synchronize their phases. However, they might still be able to oscillate with a common frequency. Thus, throughout this work we will refer to two classical systems j and j ′ being synchronized if the condition is fulfilled. Note, that if this conditions is fulfilled for all j and j ′ , the whole lattice oscillates with a single common frequency which is expected to be an eigenfrequency of the lattice (i.e. ω 0 + µ (l) ; see Sec. II A) because of the weakly nonlinear regime. Generalizing the previously introduced notion of synchronization to the quantum regime is challenging as phase space trajectories become ill defined concepts. To overcome this challenge, synchronization measures in terms of the Husimi-Q or Wigner phase space distributions [7,13,44], explicit limit cycles of system observables [45,46], or information-theoretical measures [47,48] have been proposed. In this work, we use a measure to quantify synchronization of two quantum systems j and j ′ that is based on their dimensionless quadratures as [49] The upper bound of the quantum synchronization measure arises from the Heisenberg uncertainty principle. The synchronization measure S c (j, j ′ ), also referred to as 'complete synchronization measure' [49], extends the classical concept of 'error' into the quantum domain, i.e., the smaller the variance of the quadrature differences is, the larger is the synchronization measure S c (j, j ′ ). It thus indicates how equivalent the dynamics of two quantum systems are, which in turn connects to the classical notion of two systems oscillating with the same frequency in-phase or anti-phase. Here, S c (j, j ′ ) would take on its maximum value; see also App. C.
As S c (j, j ′ ) may exhibit oscillations even in the long time limit, we perform time averages of this quantity, i.e., with initial time t i and final time t f . Throughout this work we are interested in the (periodic) steady state dynamics and discard transient dynamics, such that both t i and t f are taken to be larger than the transient relaxation time (see Sec. III A).

A. Preliminary remarks
Before we turn to the detailed analysis of two specific topological lattices in the next two sections, we would like to provide some general notions that will apply to both examples in order to minimize redundancies. Within each section, α (l) indicates an eigenstate with eigenvalue µ (l) of the respective topological (coupling) Hamiltonian, i.e., H SSH in Example (I) and H Kag in Example (II). Furthermore, the mean field amplitude of lattice site j is denoted by A j (t) and defined via A j (t) = [α j (t)+α * j (t)] 2 , which is proportional to the averaged position of each oscillator. Throughout, when random initial conditions are considered, we mean randomly distributed complex amplitudes and phases, i.e., a j (t = 0) = α r j exp (iϕ r j ) where α r j and ϕ r j are random variables chosen from the uniform distributions α r j ∼ U(0, 0.5) and ϕ r j ∼ U(0, 2π), respectively. Similarly, when random disorder with strength r is applied to the coupling between lattice sites j and j ′ we refer to λ jj ′ → λ jj ′ + δλ r jj ′ with δλ r jj ′ ∼ U(−r, r). Lastly, all numerical results shown are obtained after a significant relaxation time ω 0 t rel = 2 ⋅ 10 4 in order to discard transient effects.

The SSH model
In the following we investigate the interplay of topology and synchronization in the paradigmatic SSH model [50,51], a one-dimensional dimerized lattice with staggered nearest-neighbor hopping and time reversal, particle-hole and chiral symmetry (BDI class). Recently, implementations of this model under nonequilibrium conditions or in strongly nonlinear systems have been proposed [52][53][54][55][56] and experimentally realized [57][58][59][60][61][62]. The Hamiltonian of the model is given by where λ j = λ 1 if j is odd and λ j = λ 2 otherwise. The spectrum of the Hamiltonian (14) is symmetric with respect to the zero-energy axis (E = 0) as a result of the chiral symmetry present in the system. A phase transition of topological nature separates the trivial phase (λ 1 > λ 2 ) from the topological phase (λ 1 < λ 2 ). In Fig. 1 In addition, two degenerate zero-energy states occur within the band gap in the topological phase (δλ > 0), which are highly localized at the boundaries of the system [ Fig. 1(a)]. The edge states are topologically protected by the chiral symmetry and are thus robust against (symmetry-preserving) local perturbations as we discuss in more detail below.

Topological effects on synchronization at the mean field level
To incorporate the topological character of the SSH model into the network of vdP oscillators, the coupling matrix H top in Eq. (3) is chosen to be H SSH . In the following we compare the numerical results to our previous We start our discussion in the trivial phase (δλ < 0). In Figs. 1(b) and (c), we observe complete synchronization, i.e. all lattice sites j oscillate with the same frequency ω 0 + µ (l) . This frequency is determined by the choice of initial eigenstate α (l) with respective eigenvalue µ (l) . Hence, the oscillation frequency in panel (b) [(c)] is larger (smaller) compared to the intrinsic frequency ω 0 of the uncoupled vdP oscillators as µ (l) > 0 (< 0). Furthermore, the amplitude of the initial state α (l) j directly translates to the oscillation amplitudes A j (t) and their phases. Thus in accordance to our previous analytic considerations, the dynamics in the trivial phase is given by For δλ > 0, however, nontrivial topological effects on the dynamics of the lattice emerge as shown in Figs. 1(d) and (e): Even though the initial bulk eigenstate (shown in the left panels) exhibits only small amplitudes at the edges (j = 1 and j = 20) their oscillation amplitudes (highlighted in blue/pink) are comparable to the largest amplitudes in the whole chain. Moreover, they oscillate with the intrinsic frequency ω 0 (the vertical dashed lines are located at integer multiples of ω 0 t = 4πn, n ∈ N). In comparison to the trivial phase, where small initial amplitudes remain small for all lattice site, the edge modes in the topological phase are always excited. Note that the edge states are not strictly localized at the two boundary lattice sites, but rather extend exponentially into the bulk. Consequently, the dynamics of bulk oscillators close to the edge, e.g. j = 3 or j = 18, is a superposition of the initial state and the edge state with different frequencies. Thus, the dynamics of these oscillators is given by j and α edge j represent the initial and edge eigenstate, respectively.
Since it might be difficult in an experimental setup to initiate the system in a specific eigenstate α (l) of H SSH , the question arises how a randomly distributed initial state affects the dynamics and whether the previously discussed zero-energy edge mode synchronization persists in such a scenario. To this end, we show in the left panels of Figs. 2(a)-(d) the amplitude dynamics A j (t) for such random initial states. In the trivial phase shown in panels (a) and (b) there is no clear pattern of synchronization. We would like to note that some oscillators in Figs. 2(a) and (b) might be synchronized for short times, however, their phase difference is not constant over longer times.
In order to confirm our intuition that there are many frequencies participating in the dynamics of the vdP network, we perform a discrete Fourier transform and show the frequency spectrum in the right panels of Figs. 2(a) and (b). As expected there are multiple peaks centered around ω 0 and separated by a gap. These peaks correspond to the eigenfrequencies ω 0 + µ (l) of the system.
In the topological phase shown in Figs. 2(c) and (d) the oscillator dynamics (left panels) is similar to the one described previously for the trivial phase. However, there is an important difference: The edges (highlighted in blue/pink) do not exhibit frequency mixing but oscillate with one specific frequency, the intrinsic frequency ω 0 . Hence, the edges are synchronized with one another, yet due to the random initial conditions with a phase shift difference (e.g. different amplitudes at the dashed lines for j = 1 and j = 20). However, this phase difference remains constant over time. In the frequency spectrum of the discrete Fourier analysis, we observe now an additional sharp peak at ω = ω 0 confirming the existence of the synchronized edge modes.
In terms of the discrete Fourier analysis we are able to analyze the full spectrum of the coupling matrix H SSH for random initial conditions. In Fig. 2(e) we show the result of the frequency spectrum averaged over 10 realizations of random initial conditions as a function of δλ. Remarkably, the spectrum is identical to the eigenspectrum of H SSH , which confirms that the synchronized edge modes are only present in the topological phase (δλ > 0).
One of the hallmarks of topological insulators is that they exhibit extremely robust surface states since no local perturbation can change their global topology. For the SSH model the robustness arises from an underlying chiral symmetry. To test whether this extraordinary feature is also present in our open nonlinear system of vdP oscillators, we apply random disorder to the couplings between neighboring sites [cf. Eq. (14)]. In Fig. 3 we show the frequency spectrum of the discrete Fourier analysis for the coupling δλ = 0.8t 0 as a function of the disorder strength r for 10 realizations. While the frequen- cies within the upper and lower band spread over a wide range as the disorder strength is increased, the edge mode persists even for disorders as large as r = 1.5λ 0 before the bands start overlapping with the zero energy mode. The edge state quantum synchronization is topologically protected and robust even for large amounts of disorder. Parameters:

Quantum signatures of topological synchronization
In the previous section the focus has been the effects of topology on the mean field dynamics. Thus, the observed signatures are classical in nature even though the underlying model is quantum. By contrast, in this section we analyze the quantum fluctuations about the mean field amplitudes and investigate whether a similar interplay of topology and synchronization carries on beyond the mean field level. To this end, we use the measure ⟨S c (j, j ′ )⟩ [cf. Eq. (13)] to quantify quantum synchronization between two lattice sites j and j ′ . For the effective quantum model [cf. Eq. (8)- (9)] this measure can be calculated via the covariance matrix C and we use throughout this section the mean field amplitudes after relaxation as initial conditions. We focus here on the case of random initial conditions, that is the quantum fluctuations corresponding to the mean field dynamics shown in Figs. 2(a)-(d). In this way we are able to highlight the topological signatures while keeping the discussion precise. Furthermore, as we show in App. D, the significant results that occur for random initial conditions also carry on if eigenstates α (l) are chosen as initial conditions (similar as to the topological effects observed in the mean field dynamics). Finally, we test the robustness of the edge state synchronization when quantum fluctuations are included. In Fig. 4(e) we show ⟨S c (j, j ′ )⟩ between the two boundary oscillators (j = 1 and j ′ = 20) in the topological phase with coupling δλ = 0.8λ 0 as a function of the disorder strength r. Here, the overbar denotes that the quantity is averaged over 100 realizations of disorder. Similar as to the results of the mean field dynamics, also the quantum synchronization of the edges is robust for large amounts of disorder. As many fascinating topological effects occur in dimensions higher than one, we also investigate an example of a two-dimensional system. We consider here the breathing Kagome lattice [63][64][65][66], which is a natural extension of the SSH chain and a paradigmatic model of a so-called higher-order topological insulator (TI) [67]: While an ordinary d-dimensional TI, such as the SSH chain, exhibits d − 1 dimensional topological edge states, in a higherorder TI d − n dimensional topological boundary states emerge with n > 1 and the d − 1 dimensional topological edge states are absent. The reason for this peculiar phenomenon is that the boundary of a higher-order TI itself represents an ordinary TI. For the breathing Kagome lattice with a generalized chiral, time-reversal and particlehole symmetry (BDI class), the boundary states are zerodimensional corner states in the topological phase [68]. Recently, their robustness has been discussed and it has been realized that they are pinned robustly to zero only if the perturbations respect the generalized chiral and crystalline symmetries and the lattice connectivity [69,70]. This is the case considered in this work.
The tight-binding Hamiltonian of the lattice is given by where the two sums are over neighboring sites in the upward and downward triangles, respectively; see Fig. 5(a). Throughout we assume λ 1 ≤ 0 and λ 2 > 0 as we are only interested in the transition from trivial (λ 1 λ 2 < −1) to topological insulator phase (λ 1 λ 2 > −1). For λ 1 λ 2 > 0 there exists another phase transition to a metallic phase at λ 1 λ 2 = 1 2 [63], however, we do not consider this transition in the present work. In Fig. 5(a)

Topological effects on synchronization at the mean field level
In the following we investigate the mean field dynamics of the vdP network in the breathing Kagome lattice. To this end, the matrix governing the system dynamics without dissipation in Eq. (3) is given by where H Kag denotes the matrix corresponding to the Hamiltonian defined in Eq. (15). Following our general analytic considerations and the results of the previous Example (I), we expect that the oscillatory dynamics is dominated by the closed system dynamics given by Hamiltonian (16) , the oscillator dynamics follows the same principles as discussed previously, i.e. complete synchronized oscillations with common frequency ω 0 + µ (l) except of the three corners (index j = 1, 2, 3), which amplitudes oscillate with the intrinsic frequency ω 0 . The latter additionally exhibit a temporal amplitude modulation as a result of the finite amplitude of the initial state at the corners, i.e. A j (t) =Ā j cos{[ω 0 + µ (l) ]t} +Ā corner cos(ω 0 t), whereĀ j andĀ corner reflect the mixture of the two participating frequencies.
However, if the initial eigenstate is localized at the edges as shown in Fig. 5(d) we additionally observed excitation of all bulk lattice oscillating with a common frequency. However, their frequency does not match the corner nor the edge oscillators. Moreover, their amplitude is not constant but exhibits modulations over time. The reason for this effect is that the top band in Fig. 5(a) consists of strict bulk states without participation of any edge or corner lattice sites. Due to the nonlinearity in the dynamics, these strict bulk states may grow, yet, as a superposition of many of those leading to the observed amplitude modulations.
We are now interested whether the topological synchronization with random initial conditions also carries over to the two dimensional model. In Figs. 6(a)-(d) we show in the left panels the oscillator dynamics A j (t) as function of time for different coupling strengths with ordering of the oscillators according to Fig. 5. In the corresponding right panels we show the frequency spectrum of the particular realization of oscillator dynamics obtained via discrete Fourier transformation. In the trivial phase shown in Figs. 6(a) and (b) clear signatures of synchronization are again missing and the dynamics is governed by a randomly distributed superposition of many eigenstates oscillating with different frequencies.
The emergent complex oscillation pattern may exhibit temporal synchronized structures, however, these vanish again over longer times.
By contrast in the topological phase, shown in Figs. 6(c) and (d), the dynamics of the network appears more structured and synchronized. The reason for this is that fewer frequencies are available in the eigenspectrum of H Kag , and especially in the specific realization that there exists a dominant frequency as indicated by the large peak in the frequency spectra (see right panels). However, most notably the oscillators located at the three corners of the lattice (j = 1, 2, 3) are phase locked [cf. Eq. (11)] and oscillate with the intrinsic frequency ω 0 . Thus, as in the SSH chain, the topological character of the lattice is reflected in the amplitude oscillations even for random initial conditions. This also allows us to reconstruct the full eigenspectrum of the topological coupling matrix H Kag via discrete Fourier analysis of oscillator dynamics, which we show in Fig. 6(e) averaged over 10 realizations of initial conditions. The original eigenspectrum of H Kag has a highly degenerate eigenvalue spanning diagonally across, which is also observed in the reconstruction of the frequency spectrum as much larger amplitude (blue) than the other eigenvalues (pink).
Lastly, we test the robustness of the observed corner state synchronization. As the corner states of the breathing Kagome lattice are protected by a generalized chiral symmetry, we expect protection of the corner state synchronization. In Fig. 7 we show the frequency spectrum of the discrete Fourier analysis for the coupling λ 1 = −0.05λ 2 as a function of the disorder strength r for 10 realizations. The zero-energy corner modes are robust against disorder strengths as large as r = 0.4λ 2 before they are affected by the perturbations. In contrary, the bands spread over a wider range as the disorder strength is increased. Furthermore, this analysis also shows that the edge states are not topologically protected in the higher order TI. While the band frequencies are strongly affected by the disorder, the corner modes located at ω = ω0 are robust, even for large amounts of disorder. Parameters: κ1 = 5 ⋅ 10 −4 ω0, κ2 = 10 −2 ω0, λ2 = 0.25ω0, 10 realizations of random initial states for each data value of r.

Quantum signatures of topological synchronization
After we have observed in the previous section that topological protected synchronization of the mean field amplitudes also exists for a higher-order TI, we now analyze its quantum fluctuations quantified in terms of ⟨S c (j, j ′ )⟩ [cf. Eq. (13)] between two lattices sites j and j ′ . Similar as to the case of the SSH model we choose the mean field amplitudes after relaxation as initial conditions and focus on the case of random initial conditions corresponding to the dynamics shown in Figs. 6(a)-(d).
In Figs. 8(a)-(d) we show the time-averaged quantum synchronization measure ⟨S c (j, j ′ )⟩ for different couplings strengths (a) λ 1 λ 2 = −1.3, (b) λ 1 λ 2 = −1.2, (c) λ 1 λ 2 = −0.1, and (d) λ 1 λ 2 = −0.05. The ordering of the lattice sites corresponds to Fig. 6, especially j ′ , j = 1, 2, 3 are the corners of the lattice. We are mostly interested in the latter as they are expected to show quantum signatures of synchronization. In accordance with the classical mean field amplitudes of Figs. 6(a) and (b), the synchronization measure is almost uniform in the trivial phase shown in Figs. 8(a) and (b).
However, in the topological phase shown in Figs. 8(c) and (d) we observe that the oscillators located the corners (j ′ , j = 1, 2, 3) are significantly synchronized with one another as indicated by the large value of ⟨S c (j, j ′ )⟩ at the left bottom corner (highlighted by the pink dashed circle). As a reminder the quantity is bounded by 1 [cf. Eq. (12)]. Moreover, the quantum synchronization measure indicates that all three corner oscillators are synchronized with the edges and parts of the bulk. This feature is also present for the two edges of the SSH model in Example (I) [cf. Fig. 4]. By contrast, for any two bulk oscillators, ⟨S c (j, j ′ )⟩ remains similar to the previ- Finally, we test the topological protection of the corner state synchronization at the quantum level. In Fig. 8(e) we show ⟨S c (j, j ′ )⟩ between the corners [(j = 1, j ′ = 2), (j = 1, j ′ = 3) and (j = 2, j ′ = 3)] in the topological phase (δλ = 0.8λ 0 ) as a function of the disorder strength r. Here, the overbar denotes the average over 100 realizations of disorder. Consistent with the previous observations, also the quantum synchronization of the corners is robust for large amounts of disorder.

IV. DISCUSSION
An adequate formulation of topological insulators can be provided within the framework of solid state band theory and it is thus far from obvious whether and how their effects persist if affected by dissipation. Moreover, in the case investigated in this work dissipation represents a necessary resource to drive the system far away from equilibrium. Remarkably, despite the tremendous consequences open system conditions in combination with nonlinearities can have on the system dynamics, our examples show that topological features remain in such a scenario, which allows us to utilize them in our favor.
The observed topological boundary synchronization at the mean field level even if an eigenstate of the topological Hamiltonian is chosen as initial state is rather surprising and represents a key difference to systems without dissipation: In closed system dynamics, an initial eigenstate of the Hamiltonian will persist and will not mix with other eigenstates as they are orthogonal to each other. However, we observe in the nonequilibrium system that in the topological phase the edge oscillators are excited and their amplitudes grow, even though the initial eigenstate of the system Hamiltonian has very little amplitude at the edges. This may be understood as follows: Deep in the topological phase the oscillators located at the edges/corners are only weakly coupled to the rest of lattice. As their initial amplitude is finite it serves as a small perturbation away from their unstable fixed points, and since there exist an eigenstate localized at the edges this eigenstate will eventually grow. By contrast, oscillators in the bulk with small initial amplitudes are prevented from being excited because they require many other bulk oscillators to change their amplitudes accordingly. Interestingly, if the initial eigenstate is highly localized at the boundaries, they also serve as a perturbation for the bulk oscillators such that a combination of many bulk modes may be excited; see Fig. 5(d). We therefore argue, that the observed topological boundary synchronization at the mean field level is not restricted to the two specific examples investigated in this work, but also applies to other topological insulator lattices as the topologically protected edge modes are highly localized and all other bulk modes only have a small amplitude at the boundaries. Moreover, the synchronized edge modes inherit the topological protection known from closed systems with remarkable robust dynamics against local (symmetry-preserving) disorder and even random initial conditions. However, disorder in the natural frequencies of the oscillators on the considered lattices will destroy the observed edge synchronization. This is similar to topological insulators under closed system conditions, where the edge states are not robust against perturbations which break the underlying symmetries.
Often when fluctuations are considered synchronization is lost. Our results, however, show signatures of boundary synchronization beyond the classical mean field approximation. Furthermore, it remains unaffected for large amounts of local disorder in the couplings due to the underlying topology, which has the advantage that even if perfect fabrication of the lattice is impossible our findings can still be observed. This makes our results appealing for experimental realizations, two of which we discuss in more detail in the next section.
Lastly, let us highlight the ability to reconstruct the full eigenspectrum of the underlying topological lattice from the oscillator (mean field) dynamics alone. In our numerical studies with only 10 realizations the full spectra of the (closed) SSH and Kagome lattice could already be obtained. This holds the potential of a new experimental mechanism to measure eigenspectra of a topological systems from the dynamics of a nonlinear system far away from equilibrium and without preparation of an initial state.

V. EXPERIMENTAL PROPOSAL
There are many possible experimental platforms where topological synchronization may be observed in the mean field (classical) as well as the quantum regime. To realize the dynamics of a single vdP oscillator in a quantum system two different implementations have been suggested: The first one uses trapped ions [7,9] as an experimental platform while the second one focuses on optomechanics [9,10]. Both of these proposals bear the potential of vdP oscillator networks and thus to observe the previously discussed topological synchronization phenomena. Therefore, we summarize both of them and discuss their extensions towards topological lattices.
A motional mode of a trapped ion in the Lamb-Dicke regime and when the trapping potentials are tight with ground state g⟩ and excited state e⟩ is represented as an harmonic oscillator with frequency ω 0 . In order to fulfill the required conditions, i.e., negative and nonlinear damping, two side bands are excited simultaneously [71]: The first laser drives the transition g⟩ → e⟩ but blue detuned by ω 0 which absorbs one phonon after subsequent decay to the ground state [see Fig. 9(a)]. The second laser is double red detuned by −2ω 0 and excites to a state e ′ ⟩ such that after relaxation to the ground state two phonons have been emitted. The combination of these two processes approximately implement the dissipators of each oscillator in Eq. (2). Specific parameters for an implementation with 171 Yb + are provided in Ref. [7] and references therein. Additionally, the bosonic modes of neighboring trapped ions need to be coupled: By off-resonantly exciting an additional level e ′′ ⟩ via a blue sideband [see Fig. 9(b)] an effective Hamiltonian H = λ 1 (a † j a j+1 + a † j+1 a j ) for the phonon exchange between oscillator j and By exciting the blue sideband of the transition g⟩ → e⟩ negative damping may be realized (blue arrow). Simultaneously driving the double red sideband transition g⟩ → e ′ ⟩ allows one to realize the nonlinear damping (pink arrow). (b) Two modes may be coupled via off-resonant excitation of blue sidebands of g⟩ → e ′′ ⟩ or g⟩ → e ′′′ ⟩, which realizes the alternating couplings between nearest neighbors. (c) Optomechanical implementation: The so-called 'membrane in the middle'-setup allows to realize a vdP oscillator with nonlinear damping via a laser detuned to the red two-phonon sideband and negative damping by a laser on the blue one-phonon sideband. Mechanical coupling with alternating distances then allows to implement the topological lattice interactions.
j + 1 may be implemented. Using the same strategy but with an additional level e ′′′ ⟩ leads to a different coupling strength λ 2 . In this way the alternating interactions necessary for topological effects to emerge may be realized in trapped ion experiments.
Another experimental setup where our proposal could be implemented is provided by optomechanics as sketched in Fig. 9(c): The so-called 'membrane-in-themiddle'-setup [72][73][74] allows to realize two-phonon processes [75] as the cavity mode is parametrically modulated by the squared position of a movable membrane. The Hamiltonian of a single mechanical membrane inside a cavity driven on the red two-phonon resonance (in the good-cavity limit and focusing on the resonant terms) is given by H opt = g(a † a † b + h.c.) [75]. The dynamics of the coupled system is described by the master equatioṅ Adiabatically eliminating the optical mode and defining κ 2 ≡ g 4 γ c results in the two-phonon dissipator of Eq. (2). However, the other two linear dissipators in Eq. (17) are related via the thermal phonon occupationn, such that one-phonon gain cannot be larger than the one-phonon loss. Therefore, another laser driving the cavity on the blue one-phonon sideband is needed [9,10]. Nevertheless, compared to Eq. (2) an additional one-phonon loss term γ(1 +n)D[a] is present. If we assume γ(1 +n) ≪ κ 1 , the influence of this additional term can be completely neglected at the mean field level as it only results in an effective damping rateκ 1 = κ 1 − γ(1 +n). Also our results in the quantum regime remain unchanged if one would include such a linear loss channel (see App. E).
When two membranes share a common support they are mechanically coupled expressed via the Hamiltonian H =λ j x j x j+1 where x j = 2mω 0 (a † j + a j ) with oscillator mass m. After the rotating-wave approximation and with the definition λ j = 2mω 0λj , the Hamiltonian becomes equivalent to the coupling Hamiltonian (1). The interaction strengths between two oscillators may be modified by altering the distance between neighboring membranes; see Fig. 9(c) for an example of the SSH model with staggered distances d 1 and d 2 . Thus, also optomechanics has the potential to observe topologically protected synchronization in a network of quantum vdP oscillators.

VI. CONCLUSIONS
We investigate the interplay of topology and synchronization in a network of coupled quantum van der Pol oscillators simulating different topological insulator lattices. We show via a linear stability analysis that the dynamics of the resulting topological lattice of oscillators at the mean field level is governed by the eigenvalues of the topological Hamiltonian and thus reflects the features of the underlying topology even though the system is highly nonlinear and far away from equilibrium. Furthermore, we derive an effective quantum model which takes quadratic quantum fluctuations about the classical trajectories into account in order to investigate quantum signatures of topological synchronization beyond the mean field approximation. For two specific examples of topological insulator models in one and two dimensions we demonstrate that in the nontrivial phase, synchronization at the boundaries is always present independent of initial conditions, and that it inherits the protection against perturbations from the underlying lattice structure. In terms of a Fourier analysis of the oscillations, we are able to reconstruct the full topological eigenspectrum of the system, which not only represents a possible route to observe topological synchronization of boundary modes experimentally, e.g. in trapped ions or optomechanics, but an additional opportunity to measure topological eigenspectra solely from the dynamics of a highly nonlinear and open system.
Researchers and engineers make great efforts to fabricate dynamical systems which are nearly identical in order to facilitate the emergence of synchronized collective behavior in large networks. However, our work demonstrates a general advantage of topological lattices in the design of potential experiments and devices as fabrication errors and longterm degradation are circumvented in this way. This is especially important in networks where specific nodes need special protection. While the examples investigated in this work posses zero dimensional protected boundary states, our work can easily be extended to host higher dimensional topologically protected states for additional robust network nodes. Synchronization is desirable in situations where high oscillating power, strong coherence, or low phase noise are needed, such as lasers [76], phase-locked loops [77], Josephson junction arrays [78,79], spin-torque resonators [80], quantum heat engines [81] or power grids [82]. Even today, the originally observed phenomenon of clock synchronization remains a crucial application for modern communication networks [83,84] and has recently been extended to quantum networks and quantum key distribution protocols [85,86]. All of these examples require the synchronized behavior to be robust to fulfill their desired purpose and will benefit from the application of topology. Given the universality of the concept of combining nonlinear dynamics in open quantum systems with topological phases of matter, we expect that our approach could be successfully applied also to other systems where robust dynamics is crucial. where the transformed system Hamiltonian is given bỹ and the two dissipators by , We now evaluate the effect of applying the displacement operator on the different terms appearing in Eq. (A2). For transformed HamiltonianH S we obtaiñ For the dissipative term proportional to κ 1 we obtain and for the dissipative term proportional to κ 2 we obtain We now realize that as long as condition (3) i.e., the mean field equation, is satisfied, the linear terms vanish. If we neglect higher order terms, the quantum fluctuations are governed by master equation of Lindblad form: with effective Hamiltonian In this section we specify the entries of the matrices B and D governing the equation of motion of the covariance matrix C:Ċ The entries of the matrices are given through Eq. (8).
Specifically, B is a block matrix, where the blocks on the diagonal B jj take the form and the off-diagonal block matrices B jk for j ≠ k the form Furthermore, D is a diagonal block matrix with entries D jj on the diagonal given by

Two coupled vdP oscillators
The full quantum model of two coupled vdP oscillators is given by Eq. (2) with j = 1, 2 and system Hamiltonian Following Ref. [7], we characterize the two oscillator system in terms of the two mode Wigner function W (α 1 , α * 1 , α 2 , α * 2 ) and integrate α 1 , α 2 and ϕ 1 + ϕ 2 , such that the Wigner function is a function of the phase difference ϕ 1 − ϕ 2 alone. In Fig. 11 we show the Wigner function W (ϕ 1 − ϕ 2 ) − π 2 as function of the phase difference between the two oscillators without coupling (blue) and with coupling λ = 0.5ω 0 (pink). We choose the latter value as it is the largest coupling we consider in the main text. While the Wigner function is completely flat in the absence of interactions, it shows two peaks at ϕ 1 − ϕ 2 = 0, π corresponding to in-phase and anti-phase synchronization respectively. Due to the strong quantum noise in the regime considered in this work, the peaks are quite small, which has also been reported in Ref. [7]. For two coupled vdP oscillators the effective dynamics are given by Eq. (9) with effective Hamiltonian where the mean field amplitudes α j (t) follow the dynamicsα 1 = −iω 0 α 1 − iλα 2 + κ 1 2 α 1 − κ 2 α 1 2 α 1 , In the case without interaction (λ = 0) the mean field amplitude of each oscillator will eventually reach the steady stateᾱ j (t) =Ā exp(−iω 0 t + ϕ j ) with an arbitrary phase difference ϕ 1 −ϕ 2 set by the initial conditions. In the case of large interaction strengths the two oscillators will synchronize, however, either in-phase or anti-phase, again set by the initial conditions. It is thus necessary to average over the initial conditions in order to properly compare it to the full quantum dynamics. As discussed previously, the initial condition of the effective density matrix is given by α (t = 0) = ⊗ j 0 j ⟩ ⟨0 j . In Fig. 11(b) we show the Wigner function W (ϕ 1 − ϕ 2 ) − π 2 as function of the phase difference between the two oscillators for λ = 0 (blue) and λ = 0.5ω 0 (pink). Similar to the full quantum model shown in Fig. 11(a), we observe that if the interaction is turned off the Wigner function is completely flat, while it is peaked at ϕ 1 − ϕ 2 = 0, π with the coupling turned on, however, compared to panel (a) the peaks are smaller. The reason is that the effective model overestimates the quantum fluctuations leading to less pronounced synchronization. This is consistent with the observations of Ref. [7], where for κ 2 = 10κ 1 the height of the peaks is in the same order of magnitude.
In the main manuscript we quantify the level of synchronization between two vdP oscillators in terms of the synchronization measure defined in Eq. (13) rather than the Wigner distribution. We therefore also compare this quantity for the full quantum model and the effective model. In Fig. 12 we show ⟨S c ⟩ as a function of the coupling strength λ for two coupled vdP oscillators without (blue) and with approximations (pink). In both cases, ⟨S c ⟩ is finite even though the oscillators do not interact, which is due to the fact that both oscillators are equivalent. As the interaction increases, also ⟨S c ⟩ increases until it quickly approaches a plateau. In accordance with the previous observations regarding the Wigner function, also the synchronization measure is overall smaller for the effective model than for the full dynamics. and furthermore exhibit in-phase synchronization of the mean field amplitudes; see Figs. 1(d) and (e).

Appendix E: Including an additional linear loss channel
For the experimental realization of a vdP oscillator with optomechanics as discussed in Sec. V, there exists an additional linear one-phonon loss channel. However, we will show in this section that such an additional dissipative term does not affect the results discussed in the main text as long as it is small compared to the linear gain channel. Including this additional loss term, the Lindblad master equation is given bẏ Within in the mean field approximation, the term proportional toγ only appears as a shift of the dissipation rate κ 1 , i.e., with the definitionκ 1 = κ 1 −γ the dynamics is fully equivalent to Eq.(8) assuming thatγ ≪ κ 1 . The latter condition must be fulfilled in order to counteract the nonlinear damping (κ 2 ) and not simply decay into the ground state.
In terms of the quantum fluctuations about these mean field amplitudes, the additionalγ term appears in the matrices B and D governing the dynamics of the covariance matrix; see Eq. (10). Similar to the mean field approximation, in B this will only result in an effective dissipation rate κ 1 →κ 1 = κ 1 −γ in Eq. (B2). However, the sign is flipped for the effective dissipation rate in D [cf. Eq. (B4)], which is then given by We have checked numerically that forγ = 0.1κ 1 there are no noticeable changes in the dynamics or synchronization behavior, such that we are confident that in fact a linear dissipation channel does not alter the discussed results.