Topological localization in out-of-equilibrium dissipative systems

TL;DR

This paper demonstrates that topological protection concepts can be applied to out-of-equilibrium dissipative systems, including stochastic networks and driven particle collections, with implications for biological and synthetic systems.

Contribution

It extends topological protection notions from electronic and mechanical systems to far-from-equilibrium dissipative processes, showing their general applicability.

Findings

Topological protection applies to driven, dissipative systems.

Dissipative couplings breaking time reversal symmetry are crucial.

Boundary modes exhibit topological robustness in out-of-equilibrium systems.

Abstract

In this paper we report that notions of topological protection can be applied to stationary configurations that are driven far from equilibrium by active, dissipative processes. We show this for physically two disparate cases : stochastic networks governed by microscopic single particle dynamics as well as collections of driven, interacting particles described by coarse-grained hydrodynamic theory. In both cases, the presence of dissipative couplings to the environment that break time reversal symmetry are crucial to ensuring topologically protection. These examples constitute proof of principle that notions of topological protection, established in the context of electronic and mechanical systems, do indeed extend generically to processes that operate out of equilibrium. Such topologically robust boundary modes have implications for both biological and synthetic systems.

Equations

$$\frac{\partial \mathbf{P}}{\partial t}=\mathbf{WP},$$

$$\underset{ϵ\to 0}{\lim }\mathrm{Tr}\frac{ϵ\rho }{\mathbf{W}+ϵ}=Tr[\rho |u⟩⟨1|],$$

$$\frac{\partial \mathbf{P}}{\partial t}$$

$$\mathbf{J}$$

$$\mathbf{H}=\left(\begin{array}{cc}0& {\mathbf{W}}_1\\ {\mathbf{W}}_1^{\mathrm{T}}& 0\end{array}\right),$$

$$\underset{ϵ\to 0}{\lim }\mathrm{Tr}\left[ϵ{\sigma }_z\rho \frac{1}{\mathbf{H}+ϵ}\right]=\delta n,$$

$${\mathbf{W}}^{\mathrm{b}}={\mathbf{W}}_0^{\mathrm{b}}\cdot {\mathbf{W}}_1^{\mathrm{b}}=\left(\begin{array}{cccc}1& -1& & \\ & 1& -1& \\ & & .& .\end{array}\right)\cdot \left(\begin{array}{cccc}-v& & & \\ w& -v& & \\ & .& .& \end{array}\right),$$

$$\mathbf{H}(k)=\left(\begin{array}{cc}0& {\mathbf{W}}_1(k)\\ {\mathbf{W}}_1^T(k)& 0\end{array}\right)=\left(\begin{array}{cc}0& -v+w{e}^{ik}\\ -v+w{e}^{-ik}& 0\end{array}\right).$$

$$\nu =\frac{1}{2\pi \mathrm{i}}{\int }_0^{2\pi }dk\frac{d}{dk}\ln W_1(k).$$

$$\mathbf{W}={\mathbf{W}}_0^x{\mathbf{W}}_1^x+{\mathbf{W}}_0^y{\mathbf{W}}_1^y,$$

$$\mathbf{W}={\mathbf{W}}_0^x({\mathbf{W}}_1^x+{{\mathbf{W}}_0^x}^{-1}{\mathbf{W}}_0^y{\mathbf{W}}_1^y)\equiv {\mathbf{W}}_0^x\widetilde{{\mathbf{W}}_1^x}.$$

$${\bf H}=\left(\begin{array}[]{cc}0&\tilde{\bf W}_{1}^{x}\\ (\tilde{\bf W}_{1}^{x})^{\rm T}&0\end{array}\right).$$

$$\frac{d\mathbf{Y}}{dt}=-\mathbf{M}\mathbf{Y},$$

$$\mathbf{M}$$

$$\mathcal{R}=\int d{\bf x}\big{[}\eta(\nabla{\bf v})^{2}+D_{\Omega}(\nabla\Omega)^{2}+\Gamma(\Omega-\omega)^{2}+\Gamma_{v}v^{2}+\Gamma_{\Omega}\Omega^{2}\big{]},$$

$$\frac{\partial }{\partial t}\left(\begin{array}{c}\Omega \\ \omega \end{array}\right)=D_{\Omega }\left(\begin{array}{cc}{\nabla }^2-{\lambda }_{\Omega }^{-2}& a\\ -b{\nabla }^2& {\nabla }^2-{\lambda }_{\omega }^{-2}\end{array}\right)\left(\begin{array}{c}\Omega \\ \omega \end{array}\right),$$

$${\mathbf{S}}_{\mathbf{s}}$$

$$\nu =\frac{1}{2\pi i}\int dk\frac{d}{dk}\ln (D_s(k)).$$

$$m\ddot{{\mathbf{r}}_i}=-{\Gamma }_v{\mathbf{r}}_i-{\partial }_{r_i}\sum _{i\ne j}V({\mathbf{r}}_j-{\mathbf{r}}_i,{\theta }_i,{\theta }_j)+D_1{\eta }_1(t),$$

$$I\ddot{\theta }=\tau -{\gamma }_{\Omega }{\dot{\theta }}_i-{\partial }_{{\theta }_i}\sum _{i\ne j}\mathcal{V}({\mathbf{r}}_j-{\mathbf{r}}_i,{\theta }_i,{\theta }_j)+D_2{\eta }_2(t),$$

Full text

Theoretical and experimental studies of biophysical mechanisms such as error correction in DNA replication Hopfield (1974), adaptation in molecular motors controlling flagellar dynamics Wang et al. (2017); Tu (2008); Lan et al. (2012), and timing of events in the cell cycle Cao et al. (2015) are beginning to demonstrate the close connection between robust functioning and non-equilibrium forces. Theoretical results Barato and Seifert (2015); Gingrich et al. (2016) have also elucidated the connection between energy dissipation and fluctuations in a wide class of non-equilibrium systems and have demonstrated how dissipation can be used to tune steady states of many body non-equilibrium systems Nguyen and Vaikuntanathan (2016). However, unlike the behavior and characteristics of equilibrium systems, where no energy is dissipated, general principles governing fluctuations about a steady state or the steady state itself in far-from-equilibrium conditions are just being discovered. As such, strategies that allow for the engineering of specific robust steady states in non-equilibrium biological and soft matter systems are highly desirable. In this paper, we take an approach that is motivated by the physics of topological insulators to construct such states in non-equilibrium systems.

The discovery of robust localized edge states in mechanical systems Kane and Lubensky (2013) that resemble those found in topologically non-trivial electronic Hasan and Kane (2010) and photonic Haldane and Raghu (2008) systems has enabled the development of new design principles. For instance, it has been demonstrated that an assembly of coupled gyroscopes can support topologically protected directed modes at their boundary Nash et al. (2015). Such assemblies can function as a robust waveguides. Localized edge modes have also been discovered in mechanical lattices Kane and Lubensky (2013); Paulose et al. (2015). These edge modes can be set up either at the edge of a lattice Kane and Lubensky (2013) or at topological defects in the interior of the lattice Paulose et al. (2015). They are formally zero energy “free” modes but unlike the commonly found long wavelength zero energy modes in the bulk of marginally stable lattices, these topologically protected boundary modes are highly resistant to perturbations due to disorder and environmental fluctuations. Like their electronic and photonic counterparts, topological modes in the above mentioned mechanical systems have been characterized by topological indices such as winding or Chern numbers Ryu and Hatsugai (2002); Chen et al. (2014); Lubensky et al. (2015).

In this paper, we propose that topologically protected modes can be encoded in a variety of biological and soft matter contexts at the cost of energy dissipation. Like in the electronic and mechanical analogs, the bulk of these systems can be characterized by a topological index (a winding number, in the cases we discuss), and the associated topologically protected zero modes we discover are localized at the edge or an interface. The modes are highly robust and insensitive to perturbations. In some specific contexts, the localization indirectly implies an edge current, analogous to that in the quantum Hall effect Hasan and Kane (2010). In all the instances considered in this paper, topological protection requires that the fundamental equations of motion contain dissipative couplings.

We derive our results in two broad and apparently dissimilar contexts. In the first, we consider biochemical networks with connectivity motivated by those of networks commonly encountered in biophysical information processing and control Murugan et al. (2014). We show that the spectrum of the master equation rate matrix can support localized edge modes that are separated from the bulk via a band gap. In the second, we consider a hydrodynamic description of active matter, specifically that of collections of driven rotating particles in a confined geometry Tsai et al. (2005); van Zuiden et al. (2016), that exhibit large-scale flow localized at the boundaries.

The hydrodynamic equations have a structure very different from that of master equations considered in the first example. The two address phenomena at very different length and time scales, one being a coarse-grained phenomenological description and the other a microscopic approach. Nonetheless, they both describe dissipative phenomena characterized by the production of entropy and lack time reversal symmetry. In fact, we find that they can be characterized by a topological index in the same way as introduced in the context of topological mechanics by Kane and Lubensky Kane and Lubensky (2013), that is by mapping to the Su-Schrieffer-Heeger (SSH) model for electrons on a 1D crystal lattice with two sites per unit cell Su et al. (1979) . Our results then elucidate the design principles required for robust steady states in various biophysical and synthetic systems.

I Topological protection in Markov networks

In this section, we consider the out-of-equilibrium statistical dynamics of stochastic processes described by a Markov network Schnakenberg (1976). Such descriptions are routinely used in statistical mechanics as reduced models of chemical and biophysical processes Van Kampen (2011). For equilibrium systems, the steady states and dynamics of fluctuations about it can be described in terms of energy landscapes, but such a simple description is not available for out-of-equilibrium processes. Hence, any insight into the existence and robustness of steady states of systems driven out of equilibrium by energy-consuming processes, such as those involved in biological functions, is potentially valuable.

Here, we focus on the steady state of lower dimensional networks characterized by uniform or nearly uniform transition rates between various mesoscopic states. This is in analogy with periodically ordered lattices in electronic and meta-materials that host localized topological states at their edgesHasan and Kane (2010). We show in the following that the steady states of certain Markov networks can indeed be mapped onto the ground states of SSH like models for 1D periodic systems with topological properties–the “forward” and “backward” transition rates in the master equation play the role of the hopping rates between the two sites of a unit cell in the SSH model Su et al. (1979).

A simple illustration of this is depicted in Fig. 1: the probability flow in a one dimensional ($1$D) Markov chain with constant rate of site-to-site transition rates in the bulk regions on the left and right of an interface separating them. We aim to show –  by mapping the steady state of the 1D Markov chain to the ground state of a Hamiltonian of a topologically non trivial tight binding model (Fig. 1(b))–  that a possibly “disordered” interface connecting two “bulk” regions of the network with different transition rates can host localized topological modes depending on the transition rates in the bulk. Our results follow from the bulk-boundary correspondence inherent in topological systems, where the existence of localized zero modes at an edge can be predicted by studying the properties of the system in the bulk Kane and Lubensky (2013).

The procedure for establishing the mapping between the stochastic process and the Hamiltonian of a topological tight binding model is distinct and more direct than previous work in which we suggested that the properties of certain out-of-equilibrium Markov states can be understood in terms of topological winding numbers Murugan and Vaikuntanathan (2017). Indeed, we explicitly provide forms of tight binding Hamiltonians whose ground states are the steady states of the out-of-equilibrium stochastic processes we consider. While the 1D network in Fig. 1 is fairly trivial, the procedure outlined below can be used to construct effective tight binding Hamiltonians for more complex biophysical networks with many cycles.

We begin by recalling that the dynamics of Markovian systems can be modeled using a master equation

[TABLE]

where the vector, ${\bf P}$, denotes the probability of occupancy of various distinct states (nodes in a Markov network). This is evolved in time by a transition matrix, ${\bf W}$, with elements $W_{ij}$ indicating the rate of transition from state $j$ to state $i$. The zero right eigenvector of the master equation, $\ket{u}$, specifies the unique steady state of the dynamics Van Kampen (2011). We are interested in conditions under which this steady state zero mode is localized at the interface. Formally, this requirement can be expressed as,

[TABLE]

where $\rho$ is a diagonal matrix with elements $\rho_{i}=1$ for nodes $i$, that lie in the interfacial region and $\rho_{i}=0$ otherwise. Here, $\bra{1}$ is the zero left eigenvector of the master equation corresponding to the zero right eigenvector, $\ket{u}$. The trace in Eq. 2 counts the number of zero modes of ${\bf W}$ that are localized at the interfacial region Kane and Lubensky (2013). This is exactly equal to $1$ if the unique steady state solution, $\ket{u}$, is localized at this interface.

We now show that the condition in Eq. 2 can be related to a topological quantity calculated from the master equation in the bulk network. The state to state transition matrix, ${\bf W}$, does not itself possess the symmetries usually associated with topological protection in electronic or mechanical materials Hasan and Kane (2010); Süsstrunk and Huber (2016). The eigenvalue spectrum of the master equation necessarily has one zero eigenvalue with the rest of the eigenvalues being less than zero Van Kampen (2011). Further, ${\bf W}$ is usually non-hermitian and can have complex eigenvalues. In this form, Eq. 1 does not possess any obvious topological properties.

In order to uncover the topological properties of the master equation, we first note that the rate of change of probability can be expressed as Jack and Sollich (2009),

[TABLE]

where the first is a continuity equation that expresses the conservation of probability with ${\bf W_{0}}$ being a discrete representation of the divergence operator, and ${\bf J}$ is a vector of currents across each link in the network. The matrix ${\bf W}_{0}$ depends only on the topology of the network and not on the transition rates. The current vector, ${\bf J}$, can in turn be expressed in terms of the probability vector, ${\bf P}$, through the matrix ${\bf W}_{1}$, which depends on the transition rates in the network.

We are interested in cases when the trace count as described by Eq. 2 predicts the existence of localized zero modes. We will show below with concrete examples that Eq. 2 can be expressed as a topological invariant by using the fact that the information about the interface between two homogeneous bulk subregions is contained in ${\bf W}_{1}$ and not ${\bf W}_{0}$.

As an illustration of these steps, we first consider the $1$D Markov chain in Fig. 1. The dynamics of the random walker can be described by the master equation in Eq. 1. Using the above mentioned decomposition into ${\bf W}_{0}$ and ${\bf W}_{1}$, we will show that the steady state properties of the $1$D random walker map onto those of the well known SSH model Su et al. (1979); Kane and Lubensky (2013). The central argument is that: if ${\bf W}_{1}$ has a zero right eigenvector, this eigenvector is the unique right eigenvector of ${\bf W}$ and hence the steady state accessed by dynamics under ${\bf W}$. The topological properties of the zero eigenstate of ${\bf W}_{1}$ can be inferred by constructing the Hermitian matrix,

[TABLE]

and considering the trace,

[TABLE]

where $\sigma_{z}\equiv\begin{pmatrix}1&0\\ 0&-1\end{pmatrix}$ is the Pauli z matrix. The trace in Eq. 5, denoted by $\delta n$, is the difference in the number of zero modes of ${\bf W}_{1}$ and ${\bf W}_{1}^{\rm T}$, that are localized in the interfacial region. Since ${\bf W}_{1}$ is constrained by the properties of the master equation to have only $1$ zero eigenvector, the trace in Eq. 5 can only take values $\delta n\in\{-1,0,1\}$. The steady state is localized at the interface between the two bulk regions when $\delta n=1$ as shown in Fig. 1a, whereas it is localized at the opposite ends when $\delta n=-1$, as shown in Fig. 1c.

That $\delta n$ in Eq. 2 can be expressed as a difference of integer topological winding numbers characteristic of the left and right bulks, $\nu_{L/R}$, is well-established as an index theorem in the context of electronic and mechanical topological lattices Kane and Lubensky (2013). Systems with $\delta\nu\equiv|\nu_{L}-\nu_{R}|\neq 0$ have zero modes that are topologically protected and resistant to perturbation. These arguments show that the steady state of the master equation ${\bf W}$ is localized and topologically protected whenever ${\bf W}_{1}$ is topologically non trivial.

The topological winding numbers for the 1D random walker can be derived by considering the transition matrix in one of the bulk regions, ${\bf W}^{\rm b}$:

[TABLE]

where $w$, $v$ denote rates of forward and back hopping in the bulk. In Eq. 6, ${\bf W}^{\rm b}_{0}$ is a discrete representation of a gradient operator that encodes the topology of the network while ${\bf W}^{\rm b}_{1}$ depends on the hopping rates. Because of the repeated nature of these matrices corresponding to the periodic symmetry of the Markov chain in real space, they can be diagonalized in a Fourier basis, $\ket{e^{ikx}}$, in the bulk. They can then be compactly represented as ${\bf W}_{0}(k)=1-e^{-ik}$ and ${\bf W}_{1}(k)=-v+we^{ik}$, where the lattice constant is taken to be unity. The winding number of the bulk region can be obtained from this Fourier representation. Specifically, the bulk Hermitian matrix ${\bf H}$ can also be diagonalized in a Fourier basis,

[TABLE]

In the form of Eq. (7), one can see that ${\bf H}$ is isomorphic to the Hamiltonian of the quantum SSH model with hopping rates, $-v$ and $w$. Using the results obtained for the SSH model, the winding number $\nu$ for the bulk phase can be calculated as Kane and Lubensky (2013),

[TABLE]

Thus we have mapped the zero modes of the 1D Markov chain directly to those of a corresponding SSH model with the backward and forward transition rates playing the role of the two hopping parameters in the tight-binding model of the Hamiltonian. The polarization generated in a topologically nontrivial SSH model is simply related to current generated by the bias in the random walk model. In this sense, dissipation in the bulk of the random walk model plays a crucial role in the generation of localized states. In analogy with the SSH model, two connected chains with opposite polarizations of probability flux will naturally lead to a accumulation of probability at the interface as the system approaches steady state.

The Markov state model in Fig. 1 does not possess multiple cycles. Such cycles can allow for feedback at the cost of energy dissipation and are features common to Markov state representations of many out-of-equilibrium biophysical processes Tu (2008). To derive our results in the context of out-of-equilibrium stochastic models relevant for biological processes, we consider the minimal Markov state model shown in Fig. 2a. that was introduced in Ref. Murugan and Vaikuntanathan (2017). This ladder-like Markov network possesses two horizontal rails with transition rates, $l^{U}$,$r^{U}$, $l^{L}$ and $r^{L}$ denoting the leftward and rightward transition rates in the upper and lower rails. There is also an upward and downward transition probability along each vertical rung of the ladder denoted by $u$ and $d$ respectively.

The Markov state model is composed of two translationally invariant bulk like regions with an interface connecting them. Specifically, the rates of transitions in the bulk regions do not depend on the position along the horizontal axis. The rates in the interfacial region interpolate between the two bulk regions. The transition rates in the right bulk region are denoted by the $\,\tilde{}\,$ symbol to distinguish them from those on the left. As discussed in Ref. Murugan and Vaikuntanathan (2017), the spatial connectivity and structure of this Markov state network resembles that of networks routinely used to study adaptation Lan et al. (2012), kinetic proofreading Murugan et al. (2012); Murugan and Vaikuntanathan (2016), and cell signal sensing Mehta and Schwab (2012). These and other Markov state representations of biophysical processes can often be decomposed into bulk like subgraphs stitched together by interfaces as indicated in Fig. 2a. The subgraphs are formed by finite periodic replication of a particular module or motif.

Since we are mainly interested in networks of the form in Fig. 2, which possess translational symmetry along one (horizontal) axis and the interface spans the other (vertical) axis, we decompose rate matrix of this system as

[TABLE]

where ${\bf W}_{0/1}^{x/y}$ are square matrices and are discrete representations of the continuity operator in the (horizontal) $x$ and (vertical) $y$ directions. Since the interface spans the vertical axis, we choose the decomposition,

[TABLE]

These arguments imply that any master equation rate matrix ${\bf W}$ can be factorized as, ${\bf W}={\bf W}_{0}^{x}\tilde{\bf W}_{1}^{x}$. Again, the crucial point in this decomposition is that ${\bf W}_{0}^{x}$ does not depend on the transition rates and possesses no interfaces. Using this property and arguments similar to the case of 1D random walk model, one can show that W has topologically protected modes whenever the following Hermitian operator constructed with $\tilde{\bf W}_{1}^{x}$ is topologically non trivial,

[TABLE]

Like the master equation rate matrix, the Hamiltonian ${\bf H}$ is also composed of two bulk phases connected by interfaces. A schematic of the connections in the effective Hamiltonian is provided in Fig. 2. We provide the detailed derivation and explicit expression for ${\bf H}$ in the Supplementary information.

The polarization implicit in the Hamiltonian ${\bf H}$ reflects the currents generated by the master equation rate matrix ${\bf W}$. For instance, the effective Hamiltonian we have depicted in Fig. 2b, is composed of two horizontal rails each with a lattice structure that resembles that of the SSH model. The two rails can potentially have polarizations with the same magnitude but opposite signs. This choice corresponds to the case $l^{U}=r^{L}$ and $r^{U}=l^{L}$. In such cases, the effective Hamiltonian can have a net polarization only when the vertical links connecting the two rails break symmetry. Indeed, this condition reflects the requirement that the symmetry between the two rails of the ladder in the master equation rate matrix be broken, $u\neq d$, for a current along the horizontal axis of the bulk networks. In this sense, topological protection depends crucially on the currents generated by the dissipative fluxes in the master equation rate matrix. In the context of the models considered here, systems without dissipative fluxes in their bulks cannot support a topologically protected mode. In Fig. 2c, we show the numerically computed count of localized zero eigenmodes of ${\bf W}$ using Eq. 2, and compare it to that of the effective matrix, ${\bf H}$. We also show the gap in the eigenvalue spectrum for ${\bf W}$ and ${\bf H}$ in Figs. 2c & d respectively, by numerically computing the first two eigenvalues as the hopping rates are varied. Such a gap between the steady state and subsequent eigenvalues in the eigenvalue spectrum make the steady state robust against random fluctuations. We indeed find that as long as the Hamiltonian $H$ has a winding number mismatch that supports a localized zero mode at the interface, ${\bf W}$ is also guaranteed to have a zero mode localized to the interface.

This is the main result of the first part of the paper. It establishes that the topological properties of ${\bf W}$ can be computed by simply constructing winding numbers for the matrix, $\tilde{\bf W}_{1}^{x}$. This mapping allows us to infer the properties of a stochastic out-of-equilibrium system in contact with thermal reservoirs in terms of the topological properties of an isolated quantum mechanical system described by the Hamiltonian ${\bf H}$. The Hamiltonian reflects the polarization of the bulk master equations that support a current.

II Topological protection in many body systems

There is broad interest in dissipative, steady-state structures (such as ordered phases and topological defects therein) that emerge in collections of particles driven away from equilibrium in both synthetic Cates and Tailleur (2015) and biological contexts Marchetti et al. (2013); Kruse et al. (2005). Such steady states, if topologically protected, are likely to be robust against disorder and can be categorized into distinct topological classes – thus guiding applications that involve organization away from equilibrium Whitelam et al. (2012); Nguyen and Vaikuntanathan (2016).

The results of the previous sections are, however, specific to effectively single particle Markov state processes. A natural question is whether similar statements about topological modes can be made for collections of many interacting particles out of equilibrium. Indeed, an effective Markov state representation of the dynamics of a many particle interacting system cannot be simplified in the same way as the models in Figs 1 &  2. While a full microscopic description of a many-particle interacting system is in general intractable without doing detailed simulations, the relaxation of such a system towards equilibrium or a steady state can be described in terms of a few slowly decaying collective modes with long wavelengths. Such a hydrodynamic description in terms of long wavelength collective modes differs fundamentally from the electronic or mechanical materials with periodic order for which topological physics is typically demonstrated Hasan and Kane (2010); Kane and Lubensky (2013).

In this section, we show that the equations of motion derived for a set of $m$ hydrodynamic fields, $\{{\bf X_{i}}({\bf x},t)\}$, each corresponding to a conserved quantity, can have steady state solutions that are topologically protected. We focus on purely dissipative hydrodynamics, where the principles of linear irreversible thermodynamics derived in the seminal works of Prigogine, Onsager, deGroot and Mazur de Groot and Mazur (1969) can be used to write the equations of motion as,

[TABLE]

where ${\bf Y}=\{{\bf X_{1}},...,{\bf X_{m}}\}$ is the set of all the $m$ hydrodynamic variables in a given problem, and the matrix ${\bf M}$ contains the dissipative couplings between the various hydrodynamic variables. Since we consider purely dissipative processes, an excited collection of particles should relax monotonically in time towards the steady state which corresponds to minimal or no dissipation. This constrains the eigenvalues of ${\bf M}$, which controls the dynamics of the variables ${\bf X_{i}}({\bf x},t)$, to be positive. This property allows us to decompose a suitably symmetrized form of ${\bf M}$ as a matrix product ${\bf D}{\bf D}^{\rm T}$. We can then extract its topological properties using suitable index theorems for the square root operator given in Ref. Kane and Lubensky (2013): ${\bf S}=\begin{pmatrix}0&{\bf D}\\ {\bf D}^{\rm T}&0\end{pmatrix}$, that possesses a symmetry, $\{\sigma_{z},{\bf S}\}=0$, that leads to positive and negative pairs of eigenvalues, and $S^{2}$ has the same eigenvalue spectrum as the original dynamical operator, ${\bf M}$, including its zeros.

With this background in mind, we now demonstrate topologically protected flow in a simple fluid setup as a prelude to application to more complex chiral, active fluids Tsai et al. (2005); Fürthauer et al. (2013). We consider a thin layer (quasi 2D in the $xy-$plane) of viscous, incompressible, fluid on a solid substrate and bound between two vertical plates(positioned in the $yz-$ plane). Moving one of these bounding plates with velocity $v_{0}$ induces a flow in the $y-$direction. The hydrodynamic equation for the velocity of the fluid is determined by momentum conservation which includes the dissipative processes in both the bulk of the fluid and at its surface in contact with the substrate which acts as a “momentum sink”. In the linear response regime, the dissipation potential governing the flow in this setup is, $\mathcal{R}=\int d{\mathbf{x}}\,(\eta(\nabla\mathbf{v})^{2}+\gamma\mathbf{v}^{2})$, where $\eta$ is the fluid viscosity and $\gamma$, the substrate friction coefficient. The equation of motion obtained by extremizing the dissipation functional, $\mathcal{R}$, is: $\rho Dv/Dt=(\nabla^{2}-\lambda^{-2})v$, where $-(\nabla^{2}-\lambda^{-2})$ is the hydrodynamic operator, ${\bf M}$, for this setup, and $\lambda\equiv\sqrt{\eta/\gamma}$ is a friction length. The resulting steady state velocity profile is exponentially localized near the moving plate: $v_{y}(x)=v_{0}e^{-x/\lambda}$ assuming no slip boundary conditions.

Although, we need appropriate boundary conditions (here, no slipping at the plate surface) to obtain the exact spatial profile of the fluid velocity, topological protection —which we demonstrate in the next paragraph —implies that the fluid velocity is localized near a moving plate near the plate irrespective of the details of the boundary conditions, as long there is a finite dissipative coupling ($\gamma>0$) with the substrate.

The simple linear operator, $\nabla^{2}-\lambda^{-2}$, which occurs generically in dissipative hydrodynamics, can be characterized by a topological index that is related via an index theorem to the localized zero mode of the operator, $i.e.$ the steady state solution of the hydrodynamics. To see this explicitly, we construct the square root, $-(\nabla^{2}-\lambda^{-2})={\rm D}_{0}{\rm D}_{0}^{T}$, of the linear hydrodynamic operator. This may be represented by a finite ($N\times N$) matrix by discretizing the Laplacian operator on a $1D$ mesh in real space with a lattice constant, $h$:

[TABLE]

where the $D_{0}$ matrix is seen to be identical to the SSH Hamiltonian, with v,w=\frac{1}{2}\big{(}\sqrt{\lambda^{-2}+4h^{-2}}\mp\lambda^{-1}\big{)}. As shown in the case of an SSH chain Su et al. (1979), the interface with vacuum (in this case, a hard wall) can host a localized edge mode (in this case, of the velocity at steady state) when the bulk is characterized by a finite topological index, as it happens for $w>v$. This is satisfied in the choice of decomposition above irrespective of the fineness of the mesh, i.e. the value of $h$. The only exception to this is for the case $\lambda^{-1}=0$ i.e. there is no friction from the substrate, when the above construction results in $v=w$, the topologically trivial case corresponding to a closing in the gap of the eigenvalue spectrum of $(\nabla^{2}-\lambda^{-2})$. The details of the calculation and an argument connecting friction to an effective polarization are presented in the SI. We note that while the lattice discretization $h$ was introduced to construct the decomposition in Eq. LABEL:squareroot:1, the topological index is independent of it. Finally, the decomposition of the linear operator is valid even with $w$, $v$ interchanged. While such a choice leads to a trivial topological state, it doesn’t invalidate our conclusions. Indeed, the topological state of the SSH model is similarly changed when the unit cell is shifted by one lattice unit. For the hydrodynamic equations to have a topological mode, we simply require that there exists a decomposition that leads to a non zero topological index.

As a nontrivial application of this notion of topological protection to a many-particle system, we focus on collections of actively rotating particles (rotors) with an intrinsic “spin” angular momentum degree of freedom. The collective dynamics of such particles are more complex and relatively less understood than those that involve linear self-propulsion Cates and Tailleur (2015) though there is now a wide range of experiments where such active angular momentum injection in the bulk can be realized in a controlled manner (see Refs. van Zuiden et al. (2016); Nguyen et al. (2014)). Examples include liquid crystals in a rotating magnetic field Migler and Meyer (1991), shaken chiral grains Tsai et al. (2005) and light-powered colloids Maggi et al. (2015). Instances of naturally occurring flowing matter with actively generated rotation by molecular motors range from the rotational beating of flagella of swimming bacteria to active torque generation in the cellular cytoskeleton Fürthauer et al. (2013). Such molecular torques are potentially biologically significant – and have been implicated in the streaming chiral flows in the actin cortex that lead to left-right symmetry breaking of a developing organism during morphogenesis Naganathan et al. (2014), for example. Boundary modes where the rotors circulate around the edge of a container have been seen in Ref. Tsai et al. (2005) and more recently in simulations in Ref. van Zuiden et al. (2016) but connections to topological protection, if any, have not been explored. We now demonstrate that active spinners in confined geometry can indeed support topologically protected localized edge modes. These modes allow the system to robustly localize flows and transport to boundaries and provide a route to the breaking of chiral symmetry.

The dynamics of a collection of actively spinning rotors can be described by a coarse-grained hydrodynamic theory formulated in terms conservation laws and constitutive relations derived using the general principles of irreversible thermodynamics de Groot and Mazur (1969). The key hydrodynamic variables in the theory are the fields of intrinsic rotation rate (or “spin” angular velocity), $\Omega({\bf x,t})$, and the linear flow velocity of the rotors, ${\bf v}({\bf x},t)$, corresponding to the conservation of angular momentum and linear momentum respectively. We assume incompressibility based on the observation that the density of rotors remains nearly uniform van Zuiden et al. (2016). A crucial feature of these phenomenological equations of motion is the dissipative coupling between the angular velocity, and the vorticity, $\omega=(\nabla\times{\bf v}({\bf x}))_{z}/2$ through the rotational strain rate, $\Omega-\omega$, which is the generalized force corresponding to a thermodynamic flux: the antisymmetric component of shear stress induced by the relative rotation of adjacent fluid elements Fürthauer et al. (2012); Klymko et al. (2017a, b). The antisymmetric stress, is responsible for the dissipative transfer of angular momentum from the intrinsic spin degree of freedom to the linear bulk vorticity Fürthauer et al. (2012); Klymko et al. (2017a, b), and is always present unless the spin and vorticity become equal.

The hydrodynamic equations corresponding to the complex rotor systems can be inferred using the principles of linear irreversible thermodyamics de Groot and Mazur (1969) from a dissipation functional,

[TABLE]

derived from the proportionality of the thermodynamic forces and fluxes through the dissipative phenomenological coefficients: the spin-vorticity coupling, $\Gamma$, the rotational diffusion constant, $D_{\Omega}$, the viscosity of the medium, $\eta$ and the substrate friction coefficients for angular and linear velocity: $\Gamma_{\Omega}$ and $\Gamma_{v}$. These dissipative processes that contribute to the net rate of entropy production in Eq. (14), together with the uniform active torque, $\tau$, driving each rotor, determine the dynamics at steady state.

The coupled dynamical equations for the spin angular velocity (redefined as a difference from its bulk steady state value: $\Omega({\bf x})\rightarrow\Omega({\bf x})-\Omega_{b}$), and vorticity, $\omega=(1/2)(\nabla\times{\bf v})_{z}$, are derived by extremizing the above dissipation functional, $\mathcal{R}$, as:

[TABLE]

where the decay length scales, $\lambda_{\Omega}^{-2}=(\Gamma+\Gamma_{\Omega})/D_{\Omega}$, and, $\lambda_{\omega}^{-2}=\Gamma_{v}/(\eta+\Gamma)$, and coupling coefficients, $a=\Gamma/D_{\Omega}$, $b=\Gamma/(\eta+\Gamma)$, are defined in terms of relevant dissipative parameters Tsai et al. (2005). This is the bulk steady state to which the collection of spinning rotors decays under the action of an active intrinsic torque, $\tau$, and in the presence of substrate friction, $\Gamma_{\Omega}$: $\Omega_{b}=\tau/(\Gamma+\Gamma_{\Omega})$ Tsai et al. (2005); van Zuiden et al. (2016). In a confined geometry, the rotors are prevented from rotating freely at the walls, which induces a spatial profile in the angular velocity, and therefore in the vorticity, at steady state Tsai et al. (2005); van Zuiden et al. (2016).

We now demonstrate that Eq. (15) does indeed lead to topologically protected zero modes corresponding to its steady state solutions. The details of the derivation are given in the SI. While the full matrix ${\bf M}$ describing the rotor dynamics can be decomposed to reveal its topological properties (as we do numerically in Fig.3), the calculation of a topological index simplifies considerably at steady state. Briefly, we integrate out the vorticity field to obtain a linear hydrodynamic operator for spin angular velocity, ${\bf L}_{s}$, at steady state. This can be discretized in real space and decomposed as ${\bf L}_{s}={\bf D}_{s}{\bf D}_{s}^{\rm T}$. The resulting ${\bf D}_{s}$ is a sparse matrix with elements, $({\bf D}_{s})_{i,j}=v\delta_{i,j}+w_{1}\delta_{i-1,j}+w_{2}\delta_{i-2,j}$, that are non-zero only in the diagonal and two nearest off-diagonals. This allows us to define a matrix, ${\bf S}_{s}$, which is represented in Fourier basis (with appropriately defined periodic boundary conditions) as,

[TABLE]

such that ${\bf S_{s}}^{2}$ and ${\bf L}_{s}$ have the same eigenvalues. ${\bf S}_{s}$ corresponds to an SSH like model for 1D topological insulators with both nearest and next nearest neighbor hopping. The topological properties of this model can be characterized in terms of a winding number Kane and Lubensky (2013),

[TABLE]

Although we consider rotors confined in a 2D geometry, we discretize the rotor equations along one direction transverse to the boundary (and therefore, to the edge currents), leading to a winding number Ryu and Hatsugai (2002). We again stress that the winding number is not a function of the lattice discretization.

In the presence of finite coupling between spin and vorticity, the gap in the bulk spectrum of ${\bf L}_{s}$ closes only if $\lambda_{\omega}^{-2}=0$, i.e. if friction with the substrate vanishes, $\Gamma_{v}=0$ (Fig. 3d). This corresponds to $v+w_{2}=w_{1}$, i.e. winding number of zero (see SI for the winding number calculation). In general, in the presence of finite friction, the winding number is $\nu=2$ corresponding to two localized modes. The exact spatial profile of the spin velocity and vorticity is a linear superposition of these independent modes and depends on the specific boundary conditions imposed in a confined geometry Tsai et al. (2005).

In Fig.  3 we plot the eigenvalues of the discretized matrix for a slab geometry. The eigenvalues are real and non-positive. We calculate numerically the trace in Eq. (2) for the full rotor hydrodynamic operator defined in Eq. (15), to show that there are indeed two localized zero modes in the presence of friction and finite spin vorticity coupling. This corresponds to the winding number, $\nu=2$, calculated for the rotors at steady state. Further, we observe that in the presence of friction and finite spin vorticity coupling, there is a gap in the spectrum of the eigenvalues between the zero eigenvalue and subsequent eigenvalues. As discussed in other contexts Kane and Lubensky (2013), such a gap protects the steady state zero eigenvalue solution against perturbations. The gap and the localized modes mirror features encountered during the analysis of the master equation. Like in the master equation analysis, the hydrodynamic equations of motion in Eq. 15 are dissipative, non-Hermitian and their topological properties if any, are not immediately apparent.

As a further test of these ideas, we set up molecular dynamics simulations of a system composed of driven rotors confined in a box. It has been well established that the hydrodynamic equations in Eq. 15 are good descriptors of the dynamics of such rotor systems Tsai et al. (2005); van Zuiden et al. (2016). The equations of motion of individual rotors in our molecular dynamics simulations are:

[TABLE]

where $\Gamma_{v}$ is the substrate friction, $\gamma_{\Omega}$ is the friction imposed on the $\theta$ degree of freedom, $\tau$ is the applied torque, $\eta_{1}(t),\eta_{2}(t)$ are delta function correlated Gaussian random variables, and the interaction potentials $V$ are derived from the Yukawa pair potential van Zuiden et al. (2016). The detailed parameters used for the simulations are described in Methods.

In order to check for the existence of localized vorticity predicted by the theory, we computed the time averaged velocities of the center of mass of the rotors. As is clear from Fig. 4a , this system supports localized edge flows. This flow is immune to backscattering in the presence of obstacles and persists even if the edges of the box used in the simulation are jagged (as depicted in Fig. 4b). Further, when the rotors are actively driven only in a partial region of the simulation box (Fig. 4c), the flow is localized at the boundary of this region of driving.

According to a solution of the hydrodynamic equations with a specific choice of boundary conditions van Zuiden et al. (2016) (SI), one expects the inverse localization length $\kappa$ to scale like $\kappa^{2}\propto\Gamma_{v}$. This solution, and the theoretical analysis above, predicts a loss of topological protection at zero substrate friction. In agreement with these predictions, the boundary flow velocity continues to be robustly localized as substrate friction $\Gamma_{v}$ is reduced until very small values of friction, when the flow profile turns noisy (Fig. 4d). The loss of topological protection in a finite system at non-zero friction can be interpreted as a signature of the narrowing of the associated gap in the eigenvalue spectrum of flow solutions in relation to fluctuations as friction is decreased, resulting eventually in an unstable steady state. Similar features are apparent as topological protection is lost in finite size electronic and mechanical topological insulator models Rocklin et al. (2016). These observations taken together constitute a demonstration of the robustness of edge flows to boundary conditions.

We have shown here that hydrodynamic equations describing the dynamics of rotors can be characterized by a topological index. This quantity is calculated for the bulk of the flow and guarantees the existence of flows localized at the edge of a confined fluid as a consequence of the principle of bulk-boundary correspondence well-known in topological physics Hasan and Kane (2010); Kane and Lubensky (2013). This conclusion is reached without explicitly solving hydrodynamic equations which require one to assume hydrodynamic boundary conditions that are not always intuitive. For example, the direction of the edge flow seen in these active rotors is not obvious a priori. Although the flow velocity obtained from simulations has been fitted to the flow profile corresponding to zero tangential stress from the confining wall van Zuiden et al. (2016), we show in the SI that a no-slip boundary condition at the wall also leads to a localized velocity profile. This shows that the solutions guaranteed by topological considerations are independent of the exact boundary conditions for a given hydrodynamic problem.

Using two disparate examples of complex dynamical systems that are out of equilibrium, we have shown in this paper that topologically protected states can arise in principle in a variety of dissipative systems: both stochastic networks and active flows. Unlike other systems for which topological protection has been previously explored, such as mechanical lattices Kane and Lubensky (2013) or propagating sound modes in hydrodynamic equations Souslov et al. (2016); Shankar et al. (2017), dissipation is key to the phenomena considered here. In both cases considered, the topological properties are not readily apparent from the structure of the relevant operators that govern their dynamics. We reveal their topological properties by decomposing them suitably in order to map the properties of their steady state solutions to the zero energy states of 1D Hamiltonian models that can be characterized by a topological index. Our work indicates that both single and many particle dynamics with interactions at both microscopic and macroscopic scales can result in topological modes that are localized at boundaries. Further, our results also suggest that the chiral edge flows ubiquitously seen in synthetic and biological matter are potentially robust, topologically protected modes. The robustness in this case is explained using broadly applicable ideas from linear irreversible thermodynamics. This has applications for self-assembly of metamaterials as well as guarantee robust localized flows of both information and matter in biology.

III Acknowledgments

We gratefully acknowledge very useful discussions with William Irvine, Sid Nagel and Tom Witten. K D. and S.V. were funded by NSF DMR-MRSEC 1420709. SV also acknowledges funding from the University of Chicago and the Army Research Office under grant number W911NF-16-1-0415. KKM acknowledges support from a National Institutes of Health Grant R01-GM110066. He is also supported by Director, Office of Science, Office of Basic Energy Sciences, Chemical Sciences Division, of the U.S. Department of Energy under contract No. DE- AC02-05CH11231

IV Methods

The molecular dynamics simulations were performed by evolving the Langevin equation with Euler dynamics. The code is available upon request. The rotors in our simulations were composite particles composed of three equally spaced point particles arranged along a line. The particles interact according to a Yukawa potential with a decay length $2d$ where $d$ is the spacing between the particles in the rotor. This force field was used to derive the interatomic forces and torques on the rotor. We set $\Gamma_{v}=\gamma_{\Omega}$ as the substrate friction. The simulations in Fig 4 a, c, d were performed with $N=160$ particles in a two dimensional square box with length $L=80$. The size of the rotors is $d=2.5$. The charge on the rotors (for the Yukawa potential) was set to $q=2.5$. The simulations in Fig 4 b. were performed with $N=362$ particles in a box with size $L=120$. The applied torque $\tau$ in all the simulations was $\tau=10$. The variance of the Gaussian noise in simulations was set using $D_{1}=1$ and $D_{2}=0.5$.