Topology protects chiral edge currents in stochastic systems

Constructing systems that exhibit time-scales much longer than those of the underlying components, as well as emergent dynamical and collective behavior, is a key goal in fields such as synthetic biology and materials self-assembly. Inspiration often comes from living systems, in which robust global behavior prevails despite the stochasticity of the underlying processes. Here, we present two-dimensional stochastic networks that consist of minimal motifs representing out-of-equilibrium cycles at the molecular scale and support chiral edge currents in configuration space. These currents arise in the topological phase due to the bulk-boundary correspondence and dominate the system dynamics in the steady-state, further proving robust to defects or blockages. We demonstrate the topological properties of these networks and their uniquely non-Hermitian features such as exceptional points and vorticity, while characterizing the edge state localization. As these emergent edge currents are associated to macroscopic timescales and length scales, simply tuning a small number of parameters enables varied dynamical phenomena including a global clock, dynamical growth and shrinkage, and synchronization. Our construction provides a novel topological formalism for stochastic systems and fresh insights into non-Hermitian physics, paving the way for the prediction of robust dynamical states in new classical and quantum platforms.


I. INTRODUCTION
Why are biological functions carried out so robustly, even when the underlying components are stochastic in time and randomly distributed in space? Living systems can have stable properties that endure for time-scales much longer than the lifetime of the underlying constituents, that contribute to memory and adaptive processes [1,2]. The emergence of stable and reproducible timescales and length scales is also a key goal in the design of synthetic biological systems [3,4] or in the engineering of reconfigurable materials e.g. through dissipative self-assembly [5,6]. However, these strongly out-ofequilibrium systems often lack a comprehensive theoretical framework, which prevents us from understanding or describing these processes [7].
As a step towards answering this question, we present stochastic models that exhibit emergent dynamical phenomena reminiscent of many observations prevalent in biology [8][9][10][11][12], such as a global clock, dynamical growth and shrinkage, as well as synchronization. In our models, all these phenomena hinge on the emergence of a topologically-protected chiral edge state in the stochastic system. To illustrate this connection, we use the canonical example of a topological phase: the quantum Hall effect [13], see Fig. 1(a). Here, electrons make cyclotron orbits in the bulk that correspond to unidirectional edge states at the boundaries. The system transport is hence dominated by the propagating edge states which are exponentially localized at the boundaries and protected from disorder and perturbations. The robust-ness of these edge currents make them a desirable feature for the support of stable emergent phenomena. However, they have yet to be realized in biochemical systems, i.e. systems governed by memoryless (Markovian) classical stochastic dynamics. While topological states showing stationary polarization have been recently reported in one-dimensional stochastic systems [14,15], states with propagating edge currents have not yet been reported.
The two-dimensional networks we introduce are constructed from simple repetitive motifs, which correspond to out-of-equilibrium cycles at the molecular scale, and form the analog of cyclotron orbits in the quantum Hall system. In biochemical networks, such microscopic transitions are common due to out of equilibrium transitions that consume a fuel such as ATP or GTP. Many of these appear to leave the system unchanged while consuming energy, and having been dubbed "futile cycles" are ubiquitous in biology [16,17]. Beyond the analogy to quantum Hall physics, the chiral edge currents in our system imply motion along the boundaries of configuration space rather than real space. As such, they enable oscillations (e.g. cyclical conformational changes of a protein complex, or assembly and disassembly of a biopolymer) governed by physical constraints in the system rather than the specific timescales of the underlying microscopic transitions, which do not need to be fine-tuned [1,2,7]. Hence, these systems constitute an excellent example of the "structure determines function" paradigm of biology. Besides their possible relevance to biological oscillators, our models could guide the engineering of synthetic nonequilibrium machines that perform work (e.g. stochastic low Reynolds number swimmers [18,19]).
The out-of-equilibrium or nonreciprocal stochastic transitions render the transition matrix non-Hermitian. states of matter was initially developed for quantum electronic systems [20][21][22][23][24][25][26][27] and later extended to classical systems such as mechanics [28][29][30], photonics [31][32][33][34], acoustics [35,36], electrical circuits [37,38], active matter [15,39,40], and population dynamics [41,42], this formalism has traditionally been reserved for Hermitian systems that respect energy conservation and isolation from the environment. As our systems are strongly dissipative and break detailed balance at the microscale, we build a topological description that is explicitly non-Hermitian. Non-Hermitian physics is an area of great recent interest [22-26, 33, 34, 43-48], and can demonstrate unique properties that have no analog in Hermitian systems. A celebrated example is that of exceptional points: singularities where eigenstates coalesce [23,24,[43][44][45][46], in contrast to degenerate level crossings in Hermitian systems where eigenvectors remain distinct. A second example is vorticity, a topological winding of eigenvalues in complex space [22], differing from topological invariants in Hermitian systems which are typically defined using eigenvectors, as in these systems the eigenvalues are purely real. We show that our models support the topological Zak phase [21, 27, 31-33, 35, 37], but only support propagating edge currents in the non-Hermitian case. Fur-ther, they exhibit the above-mentioned non-Hermitian features such as exceptional points and vorticity, creating a novel theoretical formalism for stochastic systems distinct from previous proposals that do not show these features [14,15].
The paper is organized as follows. In Sec. II, we describe the minimal motifs that can be used to construct two-dimensional stochastic networks with emergent chiral edge currents (Sec. II A), and show how global cycles arise using stochastic simulations and analytical calculations (Sec. II B). In Sec. III, we generalize the model to explore the transition from Hermiticity to non-Hermiticity, and explain how the emergence of edge states can be understood as a topological transition. This is a consequence of the Zak phase (Sec. III A) and uniquely non-Hermitian topological invariants such as vorticity (Sec. III B), which results in exponential localization at the edges (Sec. III C). Finite-size effects and the effects of disorder on topological protection are discussed in Sec. III D. In Sec. IV, we apply our model to specific examples in soft and living matter, such as biomolecular oscillators (Sec. IV A) and stochastic swimmers (Sec. IV B), and further extend the model to account for other biologically-relevant features, such as asymmetry in configuration space (Sec. IV C) and shared boundaries between different subsystems (Sec. IV D), which lead to new behaviors such as dynamical instability and synchronization. Finally, in Sec. V we put our work in the context of previous proposals for topological protection in related systems. We end with a discussion of the implications and future directions for our work.

A. Minimal motifs
We consider discrete stochastic processes that operate in a two-dimensional configuration space, i.e. for which the state of the system is determined by two integers (x, y). These two numbers could represent, for example, the state of a biopolymer assembled from two types of monomers X and Y, or from monomers of a single type X but which can be modified (e.g. via phosphorylation). They could also represent two types of modifications applied to the monomers that make up a fixed-size structure such as a protein complex. Implementing transitions between contiguous (x, y) states results in a lattice-like description of the system. Such a lattice will have boundaries or "edges" representing the physical constraints in the system, for example 0 ≤ x ≤ N x and 0 ≤ y ≤ N y where N x and N y represent e.g. the number of X and Y monomers available for binding, or the number of binding sites for X and Y in a protein complex. More elaborate constraints can also arise, such as 0 ≤ y ≤ x if y describes the number of monomers in a biopolymer that have undergone some modification out of a total of x.
A simple implementation of microscopic out-ofequilibrium "futile" cycles, reminiscent of cyclotron orbits in the quantum Hall effect [ Fig. 1(a)], can be achieved in a system with four internal states (A,B,C,D) and four external transitions We have defined these transitions such that cycles are clockwise in (x, y) space, without loss of generality (reversing all arrows would give counter-clockwise cycles). In the absence of any other transitions, the system will be trapped in such cycles and will not explore the available configuration space. However, the system can break out of a cycle if internal states undergo decay, with four  Fig. 1(c)]. Initially, motion is diffusive, until the system encounters the edge at y = 6, after which it shows persistent oscillations in both x and y. That the oscillations correspond to counter-clockwise edge currents in (x, y) space is clearly seen in (b) which depicts the same trajectory, but in two dimensions (color represents time). (c) The probability distribution in (x, y) space, obtained from simulations, and (d) the steady-state probability distribution in full configuration space, obtained from direct solution of the master equation, both show strong accumulation of probability at the edges. The four arrows in (d) indicate the direction of the edge currents. (e) Edge currents and the resulting cycles are robust with respect to the shape of the boundaries. Parameters used: in (a-c) γin = 10 −3 γex, in (a-d) system size Nx = Ny = 6, in (c) 10 7 steps starting from a random initial state were used. See also Movie 1 in the Supplemental Material [49].
which enable effective diffusion over the bulk of (x, y) space over time. A possible implementation of these 8 transitions for a biopolymer is shown in Fig. 1(b). The resulting lattice can be embedded in the plane, as shown in Fig. 1(c).
Cycles with only three internal states are possible if we allow for diagonal transitions in (x, y), such as This gives cycles with only  three external transitions, such as Including three internal decay transitions for the same reason as above (see Fig. 1(d) for a possible implementation of the 6 resulting transitions for a biopolymer), we again build a lattice that can be embedded in the plane, e.g. as a Kagome lattice [ Fig. 1(e)].
Note that, except where explicitly stated, we will consider rectangular geometries in the 4-state model, corresponding to the constraints 0 ≤ x ≤ N x and 0 ≤ y ≤ N y and thus a lattice with a total number of 4(N x +1)(N y +1) states; and triangular geometries in the 3-state model, corresponding to the constraints 0 ≤ x ≤ N x and 0 ≤ y ≤ x, which gives a lattice with a total number of 3(N x + 1)(N x + 2)/2 states. In both cases, the indexing of states starts from (x, y) = (0, 0).

B. Chiral edge currents and global cycles
Inspection of the lattices in Fig. 1(c,e) suggests that persistent counter-clockwise trajectories of the system along the edges are possible if γ ex γ in , i.e. if the external transition is more likely than the internal one when both are possible (e.g. at a B state in the bottom edge) so that the system remains on the edge. These trajectories are then analogous to the "skipping orbits" in the quantum Hall effect, see Fig. 1(a). Stochastic simulations of both the 4-state (Fig. 2) and the 3-state (Fig. 3) models confirm this expectation. Starting from a state within the bulk of the lattice, the system initially displays local clockwise cycles (driven by γ ex ) interspersed with occasional sideways steps (driven by γ in ), leading to diffusive motion in the bulk. Once the system reaches any state on the edge, however, persistent motion on the edge leading to counter-clockwise cycles along the boundaries of the system is observed, see Fig. 2(a,b) and Fig. 3(a,b). Over time, the probability of finding the system at the edge is significantly larger than in the bulk, see Fig. 2(c) and Fig. 3(c), which correspond to histograms obtained directly from the combined x(t) and y(t) trajectories. Direct solution of the steady state probability of the full master equation of the system confirms this result, see Fig. 2(d) and Fig. 3(d) and further shows a more detailed structure for the probability of different internal states (or sites) on any given edge cell. This more detailed structured is also recovered from the stochastic simulations (data not shown). As long as γ ex γ in , the edge cycles are robust to variations in the system size or shape, provided that the directionality of the lattice edges is preserved, see Fig. 2(e) and Fig. 9. We note however that, in the case of Kagome lattices with nontriangular geometry, transitions cutting across concave corners must be absent for edge states to be fully pro-tected, see Appendix A and Fig. 9 for more details. Because the global cycles occur in configuration space, they can describe cyclic changes of a molecular system. In a system of variable size, such as a biopolymer composed of two types of monomers X and Y and governed by the 4-state model, a global cycle around the system boundary would imply sequential assembly of X monomers, followed by assembly of Y monomers, followed by disassembly of the X monomers, and finally disassembly of the Y monomers, leading back to the initial state. The maximum length in X or Y obtained would only be limited by the availability of monomers. A similar growth-shrinkage cycle can be obtained in the 3-state model, where now y could represent e.g. the number of monomers that have been dephosphorylated out of a total of x. A full cycle would then involve assembly of phospohorylated monomers, dephosphorylation of the monomers, and disassembly of the dephosphorylated monomers. Alternatively, the configuration space could correspond to a system with a fixed number of components, such as a protein complex, but whose components can undergo two types of conformational changes. Such a model, which can describe a biomolecular clock, will be described in Section IV A, where we also explore the period and coherence of the emergent oscillations in more detail.
The persistence of edge trajectories can be understood quantitatively. The probability of remaining L steps along the edge and then "unbinding" from it is given by P (L) = 1 − γin γin+γex L γin γin+γex , which results in an average run length L = ∞ L=0 LP (L) = γ ex /γ in . Thus, for γ ex = 10 3 γ in and N x = N y = 6 as in Fig. 2, we expect the system to perform 10 3 /(6 · 4) ≈ 42 full cycles on average before unbinding. Even then, the system is likely to encounter the edge again soon after and thus undergo a new run along the edge. Moreover, we can analytically obtain (see Appendix A) the stationary probability distribution of the system, both in the 4-state and the 3-state models, and find that probability accumulates in the edge sites that precede an internal transition [e.g. C sites at the bottom edge; see Fig. 1(c,e)], which have stationary probability p C = γin+γex γin p b , where p b is the probability corresponding to all other sites, including bulk sites as well as edge sites that precede an external transition (e.g. B sites at the bottom edge). This corresponds to the results in Fig. 2(d) and Fig. 3(d), as well as Fig. 2(e) and Fig. 9(b). Importantly, these results imply that probability is infinitely localized at the edge sites in the steady state, in the sense that there is not a gradual decay to the bulk probability as we move away from the edge sites. Summing up the probability of all edge sites, we can obtain the overall probability P edge of finding the system at the edge at any time, or equivalently, the fraction of time that the system spends at the edge. For a square 4-state system of size N x = N y = N , in the limit γ ex γ in , we find P edge γex/γin N +γex/γin (see Appendix A). In the example of Fig. 2, this implies that the system spends ≈ 99.4% of the time at the edge.

III. TOPOLOGICAL PROTECTION OF THE EDGE STATES
A. Berry connection and transition to the Zak phase The emergence of edge states that dominate the system dynamics is a hallmark of a topological phase, which we can see upon analyzing the Master equation that describes stochastic systems, d dt p = Wp. Here p is a vector of the probabilities of being in each state, and W is a real matrix specifying the transition rates [50]. The lattice structure of the bulk transitions in W allows for the calculation of the Berry connection of W, which determines the topological index of the system, i.e. whether the system is in the Zak phase. When this index is 0, the system is in the trivial phase, in which case the system will mostly remain in the bulk. In contrast, when the index is π, long-lived edge states will appear. This relation between the bulk topological phase and its edge properties constitutes the celebrated bulk-boundary correspondence [13]. The link between a particle moving in a lattice and its Berry connection was first formulated by Zak for electrons in a crystal lattice [21]. Later, this link was extended to many other systems including nonelectronic ones such as photonics [31][32][33], mechanics [29], acoustics [35], and electrical circuits [37]. While the lattice coordinates represent real space in these previous works, in ours they represent configuration space. Remarkably, W in our system is a 2d non-Hermitian generalization of the Su-Schrieffer-Heeger (SSH) model for polyacetylene [20], up to a diagonal matrix (see Appendix B). When the lattice has periodic boundary conditions (PBC) this diagonal matrix is proportional to the identity matrix, and the eigenvectors of W are exactly those of a 2d SSH model, hence they have the same Berry connection (see Appendix C). Integration of the Berry connection over reciprocal space is quantized in our system due to the presence of inversion and sublattice (or chiral) symmetries [21,26,27,33], and can support the Zak topological phase.
The 2d SSH model on the square [27] and kagome [36] lattices have been previously studied in the Hermitian limit, where the system is in the Zak phase when γ ex > γ in and exhibits edge states. To understand how these properties extend into the non-Hermitian case which characterizes our out-of-equilibrium system, we generalize the phase space. We introduce transitions in the reverse direction from γ ex and γ in , which we call γ ex and γ in respectively [ Fig. 4(a)]. In Fourier space, the spectrum of W for the 4-state model takes the form (see Appendix B) To simplify notation, we introduce two parameters, the ratio r and the chirality c, where The ratio r weights the relative strength of internal and external transitions, with r > 1/2 (r < 1/2) when external (internal) transitions are stronger. The chirality parameter c weights the relative strength of forward and reverse transitions, where we note that γ in + γ ex = cγ tot and γ in + γ ex = (1 − c)γ tot . For c > 1/2, external and internal cycles are biased in the clockwise and counterclockwise direction, respectively, and edge currents are counter-clockwise. For c < 1/2, the chiralities are opposite, with external and internal cycles biased in the counter-clockwise and clockwise direction, respectively, and edge currents being clockwise. The value c = 1/2 corresponds to the Hermitian, achiral case with equal forward and reverse rates for all transitions. Finally, c = 1 describes the fully chiral case with only forward rates studied above (see Fig. 1), and c = 0 is also fully chiral but with opposite chirality.
We can now examine if the system remains topological as chirality is tuned from the Hermitian case to the fully chiral case, by verifying if the system bandgap remains open throughout. Using Equation 1, we indeed note that the bandgap remains open as we vary c from 1/2 to 1 above a critical value r * [ Fig. 4(b)]. We find that r * = 1/2 for both the purely Hermitian case [27] (c = 1/2) and the fully chiral case (c = 1 and c = 0), see Appendix C and Fig. 10 there. However, stronger external transitions are needed for the transition in between these values, similar to when finite-size effects come into play [27]. The value r * (c) can be obtained from the maximum value of r for which a(k) 2 = b(k) has a solution. The global maximum across all values of c is r * max = 0.59, which occurs at c = 0.7 (and c = 0.3).
In the topological phase, i.e. when r > r * , the system exhibits edge states in open boundary conditions (OBC) as expected from the bulk-boundary correspondence. These states are colored in red in Fig. 4(c), forming two rings that shrink into points as r → 1 . The larger of these rings shrinks towards the ground state at E = 0, creating a large number of edge states that are increasingly close to the ground state as r → 1. Note that

Right edge
Left edge the total number of edge states scales with the number of sites on the edge, e.g. 4(N x + N y − 1) in a rectangular geometry.

B. Uniquely non-Hermitian features: propagating edge states, exceptional points and vorticity
We find that our system exhibits several uniquely non-Hermitian properties in the chiral case, which is also when the edge states change qualitatively from being stationary polarization to chiral currents. Such non-Hermitian features include exceptional points (EPs) [23,24,[43][44][45][46] or topological invariants without a Hermitian counterpart [22,48]. These unique non-Hermitian properties emerge only in the finite Zak phase and most notably in the case of OBC. In the fully chiral case, the spectrum of the system with PBC (grey lines) is symmetric with respect to r * = 1/2 (where γ ex = γ in ) [ Fig. 4(d)].
In OBC (blue lines), the spectrum changes radically past r * , with many states coalescing towards E = 0 and some towards E = −γ tot as r → 1, indicating the existence of EPs [24,43] at r = 1 (see dashed circles in Fig. 4(d)). The transition at r = r * coincides with when the system is just as likely to unbind as to remain on the edge at every step, with the average run length of L = 1, whereas the limit of r → 1 corresponds to when the system spends all of its time on the edge. Lastly, varying the edge transition rate γ 0 ex (see Appendix B) interpolates between PBC (γ 0 ex = γ ex ) and OBC (γ 0 ex = 0). In the Zak phase, exceptional points emerge in the spectrum at γ 0 ex = γ in [Fig. 4(e)]. Another uniquely non-Hermitian feature emerges in the edge states in the strongly chiral limit, that of vorticity [22]. Using a half-periodic geometry with OBC in x and PBC in y (Fig. 5), we can calculate the band structure (eigenvalues) of W along the reciprocal lattice index k y . In the Hermitian limit of c = 1/2, all the bands have similar amount of group velocity d(ReE)/dk y and are completely real, see Fig. 5(a,b) (left). When Hermiticity is broken and the system becomes increasingly chiral, bands localized on the edge emerge, which also exhibit the largest group velocity, see Fig. 5(a) (center and right). Many bands also have an imaginary component which increases in magnitude with c, see Fig. 5(b) (center and right). As the group velocity describes the propagation speed of probability disturbances, the increase of group velocity with chirality implies a decrease in the period of oscillations around the edge of the system with increasing chirality, which is supported by the explicit calculations in Section IV A.
These edge states demonstrate a topological transition with increasing chirality. Part of these bands are completely real (red) and distributed in equal amounts on the left and right edges [red in Fig. 5(c)]. As chirality increases, they develop growing imaginary components to become localized on the left (green) and right (orange) edges respectively. In the strongly non-Hermitian limit, a transition occurs when the two pairs of edge bands touch in real space, exhibiting a doubled periodicity of 4π and vorticity of ν = 1/2. Vorticity is a uniquely non-Hermitian topological invariant describing the winding number of a pair of bands in the complex plane [22]: where Γ is a closed loop in reciprocal space, and m, n are band indices. In our system, ν indicates the strongly propagating nature of the edge states, taking the value of 0 at c = 0.8 and 1/2 at c = 1.
To summarize the topological properties of the system, we find that (i) topological edge states occur when r > r * (c) (the shape of this boundary can be seen in Figs. 6(f) and 12); (ii) when c = 0.5, the system is Hermitian and no chiral probability currents are observed at the boundary; and (iii) for c > 0.5 probability currents at the edge emerge in a counter-clockwise direction; whereas for c < 0.5 currents emerge in the clockwise direction. Topologically, the emergence of protected chiral edge currents is manifested in the transition from zero vorticity ν = 0 to non-zero vorticity ν = 1/2 and ν = −1/2 as chirality increases or decreases away from c = 0.5, respectively.

C. Edge localization in steady-state
Properties of the non-Hermitian system can also be analyzed using a transfer matrix approach, which probes the steady state of the full transition rate matrix W (see Appendix D), in the general case of arbitrary c and r. This yields the probability accumulation and current along the edge in the full phase space of c and r, see Fig. 12 in the Appendix. In particular, we find that both the probability accumulation and the current vanish in the achiral, Hermitian case with c = 0.5, as well as in the limit r → 0. The current also vanishes in the limit r → 1, at which point the sites along the edge become disconnected because γ in , γ in → 0.
Moreover, the transfer matrix analysis shows that the edge states are exponentially localized, in the sense that the probability accumulation with respect to the bulk probability in a cell situated n cells away from the edge δP n decays as δP n = δP 0 α n , where α is a decay constant satisfying 0 ≤ α < 1 that depends on both c and r. Importantly, in the limit of a fully chiral system with c → 1 or c → 0 we find α = 0, and thus recover the result of infinite localization at the edges described in Sec. II B.

D. Identifying topology in realistic systems
While some of the results above are obtained in a periodic setting, i.e. representing an infinite system without disorder, they are directly applicable to realistic systems which are typically finite and may have some amount of disorder in the transition rates. Thanks to the bulkboundary correspondence [13,51], our results for the topology of W in a periodic setting will manifest as edge state properties when a finite system is considered, or when an edge is introduced in an otherwise open system. Further, since our edge states exponentially decay into the bulk as discussed in the previous section, finite size effects would not alter the bulk properties in systems that are larger than the exponential decay length. In practice, this allows the bulk-boundary correspondence to hold even for small systems with just a few unit cells, for instance in the fully chiral limit where the edge states are infinitely localized. This permits the calculation of the topology of W for smaller systems, where we expect an edge state at the boundary of the system in the topological phase, or at the boundary between two subsystems that are in topological and trivial phases respectively.
The bulk-boundary correspondence also remains robust to weak disorder, as has been found from both analytical [51] and numerical [52] calculations. Specifically, this is true when the disorder does not does not cause the bandgap to close, which preserves the topological phase [52], e.g. if disordered transitions are too slow to cause scattering between the edge and bulk states. This can be checked by direct calculation of the disordered bulk W. Alternatively, inspection of the eigenstates of W in a finite geometry in the presence of disorder can reveal whether the long-lived states remain on the edge or have merged with the bulk states. The presence of strong disorder in such systems merits systematic exploration in future work. showing that coherence is maximized in the fully chiral limit (c → 0 or 1) and deep in the topological phase (r → 1). The white curve represents the boundary r * (c) separating the topological and trivial phases of the system. For each curve in (b) and each data point in (d,e) we used 10 7 stochastic simulation steps starting from a random initial state.

A. Global clock
As anticipated in Section II B, the two-dimensional configuration spaces described here may represent a protein complex with a fixed number of components that can undergo two types of conformational changes. As an example, in Fig. 6(a) we show how our 4-state model with N x = N y = 6 would support the dynamics of an allosteric model for a hexameric biochemical oscillator such as the KaiABC system. Typical models for such conformational oscillations [8][9][10] rely on the concerted or Monod-Wyman-Changeux paradigm [53] of allosteric regulation which restricts the configuration space to states in which either all or none of the monomers in a complex have undergone a conformational change. The latter traditionally stands in opposition with the sequential or Koshland-Némethy-Filmer paradigm [54], which allows individual monomers to separately undergo conformational changes. Interestingly, our system shares features of both paradigms: while all states in the full two-dimensional configuration space are allowed, as in the sequential model, the system only visits a strongly limited subset of these states, namely those at the edge, as in the concerted model. Moreover, this property arises from the repetition of simple, local reaction motifs without the need of fine-tuning, e.g. of some transitions to be markedly different from others at particular points in configuration space.
In order to obtain a functional clock, it is important that oscillations have a reproducible period. We have characterized the period of oscillations for our 4-state model with N x = N y = 6 in the fully chiral case (c = 1, i.e. γ ex = γ in = 0), as a function of the ratio γ ex /γ in . Note that γ ex /γ in = r/(1 − r), where r is the ratio parameter introduced in the previous Section III A. In Fig. 6(b), we show the distribution of the times τ 1 needed to complete a full cycle, defined as the time between two consecutive crossings of a reference line which is marked as the dashed line in the inset of Fig. 6(b), obtained from stochastic simulations, for several values of γ ex /γ in . In the topological regime with γ ex /γ in 1, we find a strongly peaked distribution, implying a well-defined oscillation period. The average oscillation period can be taken to be T ≡ τ 1 = J −1 | , where J | is the average current through the same line that is used to define the cycle time. As the ratio γ ex /γ in is decreased, this distribution broadens until, for γ ex /γ in 1, the distribution is no longer peaked and instead shows a standard deviation comparable to or larger than its mean. From the mean and variance of the cycle time we can define a measure of coherence N ≡ τ 1 2 /var(τ 1 ) (inverse to the Fano factor) [55,56], which represents the average number of oscillations that the system will undergo before an error of plus or minus one oscillation has accumulated and coherence is lost. The period and coherence obtained from stochastic simulations in this way are shown as the red circles in Figs. 6(d) and (e), respectively. Moreover, we can calculate the average current J | in steady state explicitly for a square 4-state model with N x = N y = N , see Appendix A, and thus the oscillation period as which is plotted as the red solid line in Fig. 6(d) and matches the simulation results. In the limit γ ex /γ in N , the period tends to T = 4(N + 1)/γ in (or T = 28/γ in for N = 6), which represents the fact that in this limit the edge states become isolated from the bulk, the waiting times for external transitions become negligible with respect to those for internal transitions, and the system thus becomes equivalent to a unidirectional cycle with transition rate γ in between 4(N + 1) = 28 states.
An alternative way to define the period and coherence of oscillations, which bypasses the need for stochastic simulations, is to consider the time correlation function of a given state, which can be obtained directly at the level of the master equation for the probabilities d dt p = Wp. [55] In Fig. 6(c), we show the correlation function C (0,0) C (t), i.e. the probability of finding the system in state (0, 0) C (which corresponds to the bottom left corner of the system, see Fig. 1(c)) at time t given that the system was in that state at time zero, again for several values of γ ex /γ in . The correlation function shows a characteristic behavior, with damped oscillations that ultimately asymptote to the steady-state probability of state (0, 0) C . At long times, the decay time of the exponential damping as well as the oscillation period are governed by the first non-trivial eigenvalue of the transition matrix W, that is, the eigenvalue λ = −X R ± iX I with smallest non-zero X R , which corresponds to the slowestdecaying mode in the system. The real part X R corresponds to the inverse of the decay time, whereas the oscillation period T is related to the imaginary part X I through T = 2π/X I [55,57]. The coherence of the oscillations can thus be measured as the number of oscillations that fit within a decay time, and is directly proportional to the ratio R ≡ X I /X R . The oscillation period and coherence obtained in this way are plotted as the black dashed lines in Fig. 6(d,e).
Interestingly, while in the topological regime γ ex /γ in 1 we find that both measures for the oscillation period coincide, the two diverge for γ ex /γ in 1, see Fig. 6(d), once again pointing to a qualitative change in the nature of the oscillations in the system between the topological and trivial regimes. Indeed, at low γ ex /γ in , the inter-nal cycle at the center of the bulk of the system starts to contribute significantly to the current J | , thus lowering the average period T = J −1 | , which in any case becomes an uninformative quantity given its large standard deviation, see Fig. 6(b). On the other hand, the two coherence measures N and R show very similar behavior as a function of γ ex /γ in , see Fig. 6(e). Note that the two differ roughly by a proportionality constant: in the limit γ ex /γ in 1, we can again exploit the fact that our system becomes a unicyclic network with 4(N + 1) states. For a unicyclic network with S states, it is known [55] that the two coherence measures are related through N /R = S/ cot(π/S), which for S = 4(N + 1) = 28 results in the proportionality constant N /R ≈ 3.15 used in Fig. 6(e).
So far, we have focused on the fully chiral case with c = 1. In Fig. 6(f), we show the coherence in the full (c, r) parameter space, again calculated from the first non-trivial eigenvalue of the transition matrix as R ≡ X I /X R . As expected, oscillations are most coherent in the limit of full chirality c → 1 (or c → 0) and large γ ex /γ in , i.e. r → 1. Conversely, coherence vanishes in the achiral case c = 0.5 and in the limit r → 0. In the limit r → 1 (γ ex /γ in 1), we can once more use the equivalence to a unicyclic network to calculate the oscillation period as T = 4(N + 1)/|γ in − γ in | and the coherence as N = 4(N + 1)|γ in − γ in |/(γ in + γ in ) [55]. This shows that the maximal coherence obtainable, which occurs in the fully chiral case, is N = 4(N + 1), and is thus purely limited by the number of states in the system. This results in the maximum value of 28 obtained for the coherence in Figs. 6(e,f).
Returning to the comparison between our topological model for the KaiABC clock and the usual concerted model [8][9][10], it is worth noting that, by including the intermediate states that are disallowed in the concerted model, we raise the upper bound that can be achieved for the coherence of oscillations in the clock. Indeed, the concerted model has 7 × 2 = 14 states and thus an upper bound of N = 14 for its coherence for the case of a unidirectional cycle [55]. Instead, our topological model has 7 × 7 × 4 = 196 states, but in the optimal regime becomes equivalent to a unidirectional cycle with 4 × 7 = 28 states, leading to an upper bound N = 28 for its coherence.

B. Stochastic low Reynolds number swimmers
In another example for a broad class of systems for which the emergence of chiral edge currents in a twodimensional configuration space is desirable, the global cycles we predict would constitute an ideal driving mechanism for synthetic nanoscale swimmers [18,19]. At such small scales, hydrodynamic flows occur at low Reynolds number, and the "scallop theorem" implies that swimming is only possible if the swimming strokes (shape deformations) exerted by the swimmer are nonreciprocal,  7. A topologically-protected three-sphere swimmer, implemented within our 4-state model. In this case, the coordinates x and y describe the extension of the left and right arms, respectively. The internal states A, B, C, D represent whether the left or right arm is primed for increasing or decreasing its length. The emergent topologically-protected edge current results in cyclic shape deformations that enclose the maximally-available area in shape space (in this example, we set the minimal and maximal extension of the arms to be 1 and 4, respectively). thus breaking time-reversal symmetry [58]. In particular, the requirement of nonreciprocity implies that the cyclic shape deformations undergone by a low Reynolds number swimmer must be described by at least two parameters, and that these degrees of freedom must describe cycles that enclose a finite area in the (at least two dimensional) shape parameter space. The swimming capabilities of the swimmer directly increase with the area of said cycle.
The topologically-protected edge currents that we predict naturally describe such cycles in a two-dimensional space, and in fact always describe the largest possible cycle within the available parameter space. As a minimal example, in Fig. 7 we show how to implement a topologically-protected version of the three-sphere swimmer [59] within our 4-state model. The parameters (x, y) correspond to the length of the left and right arms of the swimmer, and we take the minimal and maximal elongation achievable by each arm to be 1 and 4, respectively. Even if all the possible intermediate (x, y) states are accessible, the swimmer spontaneously cycles around the edge states, thus maximizing the area enclosed in parameter space. This cycle emerges without a need for each arm segment to be controlled in the right sequence, reducing the number of external inputs required and/or the internal complexity of the device.
The topological protection and emergent nature of the edge states makes the swimmer robust under external perturbations and constraints. For example, even if an internal malfunction, or the presence of an external obstacle, prevented one of the arms from reaching its full length, the system would still spontaneously describe maximally large cycles in the reduced available parameter space. Moreover, even if an external influence or large fluctuation drove the system into a bulk state, which in this context would imply making both arms "out of step" with each other, the system will diffuse through the bulk of the state space until it reaches an edge state and resumes its global shape cycles. Such occasional unbinding from the edge states will result in run-and-tumble-like behavior for the swimmer.

C. Asymmetric rates and dynamic instability
The inclusion of further features inspired by biology into these models reveals striking observations and directions for future research. For instance, we have until now considered systems characterized by identical rates for the transitions between different states, i.e. with a single value of γ ex and γ in for all the external and internal transitions, respectively. This symmetry need not exist in general, and indeed, in real systems we expect that the transition rates between the different states will be different from each other. We introduce superindices to denote the transition rates between two specific states such that, for example, γ BC ex is the rate of the external transition from B to C, and γ CB in of the internal transition from C to B. In general, there are thus 8 transition rates in the fully chiral 4-state model and 6 transition rates in the fully chiral 3-state model. Robust edge currents will survive as long as the external transitions are significantly faster than the internal transitions with which they compete.
An interesting consequence of having asymmetric transition rates is that they affect the shape of the system oscillations over time. In particular, the typical timescale for moving along the edges is governed by the slower internal transition rates γ in , which constitute the bottleneck. As an example, in Fig. 8(a,b) we show how oscillations in x change in the 3-state model when we increase the rate for the upwards internal transition γ BA in such that γ BA in γ AC in = γ CB in while keeping γ ex γ BA in . The apparent "waiting times" for which the number of subunits x remains constant (vertical edge) are strongly reduced, and we obtain a system for which growth appears to be immediately followed by shrinkage. Moreover, the enhanced upwards internal transition leads to more frequent unbinding from the bottom edge, resulting in a stochastic switching to shrinkage before the right corner x = 100 has been reached. When this 3-state model describes addition of GTP-bound monomers, conversion to GDP-bound monomers, and removal of GDP-bound monomers, this behavior is reminiscent of the dynamic instability of microtubules. Notably, our model is capable of demonstrating phases of growth followed by sporadic phases of shrinkage [11,12] using just three main timescales.
Previous work has suggested that topological phonon We observe "waiting times" between growth and shrinkage. (b) When the upwards internal transition rate is increased to γ BA in = 10 −2 γex, making the system asymmetric, waiting times between growth and shrinkage become negligible, and it is more likely for the system to stochastically unbind from the edge state during growth, switching to shrinkage. (c) Trajectories for two systems, each described by the 3-state model, coupled through the constraint x1 + x2 ≤ Nx describing competition for the same pool of monomers. Here, both systems have identical, symmetric transition rates γin = 10 −4 γex and size Nx = 10. (d) Probability distribution of finding a given x1 and x2 simultaneously, obtained from the same simulation, using a total of 10 7 steps starting from a random initial state. We find that symmetric systems clearly show anti-phase synchronization, see also Movie 3 in the Supplemental Material [49]. modes may play a role in the dynamic instability of microtubules, describing vibrational modes that concentrate near the microtubule cap [60]. The topological states described therein relate to the elastic properties of a microtubule, and are independent of the polymerization or depolymerization dynamics. We propose, on the other hand, topological edge currents in the state space of a Markov model for the dynamics of microtubule polymerization and depolymerization. Hence, our models constitute new candidates for direct comparison with experimental data of microtubule lengths during catastrophe.
Intriguingly, the internal states (A, B, C) in our model could be related to distinct topological states in the elastic properties of the microtubule. For example, our state C, in which the microtubule is primed for depolymerization, could correspond to the topological state described previously [60] for which vibrational edge modes tend to open the microtubule end cap and induce depolymerization, whereas states A and B would correspond to states without protected vibrational edge modes, in which case the end cap remains intact and polymerization proceeds.

D. Coupled systems and synchronization
In contrast with quantum topological systems, in which the boundaries represent real-space edges of a twodimensional material and are thus fixed, the boundaries in stochastic systems represent constraints in configuration space, for example determined by the availability of subunits of a certain type in solution. This implies that the boundaries can dynamically change in time. In particular, if we have two systems (1 and 2), which are determined by their states (x 1 , y 1 ) and (x 2 , y 2 ), a global constraint on the number of subunits of type X would result in the constraint x 1 + x 2 ≤ N x . The boundaries for one system then depend on the state of the other system, i.e. we have 0 ≤ x 1 ≤ N x − x 2 for system 1 and 0 ≤ x 2 ≤ N x − x 1 for system 2.
Remarkably, this boundary coupling can lead to synchronization (or entrainment) between the two systems. Stochastic simulations for the symmetric 3-state model, with constraints y 1 ≤ x 1 and y 2 ≤ x 2 for the second coordinate, show anti-phase synchronization between the x-coordinates of the two systems, see Fig. 8(c,d). This coupled behavior emerges simply from shared physical constraints between the systems, without the need for direct interaction between them. Situations in which two or more systems compete for a limited pool of constituents arise naturally in the biological context e.g. when multiple protein filaments are simultaneously polymerizing. Competition for a limited number of KaiA molecules is also hypothesized to underlie the synchronization of multiple KaiABC oscillators, through the phenomenon of differential affinity [8].
Our results serve as a first demonstration of synchronization due to shared boundaries in topological systems. Nevertheless, much remains to be explored about the nature and mechanisms of synchronization in these systems. We anticipate that, by further studying the parameter space of these models, including asymmetric and weakly chiral systems, different kinds of synchronization (e.g. in-phase synchronization) may be found. As a preliminary result, we find weak in-phase synchronization for an asymmetric system (see Appendix E and Fig. 13 there). Coupling of n > 2 subsystems (through a constraint n i x i ≤ N x ) is straightforward, and it will be interesting to explore the behavior of such systems in the large n limit. More complex constraints may also be imposed, such as nonlinear constraints or constraints that mix the x and y coordinates of different subsystems. While here we have only considered coupling between two identical subsystems, studies of coupled non-identical systems with different transition rates will shed light on whether shared boundaries are capable of synchronizing not just the phases but also the oscillation frequencies of the different subsystems.

V. COMPARISON WITH PREVIOUS PROPOSALS
We now put our results in the context of recent developments regarding topologically-protected states in related systems, namely classical systems with discrete states on a lattice [14,15,41,42]. These pertain to two main classes: (i) stochastic systems with linear dynamics described by a master equation, which includes Refs. 14, 15 as well as the present work, and (ii) nonlinear dynamical systems of the Lotka-Volterra type, which includes Refs. 41, 42. With regards to stochastic systems, previous proposals [14,15] only demonstrated topological states in onedimensional lattices. In particular, these states showed static localization of probability at the system's edge, or at the boundary between two systems with distinct topological properties. Due to their low dimensionality, these systems obviously cannot support edge currents. Our work is thus the first demonstration of the emergence of topologically-protected chiral edge currents in such systems, which are of the highest relevance to biomolecular and biochemical processes, given that these are most often described by a stochastic master equation formalism.
On the theoretical front, these previous proposals [14,15] mapped the stochastic transition rate matrix (or more precisely, the part of it that corresponds to the generator of the probability flux along the periodic direction of the system, after separating it from the divergence operator) onto a corresponding Hermitian Hamiltonian, using a construction similar to that previously used for mechanical lattices [28]. These mappings can change the resulting lattice structure, e.g. Ref. 15 shows a resulting Hamiltonian with higher-order couplings between next nearest neighbors, which are absent in the original transition matrix. The physical relation between the Hamiltonian and transition matrix, or indeed, any general connection they might have, is not entirely transparent. Instead, our work offers a general mapping of any transition matrix into a corresponding non-Hermitian Hamiltonian using the similarities between the master equation and the Schroedinger equation (see Appendix C). Our identification of the transition matrix without its diagonal part (which we show to be irrelevant in a periodic setting) as a typical electronic model allows access to the wealth of previous work in condensed matter theory. Here, our transition matrix retains the same symmetries and connections between neighbors as the electronic model, providing a clear physical relation between them, and a pathway to further progress.
Even more recently, concepts of topological protection have been applied to nonlinear dynamical systems of the Lotka-Volterra type [41,42], representing massconserving population dynamics. In these systems, multiple three-state loops representing predator-prey interactions of the rock-paper-scissors type are linked together to form a lattice. Similarly to the case of stochastic systems, for one-dimensional lattices [41] edge localization has been reported. On the other hand, chiral edge modes have been shown when the rock-paper-scissors loops are assembled into a two-dimensional Kagome lattice [42].
We note, however, some key differences (beyond the fundamental differences between a linear-stochastic and a nonlinear-deterministic system) between the chiral edge modes reported in Ref. 42 and the chiral edge currents reported here. First, the model described in Ref. 42 shows a uniform stationary state, and edge modes are only apparent in the dynamics of perturbations about this stationary state. Our model, on the other hand, can display arbitrarily strong localization at the edges of the system, even in stationary state. A second subtle, but significant difference regards the analogy to the quantum Hall effect. In Ref. 42, it is claimed that the individual rockpaper-scissors cycles behave like cyclotron orbits and give rise to the edge modes. However, one may note that, in their system, the chirality of the rock-paper-scissors cycles matches the chirality of the edge modes, so that e.g. clockwise cycles result in clockwise edge modes. This is in contrast to the actual quantum Hall effect, in which clockwise cyclotron orbits give rise to counter-clockwise skipping orbits at the edges, see Fig. 1(a). This behavior is preserved in our system, in which the chirality of external cycles (analogous to cyclotron orbits) is opposite to the chirality of the edge currents.
Interestingly, both our 4-state and 3-state lattices [see Fig. 1(c,e)], including the different strengths of the internal and external transition rates γ in and γ ex , may be implemented within the same anti-symmetric Lotka-Volterra formalism used in Ref. 41. Future work may thus explore how the edge currents we find here are realized in contexts relevant to population dynamics and game theory.

VI. CONCLUSION AND OUTLOOK
In summary, we show how stochastic systems with out-of-equilibrium cycles at the microscopic scale support chiral edge currents and thus global cycles at the macroscopic scale, along the boundaries of the system's configuration space. The emergence of edge currents can be understood as a topological transition showing unique non-Hermitian properties, which highlight the nonequilibrium character of the system. The flexibility of the underlying configuration space affords exotic features not typically seen in topological systems, such as asymmetric transitions leading to dynamical instabilities, or coupling between subsystems through shared boundaries, leading to synchronization.
The models we propose here are not only interesting due to their application to biochemical systems, but also introduce novel topological phases. In this regard, we note that they are qualitatively different from previous extensions of the original 1d Hermitian SSH model. For instance, our model contains propagating edge currents, whereas other extensions such as the 2d Hermitian [27], 1d non-Hermitian [23,26,33] or 2d non-Hermitian where non-Hermiticity comes directly from complex terms [34], only contain edge localization. We note that this is also the case with regards to the recent attempts to identify correspondence between stochastic systems and topological phases: such 1d models describe stationary localization without global currents [14,15].
Our models exhibit rich phenomenology to be compared with various biochemical processes and other stochastic systems. In contrast to usual models for oscillatory processes at the biomolecular scale, which hardcode the desired global cycle into the structure of the reaction network, our models are constructed from identical repetitive motifs representing simple reactions at the level of single constituents. The global cycle arises as an emergent property of the system, and can occur over widely varying time scales and length scales simply by changing the number of constituents in the system, directly linking structure to emergent function.
In future work, it will be important to investigate and catalogue different oscillatory biochemical processes [61], to elucidate whether topologically-protected edge currents have already been put to use in biological systems. This may prove a difficult task, however, as in practice a system deep in the topological phase (γ ex /γ in 1) behaves very similarly to a unicyclic oscillator that only runs along the edge states. A promising sign that a system may be exploiting topologically-protected edge currents will be the observation of occasional excursions into bulk states, and in particular the observation during these excursions of microscopic "futile" cycles that leave the system unchanged (e.g. repeated single polymerization-depolymerization or phosphorylationdephosphorylation steps). In the example of the Ka-iABC system, by excursions into the bulk we mean the observation of states in which any intermediate number of monomers (i.e. between 1 and 5) have undergone a conformational change and been phosphorylated, see Fig. 6(a). For our fully chiral model of the KaiABC system, our results from Sec. II B imply that excursions into the bulk will be observed every ≈ 1 24 γex γin cycles on average, which highlights again how observation of such excursions becomes increasingly difficult when the system is deeper in the topological phase. In this respect, it would be particularly useful if the ratio γ ex /γ in could be adjusted experimentally within a system of interest, which would be possible e.g. if the external and internal transitions have distinct dependencies on a control parameter (such as concentration of ATP or GTP, pH, salt, or specific molecules known to affect function such as microtubule-targeting drugs [62]). Observation of a system that transitions from incoherent oscillations with frequent excursions into the bulk, to coherent oscillations with rare or no visits to the bulk, as a function of such a control parameter would constitute a strong indication of topologically-protected edge currents at play.
Lastly, our work opens up new directions for the theoretical exploration of new topological phases with in-teresting properties. For instance, our introduction of biologically-inspired features such as dynamical boundaries, new geometries and coupled systems suggest avenues for the realization of novel states. The consideration of heterogeneous transition rates or systems with disorder likewise merits further investigation. More broadly, the relation between general non-Hermitian features like exceptional points and vorticity, and their physical consequences in various systems, forms an emerging area of research. After all, our models suggest new topological states that can implemented in other, more easily tunable systems such as photonics [31][32][33][34] and topoelectronics [37,38]. By producing a blueprint for the emergence of protected edge currents, our formalism can be used to guide the design of new pathways in synthetic biology and materials self-assembly, e.g. towards desired oscillatory assembly or reaction processes. These rapid and continuing developments hold promise for the prediction of new states of matter in both classical and quantum systems.

ACKNOWLEDGMENTS
We thank James Sethna and Benoît Mahault for helpful discussions.
Talks and interactions during the KITP Program Symmetry, Thermodynamics and Topology in Active Matter were also stimulating and thought-provoking. This study was supported by the Max Planck Society.
Appendix A: Stationary state in the fully chiral case We first directly analyze the steady state of the system in the fully chiral case, with γ in = γ ex = 0. We take the bottom edge of the lattice, be it square or Kagome, as an example without loss of generality. We call the probability for a B site on the edge p B , and for the C site p C . The probability of bulk sites away from the edge is called p b . The stationarity condition for the B site reads γ in p C − γ in p B − γ ex p B = 0 whereas for the C site it reads γ in p b + γ ex p B − γ in p C = 0. Using both to solve for p B and p C , we obtain p B = p b and p C = γin+γex γin p b . The fact that p B = p b ensures that the bulk site contiguous to the edge site B is stationary as well, with probability p b , as are all other bulk sites.
Stationarity further implies that sites on convex corners (which have an incoming and an outgoing internal transition) have probability p C as well, see Fig. 2(d) and Fig. 3(d). On the other hand, sites on concave corners (which have an incoming and an outgoing external transition, and appear on non-rectangular 4-state and nontriangular 3-state geometries) have probability p b , see Fig. 2(e) and Fig. 9(b). Note, however, that the latter is true for the 3-state (Kagome) model only in the absence of corner-cutting transitions, see Fig. 9. If such (a) (b) FIG. 9. Edge currents and the resulting cycles are robust with respect to the shape of the boundaries, also in the case of the 3-state (Kagome) lattice. However, in contrast to the 4-state model, in this case one must make sure that corner-cutting transitions across concave corners of the lattice are absent, as otherwise such corners scatter probability into the bulk. In (a), corner-cutting transitions are present (red arrows), and the three concave corners scatter probability into the bulk, making the steady-state probabilities sensitive to the shape of the system. In (b), corner-cutting transitions are absent, the edge current does not scatter, and the steady-state probabilities are independent of the shape of the system. The steady states have been calculated from direct solution of the master equation, for a fully-chiral system with γin = 10 −3 γex.
transitions are present, the edge probability scatters into the bulk, the stationarity condition can no longer be determined using local balance only, and the stationary probability distribution becomes dependent on the system shape and size (and no longer given by just a combination of two values p b and p C = γin+γex γin p b ). The probability current along the edge can be calculated as J edge = γ in p C − γ in p b = γ ex p b . To obtain the global probability of being at the edge in a square 4state model with N x = N y = N , we note that there are n e,C = 4(N + 1) sites with probability p C on the edge, n e,b = 4N sites with probability p b on the edge, and n b = 4N 2 bulk sites, all with probability p b . Normalizing the total probability in the system by setting p C n e,C + p b (n e,b + n b ) = 1, we can calculate p b as The global probability of being at the edge is then calculated as P edge = p C n e,C + p b n e,b = 1 − p b n b or, explicitly, P edge = γin+γex γin (N + 1) + N N + γin+γex γin (N + 1) (A2) which in the limit γ ex γ in results in the expression quoted in Section II B.
To calculate the current J | , i.e. the steady-state current going through the dashed line in the inset of Fig. 6, we note that it can be written as J | = γ in p C = (γ ex +γ in )p b = J edge +γ in p b , where the second term represents the contribution from the internal loop at the center of the lattice. Introducing the explicit expression for p b , we obtain the result in Eq. (4) in the main text.
Appendix B: Symmetries and band structure of W The properties of the transition matrix W governing the Master equation can be analyzed using the decomposition W = A − D, where A ij = i|j is the transition rate from state p j to p i and D ij = δ ij k k|i [50]. A and W are also the adjacency and Laplacian matrices, respectively.
For the 4-state model in Fig. 1(c), we can write the adjacency matrix A explicitly, where we denote the four internal states of cell (x, y) as |(x, y) A to |(x, y) D , as The edge transition probability γ 0 ex interpolates between periodic boundary conditions (PBC) and open boundary conditions (OBC). PBC occur when γ 0 ex = γ ex , while OBC occur when γ 0 ex = 0. The diagonal matrix D can then be easily calculated using its definition above.
This describes a 2d non-Hermitian version of the SSH model [20], which we analyze here for simplicity (the more general case will follow). For a system with PBC, this transition matrix can be expressed in Fourier space as where γ tot = γ ex + γ in and k = (k x , k y ), the reciprocal lattice vector. W k obeys inversion and time-reversal symmetries, IW k I −1 = W −k and W k = W * −k respectively. I is a unitary operator which can be represented as I = σ x ⊗ 1, where σ x is a Pauli matrix and 1 is the identity matrix. In addition, using W k = A k −γ tot 1, we note that A k further obeys sublattice symmetry SA k S −1 = −A k , where S can be represented as S = 1 ⊗ σ z .
With regards to the spectrum, these symmetries suggest that the eigenvalues are either real or complex with conjugate pairs. We indeed see this when analyzing the spectrum of W k : where a(k) = γ in γ ex (cos k x + cos k y ) and β = (γ 2 in − γ 2 ex ). Each pair is illustrated in yellow and green respectively in Fig. 10. The spectrum is even about γ in = γ ex , where bandgaps open for γ in = γ ex . At γ in = γ ex , the bandgap closes to yield degenerate solutions at E(k) = −γ tot .
We can similarly obtain the spectrum for the 3-state system [ Fig. 1(e)] which obeys the expression Upon generalizing the phase space to include the reverse transitions, the spectrum of W k still holds a similar form as Eq. (B2) and is given in the main text. The symmetries of W k and A k that were previously discussed also remain in the general case. SA k S −1 = −A k , which has been shown to quantize the integral of the complex Berry connection Q c m (k) across reciprocal space [26,33]. When the integral is 0, the system is in the trivial phase and will mostly remain in the bulk. When the integral is π, the system is in the topological Zak phase and edge states will dominate. Since W inherits this Berry connection quantization, it therefore supports the Zak phase and long-lived edge dynamics.
W also has inversion symmetry, which similarly quantizes the integrated Berry connection in the Hermitian limit [21,27]. This limit has been well-studied, where previous results show that A (the 2d SSH model [27]) and therefore W exhibit the Zak phase above r * = 1/2. We verify that the bandgap does not close as we interpolate from c = 1/2 into the fully non-Hermitian limit c = 1 above r * (c), using our expression for the band structure of W given in the main text. As the system maintains sublattice symmetry throughout and the bands remain separated above r * (c), it exhibits the Zak phase and edge states.
Appendix D: Stationary state of W in the general case We now analyze the steady state of the general system, including the reverse transitions γ ex and γ in . We focus on the 4-state system, and consider a ribbon periodic along the vertical dimension but open along the horizontal direction (Fig. 11). Defining the vectors P − n ≡ [p n D p n C ] T and P + n ≡ [p n A p n B ] T , the stationarity conditions can be written as where we have defined the matrices Plugging in the equation for P − n into the one for P + n , we obtain a transfer matrix M in the rightwards direction for the probability of the 4-site cells P n ≡ [p n D p n C p n A p n B ] T , that is, the 4 × 4 matrix M that gives P n = M P n−1 (D7) and has the form Two of the eigenvalues of M are always equal to 1, and have identical associated eigenvectors V 1 = [1 1 1 1] T . This reflects that the steady state is uniform in the bulk. However, we also find two other eigenvalues α and 1/α, with 0 ≤ α < 1. The corresponding eigenvectors, V α and V 1/α , are non-trivial, and they are related to each other by i.e., V α is identical to V 1/α , except for a left-right, updown reflection (parity symmetry). These properties strongly suggest that they correspond to the perturbations to the bulk behavior induced by the presence of the left edge (V α ) and right edge (V 1/α ) of the system. The perturbation decays geometrically, with rate α, as we move away from the edge. To ensure that these perturbations indeed correspond to stationary solutions at the edges, we try a solution of the form P 0 = p b (V 1 + ξV α ) at the left edge, where p b corresponds to the probability in the bulk, far away from the edge. We find that the stationarity conditions at the left edge are satisfied if ξ = γ ex − γ ex [(γ in + γ in + γ ex ) − (γ in + γ ex ) − γ in 0] · V α (D10) The excess (or lack) of probability at the edge is therefore δP 0 = p b ξV α , while for the n-th cell away from the edge it is δP n = p b ξV α α n . Assuming that the system size N is large enough such that the probability disturbance decays away from the boundary, i.e. α N 1 or N The black curves represent the boundary r * (c) separating the topological and trivial phases of the system. −1/ log α, the total probability disturbance due to the presence of the edge can be calculated as which is positive if probability accumulates at the boundary, and negative if probability is depleted at the boundary. The total probability flux along the edge can be directly calculated from the steady state probabilities as and is positive for net counter-clockwise edge flux (net flux downwards at the left edge) and negative for net clockwise edge flux (net flux upwards at the left edge).
In the limit of a fully chiral system, we find α = 0 and we recover the results obtained in Appendix A. It is also interesting to note that, according to Eq. (D10), the effect of the boundaries completely vanishes (both in terms of probability disturbance and probability flux) when γ ex = γ ex . Thus, chirality in the external transitions is essential to obtain boundary effects at steady state.
In this way, we can characterize the steady states of the system by simply studying the eigenvalues and eigenvec-tors of a 4 × 4 matrix. Notably, the results are independent of the system size or the shape of the boundaries, even if they give us information about probability accumulation and fluxes at the edges, see Fig. 12.
Appendix E: Synchronization in systems with asymmetric rates Due to the high dimensionality of the parameter space in asymmetric systems (8 transition rates in the fully chi- Probability distribution of finding a given x1 and x2 simultaneously for the same simulation, using a total of 10 7 steps starting from a random initial state. We find that the two systems show weak in-phase synchronization, particularly at initial growth. The system size is Nx = 10. ral 4-state model and 6 transition rates in the fully chiral 3-state model, which double when we consider weakly chiral models), a systematic exploration of the possible types of synchronization in these models is a difficult task which we do not attempt in this work. As a proof of concept, however, we have simulated an asymmetric 3-state model, with symmetric external transitions but with internal transition rates fastest along the vertical direction, slower along the diagonal direction, and slowest along the horizontal direction (γ BA in > γ AC in > γ CB in ), see Fig. 13. We find weak in-phase synchronization, in particular during the initial stages of growth, i.e. at low x.