Topological Sound and Flocking on Curved Surfaces

Active systems on curved geometries are ubiquitous in the living world. In the presence of curvature orientationally ordered polar flocks are forced to be inhomogeneous, often requiring the presence of topological defects even in the steady state due to the constraints imposed by the topology of the underlying surface. In the presence of spontaneous flow the system additionally supports long-wavelength propagating sound modes which get gapped by the curvature of the underlying substrate. We analytically compute the steady state profile of an active polar flock on a two-sphere and a catenoid, and show that curvature and active flow together result in symmetry protected topological modes that get localized to special geodesics on the surface (the equator or the neck respectively). These modes are the analogue of edge states in electronic quantum Hall systems and provide unidirectional channels for information transport in the flock, robust against disorder and backscattering.


I. INTRODUCTION
Flocking, defined as the self-organized ordered motion of collections of self-propelled units [1,2], has been the flagship of active matter for some time now [3,4]. Active entities dissipate energy to perform work and generate motion, leading to sustained and local breaking of detailed balance. Examples abound in the living world, ranging from bird flocks [5] to bacterial suspensions [6] and migrating cells [7], and include synthetic analogues, such as reconstituted cytoskeletal extracts [8][9][10], vibrated granular media [11,12] and self-propelled colloids [13]. Active matter as a field aims to provide broad organizational principles applicable to a wide class of these non-equilibrium systems over many scales.
Collective cellular motion on curved surfaces is ubiquitous in developmental processes, such as morphogenesis and embryonic development [14][15][16], or when cells migrate in the gut [17,18] or on the surface of the growing cornea [19], and also affects cell division [20]. Recent in vitro work has demonstrated a direct effect of substrate curvature on cytoskeletal alignment and cell motility in epithelial cells [21]. Understanding the behavior of active matter on curved surfaces or confined by curved boundaries. is therefore timely. There has been growing recent interest in understanding this at a fundamental level, with the focus divided between the effect of curved confining walls on so-called scalar (non-aligning) active matter [22][23][24][25][26] and on aligning active matter systems, of either nematic [9,10,27] or polar [28,29] symmetry. Even at the level of non-interacting self-propelled particles, the curvature of confining walls can yield surprising features, such as inhomogeneous density and pressure profiles [22,25] and the breakdown of an equilbrium interpretation [23,24]. In the presence of aligning interactions that promote orientational order, curvature has an even more profound effect since it frustrates such or-der, often requiring topological defects [30] that in active systems become dynamical and are capable of driving spatio-temporal patterns and complex motion [10]. With flexible walls or membranes present, activity can lead to spontaneous motion and rectification [31,32].
A generic property of the ordered state of polar active matter is spontaneous flow and thus the breaking of time-reversal symmetry [4]. It is known that carefully engineered lattice structures with flows induced either spontaneously by activity [33], or through an external drive [34], can host exotic unidirectional sound modes that are localized at the edges of the sample and topologically protected. The presence of topologically protected edge states in classical phononic [35,36] and photonic [37] systems has lead to extensive exploration of topological meta-materials, with properties akin to electronic topological insulators and quantum Hall states [38]. Here we show that a polar active fluid on a curved substrate supports similar topologically protected modes, even in the absence of any underlying periodicity or lattice structure. This should be contrasted with many of the systems considered previously which required a carefully-designed meta-material structured on an artificial lattice. The phenomenon reported here is akin to the one recently found in geophysical flows of oceans or the earth's atmosphere, where equatorially trapped Kelvin waves were reinterpreted as topologically protected modes [39]. In that case the flow is imposed externally by the earth's rotation. In our active fluid, in contrast, flow occurs with no external drive, resulting in spontaneous topologically protected modes.
The presence of these topological modes relies on three important ingredients: • the spontaneous polar order and associated flow that breaks time reversal symmetry; • the fact that in polar active fluids the order parameter also plays the role of a flow velocity, re-sulting in distinctly non-equilibrium self-advection not present in equilibrium polar fluids [4]; • the non-zero gaussian curvature of the underlying substrate.
We emphasize that the long-wavelength topologically protected modes discussed here are generic, in the sense that they occur for active polar flow on any curved surface of non-vanishing gaussian curvature. Recent work has considered active polar patterns on a cylinder [40]. In this case the gaussian curvature vanishes and there are consequently no topologically protected sound modes. In the following we demonstrate the phenomenon explicitly for flocking on the sphere, which has constant positive curvature, and on the catenoid, which has negative, spatially varying curvature. In Sec. III, we analyze the continuum Toner-Tu model of flocking on a sphere and analytically compute the steady state configuration of the polar ordered flock. Due to the curvature, the ordered state is forced to be inhomogeneous in general. On a sphere polar order additionally requires topological defects or vortices (the "hairy-ball theorem" [41]) in order to satisfy the constraints imposed by the global topology of the surface. With a well motivated ansatz, we show that the covariant hydrodynamic model is capable of predicting generic inhomogeneous steady ordered phases that accommodate the curvature and topology of the underlying surface in a natural fashion. The steady flocking state on the sphere corresponds to the rotating band seen in recent particle simulations by Sknepnek and Henkes [29] and is a novel find in itself as it is peculiar to the active system (a passive polar liquid crystal on a sphere would have a very different equilibrium profile).
Having obtained a steady ordered state, in Sec. III B we examine its excitations. Even though our system is overdamped due to the presence of a substrate, the ordered polar flock supports long-wavelength propagating sound modes [42]. The presence of curvature introduces an additional length scale in the problem and gaps the sound spectrum at long wavelengths. In other words the propagation rate of long-wavelength density fluctuations is finite. This is a distinct property of active polar fluids and it arises because the polarization field plays the dual role of the order parameter and flow velocity and is therefore subject to the same lensing effect that forces flow to move along geodesics on curved surfaces [29,39]. With a spectral gap opened, in Sec. IV we show that a polar ordered flock on a curved surface supports topologically protected sound modes that are localized to special geodesic curves on the surface (at which the gap in the spectrum closes). In Sec. V, we compute the steady state of an ordered polar flock on a catenoid and show that topological modes are also present on this negativecurvature surface. In contrast to the setup of Souslov et al. [33], where the lattice structure was instrumental (along with the active flow) in generating a gapped spectrum at intermediate frequencies, here the curvature does the job though now at long wavelengths, implying that the result is quite general.

II. TONER-TU EQUATIONS ON A CURVED SURFACE
We consider an active polar fluid on a 2d surface. To make generic predictions independent of specific microscopic realizations, we work in the continuum limit and use the well-tested hydrodynamic description of a fluid of overdamped self-propelled particles provided by the Toner-Tu equations [2,43,44], appropriately modified to account for the curvature of the sphere [28]. Mass conservation implies a continuity equation for the density field, ρ, with µ = θ, ϕ and p = ρu the polarization density of the active fluid. Due to the presence of a substrate, momentum is not conserved and the particles velocity is assumed to be aligned with their direction of self-propulsion, leading to the identification of u with the flow velocity of the active fluid. Note that on a curved surface parallel transport of vectors requires the use of covariant derivatives [45], where Γ ν αµ are the appropriate Christoffel symbols. The equation for the polarization density is given by Note that u here plays the dual role of an order parameter field (polarization) and velocity, as discussed in [4]. The transport coefficients ν and ν ′ are the shear and bulk viscosities (or anisotropic elastic constants), respectively, a, b > 0 are coefficients setting the magnitude of the mean field polarized state for ρ > ρ c (the critical density for the flocking transition), λ is a kinematic convective parameter and v 1 > 0 is a compressional modulus. The last term on the right hand side of Eq. (3) is the first term in a density expansion of the gradient of the swim pressure [46,47] that describes the flux of propulsive forces across a unit plane of material. There are other nonlinear terms in the original Toner-Tu equations, but we only retain the most important ones here. In particular, we keep the convective nonlinearity λp · ∇p that is responsible for longranged order in 2d [43,48] and the leading density dependence in the symmetry breaking (aρp) and pressure like terms (v 1 ∇ρ) that lead to dynamical self-regulation [49], phase-separation [50][51][52] and long-wavelength instabilities of the ordered phase [53,54]. The absence of Galilean invariance means that λ = 1/ρ. Additional nonlinear advective terms ∼ λ 2 p∇ · p, λ 3 ∇|p| 2 are also present in general, but do not qualitatively change our results below (see Appendix (A) for an analysis with λ 2 , λ 3 = 0). Curvature enters Eq. (3) in two crucial places (apart from the covariant derivatives): (i) the cubic term setting the magnitude of the polarization explicitly involves the metric tensor g (|p| 2 = g αβ p α p β ) and (ii) the gaussian curvature K G explicitly appears in the viscous term because the strain rate tensor is a symmetrized derivative of the velocity and the covariant derivatives do not commute. The presence of K G = 0 is a direct dynamical consequence of the Poincaré-Hopf theorem [41] from which it follows that topological defects or vortices are required to accommodate vector order on a curved closed surface like the sphere. A covariant hydrodynamic treatment of active fluids on a curved surface has also been developed by Fily et al. [28]. These authors derived the continuum equations by coarse-graining a microscopic model of selfpropelled particles, which allowed an explicit computation of the transport coefficients in terms of microscopic parameters. The form of the continuum equations obtained in Ref. [28] is identical to those used here, the only distinction being that O(∇ 2 ) terms are neglected in that work, including the explicit K G term. In the following, we shall similarly neglect ∇ 2 terms.

III. POLAR FLOCK ON A SPHERE
As an example of a curved surface with constant positive curvature, we consider an active polar flock on the surface of a sphere of radius R. In local spherical polar coordinates {θ, ϕ}, the canonical metric and curvature on S 2 are The only non-vanishing Christoffel symbols are Γ θ ϕϕ = − sin θ cos θ and Γ ϕ θϕ = cot θ .
A. The steady state of a polar flock on the sphere At low mean density (ρ 0 < ρ c ), the isotropic phase with constant density and p = 0 is stable. For ρ 0 > ρ c , where the mean-field solution in flat space is a state of constant density and finite, but uniform polarization, on the sphere one obtains polar, spatially varying states. The density and polarization profiles for η = 2, now shown on the sphere. The color describes the density from the maximum (red) at the center of the polar band to ρc (blue) at the poles of the sphere. The polarization also vanishes at the poles.
Since the particle number is conserved, there can be no sinks or sources of flow. The simplest configuration allowed by the required conservation of topological defect charge that must sum up to the Euler characteristic χ = 2 of the sphere is then a circulating band wrapping around an equator, with two vortices of charge +1 at the two opposing poles. This yields a density band with polarization in the azimuthal direction and both density and polarization vanishing at the poles, consistent with the band state reported recently in simulations of polar particles on the sphere [29].
An explicit solution can be found analytically by assuming azimuthal symmetry, with ρ = ρ ss (θ), p θ = 0 and p ϕ = p ϕ ss (θ). The continuity equation is then satisfied identically. To simplify the polarization equation, we neglect the viscous terms as they are higher order in gradients (suppressed by 1/R 2 ) compared to the other terms arising from self-propulsion. In the microscopic realization of self-propelled polar particles with repulsive short range forces and aligning interaction studied in Ref. [29] this approximation corresponds to the regime where inter-particle repulsion (contributing to pressure) and active self-propulsion dominate over viscous and elastic stresses. We then neglect the laplacian terms entirely by setting ν = 0 (the bulk viscosity ν ′ drops out with our assumption of azimuthal symmetry). This leaves us with Writing X(θ) = ρ ss (θ) − ρ c , and seeking a solution with p ϕ ss = 0 we can eliminate p ϕ ss from the two equations to obtain is a dimensionless parameter that controls the shape of the solution, with η > 0 for the density profile to be a physical solution (note that sin θ > 0 over the entire range θ ∈ [0, π]). By symmetry the density will be maximum at θ = π/2 (the equator). Letting ρ ss (π/2) = ρ max we find We stress that the dependence on R has dropped out from Eq. (10), which therefore represents a universal density profile for an ordered flock on any size sphere. Finally, we express ρ max in terms of the average density ρ 0 ≡ ρ ss by requiring to obtain the final expression for the density profile as In order for this density profile to exist, we additionally require that |p ss | 2 > 0 and obtain As expected, an ordered flock only exists for ρ 0 > ρ c , and the magnitude of the steady state polarization and the density have the same inhomogeneous profile as shown in Fig. 1a, with the direction of polarization chosen spontaneously. This band solution is unrelated to the traveling bands found in flat space [51][52][53], which occur close to the meanfield transition and are absent deep in the ordered phase. The inhomogeneous solution obtained is simply the ordered flocking state on a sphere. The spatially inhomogeneous profile arises from the interplay of mass fluxes (∼ v 1 ∇ρ) and convective fluxes (∼ λp · ∇p) that cannot be driven to zero on a curved surface. Hence the spatial inhomogeneity is made inevitable by curvature.

FIG. 2:
The peak (maximum) density at the equator on a sphere, as we vary the mean density ρ0. For ρ0 < ρc, we are in the disordered phase with ρss(π/2) = ρ0. For ρ0 > ρc, we have a polar band with the density profile given in Fig. 1. The density reaches its maximum at the center of the band and grows with ρ0, with a slope Aη > 1. When the convective parameter λ → 0, η → 0, resulting in Aη → 1 and we go back to the homogeneous profile as in the flat plane.
The present solution is expected to break down within a region of angular width θ m ∼ exp(−aρ c R 2 /ν) (hence exponentially small on a large sphere) around the poles of the sphere, at the core of the vortices, as the elastic stresses will become important at short scales. Therefore the profile obtained is a robust and universal prediction of the continuum theory, similar to the rotating band seen in particle simulations of Ref. [29].
For an equilibrium polar or ferroelectric liquid crystal (say, a compressible lyotropic smectic-C film [55]), where the polarization is strictly an order parameter field and does not play the role of a velocity, the important convective nonlinearity in Eq. (3) is absent (λ = 0, though λ 2 and λ 3 can be present [56]) and the band solution we have here is absent (η = 0). In this case, even in the ordered phase, the density remains homogeneous and on a large enough sphere, we have nearly uniform polar order everywhere (|p ss | ≃ const.), but for two isolated defects at the poles, whose core size ξ ∼ ν/a(ρ 0 − ρ c ) is a microscopic length scale deep into the ordered state.

B. Linearizing about the steady state
Well below the mean-field transition (ρ 0 < ρ c ), the isotropic disordered state (ρ ss = ρ 0 and p ss = 0) is linearly stable at long wavelengths with fluctuations in the polarization relaxing quickly and density perturbations relaxing diffusively at long time, just as in the plane [57]. The curvature does not affect the disordered phase. It is only in the ordered state that we find novel excitations with non-trivial topological properties.
Here we consider the linear dynamics of small amplitude perturbations about the steady ordered state flock, letting ρ = ρ ss (θ) + δρ and p = p ss (θ) + δp. We focus on the long-wavelength propagating sound modes that are present even in the plane for an ordered flock [42], and continue to neglect all the viscous and elastic coupling. These are higher order in gradients and only give rise to damping of the sound modes. As the base state we are linearizing about is inhomogeneous, we additionally confine ourselves to a tangent plane linearization about a fixed latitude away from the poles (a preferred local coordinate system is picked out spontaneously by the polar order allowing for an unambiguous notion of latitude). Setting θ = θ 0 + y for a given latitude θ 0 , with θ 0 < π/2 corresponding to the northern hemisphere and θ 0 > π/2 to the southern hemisphere, relabelling ϕ as x, and letting δp θ → v and δp ϕ → u, we obtain with p 0 = p ϕ ss (θ 0 ) the azimuthal polarization at latitude θ 0 , which is finite as long as we are away from the poles, θ 0 = 0, π. Note that stability of the steady state requires η > 0, hence λ > 0. The only terms that can be negative are those proportional to cot θ 0 , arising from the Christoffel symbols, which changes sign as one crosses the equator at θ 0 = π/2.
Next we perform a Galilean boost to a comoving frame by letting x → x − λp 0 t (comoving with the logitudinal sound and not the flock itself), and relabelλ = λp 0 , The redefined parameters are summarized in Table I. Since the flock breaks Galilean invariance, this is not a symmetry operation, and yields Here m is a constant of fixed sign at any given nonequatorial latitude and changes sign across the equator, with m < 0 in the northern hemisphere and m > 0 in the southern half, vanishing only at the equator where θ 0 = π/2. We show below that a non-vanishing value of m leads to a band gap ( Fig. 3b and 3c) in the sound mode spectrum that acquires the necessary structure for non-trivial band topology. This, along with the vanishing of m at the equator, naturally suggests that the equator behaves as a "boundary" between two different "bulk" media (the northern and southern hemispheres), thereby allowing for localized topologically protected excitations on it. To simplify the notation we let |δΨ ≡ (δρ, u, v) and recast the linearized equations that control the linear stability of the steady state in the form of a Schrödinger like equation (in Fourier space, with Ψ(q) = d 2 r e −iq·r Ψ(r)), as The eigenvalues of H(q) directly give the sound mode frequencies (|δΨ ∝ e −iωt ). An important distinction compared to the Schrödinger equation is that the matrix H is not hermitian and therefore the linearized mode spectrum is not purely real, due to dissipative terms describing the overdamped dynamics and the absence of Galilean invariance. For m = 0 (θ 0 = π/2) the equations reduce to those of the planar case. Fluctuations in the polarization magnitude (u) are controlled by a fast mode that decays on microscopic time scales ∼ β −1 , The density (δρ) and the transverse goldstone mode (v) are the only slow modes that remain propagating at long wavelengths (as q → 0), where we have only kept terms to leading order in q.
For non-zero but small m = 0 (θ 0 = π/2), corresponding to the regions close to the equator in either hemisphere, the dispersion relations can be written as In the next section we analyze this mode structure.

IV. TOPOLOGICAL SOUND
We immediately see that the two branches of the propagating modes given by Eq. (26) have a gap at q = 0 of width ∆ = |ω + (0) − ω − (0)| proportional to |m|, with The terms explicitly involving m in the dispersion relations, responsible for the opening of the gap, are obtained only in the presence of both curvature and spontaneous active flow. In the plane static long-wavelength deformation of both the density and the broken symmetry mode leave the system unchanged. On a curved surface, in contrast, spatially uniform deformations of either "slow" field (δρ and v) cannot be static and invariably lead to dynamics in the system. As a result of curvature-induced forces, long-wavelength deformations of would be slow modes are required to have a finite frequency, resulting generically in the q = 0 gap of the sound spectrum, in sharp contrast to the conventional behavior of hydrodynamics in flat geometry [58]. It is useful to compare the effect at hand with one that occurs in geophysical flows. In a frame comoving with the flock, the finite curvature of the sphere plays a role similar to the Coriolis force that would be present for a passive fluid on a rotating sphere. In the case of the earth's atmosphere, this has recently been shown to give a gapped sound spectrum and equatorially confined Kelvin and Yanai waves [59] that are topological in origin [39]. In our active system no external flow or rotation needs to be imposed and the absence of Galilean invariance allows for independent tuning of the material parameters (such as λ) in order to probe regimes that are not accessible to passive fluids.
On times scales t ≫ β −1 we can slave the fast mode u to the slow fields, u ≃ αδρ/β + O(∇δρ) [85]. Upon eliminating u we get a reduced set of dynamical equations involving only δρ and v. After rescaling the wave-vector |αβ −λ|q x /β → q x , the linear matrix controlling the dynamics of δρ, v fluctuations is given by where µ = 2λα sin 2 θ 0 /β > 0 (θ 0 = 0, π) and s = sgn(α − λβ) (s = 0, ifλ = α/β). One can easily check that the eigen-frequencies of D(q) are exactly given by ω ± (q) (Eq. (26)) modulo the appropriate rescaling of q x . As D(q) is still non-hermitian, we need to evaluate right and left (adjoint) eigenvectors with the biorthogonality relation χ i |ψ j = δ ij . It is important to keep in mind the regime in which D(q) provides a valid approximation to the complete dynamics.
we can neglect the fast u mode and not worry about higher order terms in both q and m. This can be achieved deep in the ordered phase on a large enough sphere, in which case β is large, allowing for a large window of time in which the dynamics is dominated by D(q). With this set of simplifications, the linear dynamical matrix is always diagonalizable and ω + = ω − as long as m = 0 or q = 0 allowing one to adiabatically deform our model to have purely real eigenvalues by smoothly taking µ →v 1 (for µ = 0). In the process, the spectral gap remains open as long as m = 0. In order to establish the topological nature of the band structure, we compute the associated U(1) Berry gauge connection and curvature [60] along with the Chern numbers [61] for each band where s = sgn(α −λβ). The Chern number here is only quantized to a half integer as we work directly in the continuum long-wavelength approximation and the closing of the gap (for m = 0) only gives rise to one Dirac cone like structure (this is the contribution to the "parity anomaly" or Hall conductance associated to a single Dirac cone in a Chern insulator [62]). An appropriate regularization for large q, guarantees the Chern number to be an integer [38,63]. This calculation though still does have worth in predicting the correct number of topologically protected edge modes present when we stitch two of these regions with different Chern numbers together, via the bulk-edge correspondence. So each gap closing (change in sign of m), leads to a single localized edge mode, which as we shall see is unidirectional. As anticipated earlier we find that the Chern number of the acoustic band is different in the northern (m < 0) and southern (m > 0) hemispheres, vanishing at the equator (m = 0). Hence going across the equator, we have one gap closing (at q = 0) with a band inversion, leading to a single topological sound excitation localized at the equator. Note that the Chern number also vanishes for s = 0 which is obtained either when p 0 = 0 orλ = α/β, leading to a topologically trivial band structure. The first case corresponds to the absence of spontaneous active flow and the second to a partial restoration of Galilean invariance (in this limit, both density and goldstone mode excitations propagate with the same longitudinal speed). So the vanishing of the Chern number and associated band triviality for s = 0 is not due to the closure of a gap, but instead due to a restoration of symmetry. As β = 2bR 2 p 2 0 sin 2 θ 0 depends on the latitude at which we are, the vanishing of s (forλ = α/β) is a condition on the polar angle θ 0 . There is a critical density such that for ρ c < ρ 0 < ρ * , s = 0 on the entire sphere. Deeper into the ordered state (ρ 0 > ρ * ), there are two latitudes at angles θ ± such that sin η θ ± = (ρ * −ρ c )/(ρ 0 −ρ c ), at which s = 0. Even though the band topology changes as we cross the latitudes at θ ± (Eq. (33)), the spectrum remains gapped throughout and hence we do not have any gapless excitations localized at θ ± . The change in the Chern number across these special latitudes is due to an accidental additional symmetry (Galilean invariance) instead of the gap closing, thereby circumventing the bulkboundary correspondence. This is a well known point in quantum topological insulators, only realized here in a peculiar fashion as the protecting "symmetry" varies spatially in a single sample.
To summarize there are three crucial ingredients in this system that lead to and protect the topologically nontrivial band structure -• Breaking of time-reversal symmetry by the active polar flow. Changing the direction of spontaneous polarization (flow) changes the sign of s (as p 0 → −p 0 leads to α → −α,λ → −λ, β → β and µ → µ).
• The presence of the convective nonlinearity λ = 0; an equilibrium passive polar liquid crystal will not therefore exhibit these modes.
• The curvature of the base surface which opens up a gap in the sound spectrum. Changing the gaussian curvature exchanges the regions with positive and negative "m". This is entirely analogous to the Haldane model [62] where the closing of the gap at the Dirac point is protected by time-reversal symmetry, which when broken by the local magnetic field leads to a band structure with a non-trivial topology. Though the active system is not hermitian with purely real frequencies, the structure of the localized equatorial mode for varying m(y) is adiabatically connected to its hermitian analogue [64], the Jackiw-Rebbi soliton [65].
where ψ 0 is a normalization constant. The edge mode spectrum ω edge (q x ) = sq x corresponds to a pure oneway density wave that connects the two bulk bands (see Fig. (4)). This edge mode is valid when m(y) → ±m 0 (m 0 > 0) for y → ±∞.
On the sphere, reverting back to angular coordinates {θ, ϕ}, m(θ) = − cot θ, which is positive in the southern hemisphere for θ > π/2 (y > 0). This gives a chiral [a n e in(ϕ−αt/β) + c.c], (36) where a n are complex constants depending on the initial perturbation applied. A snapshot of this density mode is shown in Fig. 6 for n = 6 (and all other a n vanishing). Equation (36) defines a localization length ℓ loc = R/η set by the curvature and material parameters of the active fluid and essentially given by the ratio of longitudinal to transverse sound speeds. Note that the result given in Eqs. (35) and (36) applies for λ 3 = 0. The general case of λ 2 , λ 3 = 0 is given in Appendix A and yields a different localization width for the equatorial mode. This topological edge mode propagates unidirectionally in the direction of the flock and is robust to disorder because there are no reverse channels into which it can scatter (though it will eventually dissipate due to viscous and elastic damping). Unlike a Galilean invariant fluid for which η = 1, here η and therefore the localization length can be tuned by varying the system's parameters, although the shape of the steady state profile remains unchanged.

V. POLAR FLOCK ON A NEGATIVE CURVATURE SURFACE
The presence of such topological excitations is generic in that they will always be present when one has a polar flock on a curved surface that looks locally like a surface of revolution [86]. We illustrate this point on a catenoid, a surface with non-constant negative gaussian curvature.
In local coordinates {y, ϕ} (ϕ once again being the periodic azimuthal direction), the metric and gaussian curvature on a catenoid are g = R 2 cosh 2 y(dy ⊗ dy + dϕ ⊗ dϕ) , where R is the radius of curvature at the neck of the catenoid. In contrast to the sphere, the gaussian curvature here is both negative and spatially varying. The only non-vanishing Christoffel symbols are Γ ϕ ϕy = Γ y yy = −Γ y ϕϕ = tanh y .
Taking the same approach as for the sphere, neglecting viscous and elastic stresses, we consider an azimuthally symmetric ansatz for the steady state polar flock : ρ = ρ ss (y), p y = 0 and p ϕ = p ϕ ss (y). One can easily verify that for ρ 0 > ρ c the steady state density profile is then where B η < 1 is a constant that depends on η and the height of the catenoid (which unlike the sphere is not compact and has to be taken finite). The details of the computation are given in Appendix. B. In contrast to the sphere, which had a polar band with maximum density at the equator, the polar flock density is lowest at the neck of the catenoid (y = 0), increasing on either side as one moves away from it. The corresponding polarization profile is given by The density and polarization profiles are plotted in Fig. 5a. Below the mean field transition (ρ 0 < ρ c ), we recover the isotropic disordered phase (ρ ss = ρ 0 , p ss = 0). Linearizing about this steady state, one finds that the equations governing the propagation of sound modes on the catenoid are essentially identical to that on the sphere (Eqs. (14), (15), and (16)), but with modified parameters. As a consequence of the negative curvature, the most important change is that m = −2 tanh y is positive below the neck of the catenoid (y < 0) and negative above (y > 0), vanishing right at the neck (y = 0). This leads to a chiral mode of goldstone fluctuations localized at the neck of the catenoid (δρ = 0), which written in the lab frame is given by δv edge (ϕ, y; t) = sech 2 y n≥0 b n e in(ϕ−λt) + c.c , (42) where as beforeλ = λp 0 and b n is a complex coefficient determined by the initial perturbation. A snapshot of this mode is shown as well in Fig. 6. This mode propagates in the same direction as the flock, but with a different speed and is topologically protected. The localization length ℓ loc = 2R is controlled by the scale of the curvature in the system and seemingly independent of the material parameters of the flock.

VI. CONCLUSION
The frustration associated with the interplay of curvature and order [30,66,67] has many consequences for crystals [68][69][70][71][72], tethered membranes [73,74], liquid crystalline membranes [75][76][77], and jammed and glassy systems [78,79]. The presence of activity adds an entirely new non-equilibrium dimension to the whole story, allowing for completely new physics arising from competing order, curvature and the active drive. The active polar fluid is peculiar as the polarization order parameter is also a velocity that advects the fluid around [4], with this dual role being at the heart of many of the phenomena we have explored in this article. In particular, unlike a superfluid film on a curved substrate (where the order parameter can be parallel transported trivially) [80], the fact that the order parameter of the flock is a physical velocity implies that it advects itself non-trivially in the presence of curvature. This is nothing but a restatement of the physical fact that self-propelled particles move persistently along geodesics (in the absence of interactions), which get "lensed" by curvature, whereas passive polar particles don't do any such thing.
In order to handle these points we generalized the continuum Toner-Tu model for an active polar fluid, to an arbitrary curved surface and found new terms that are absent in flat space. In general, studying ordering phenomena on curved surfaces is rather complicated, even in equilibrium. This is when symmetry is often a very useful guiding tool, using which we explicitly computed the ordered phase of an active polar flock on two surfaces, the sphere and the catenoid. The continuum model affords us the privilege of having model independent predictions and we find many of the features of the steady ordered state of a polar flock are quite generic, with positive curvature surfaces having density profiles with a maximum, decreasing on either side of the peak, while negative curvature surfaces have the density profiles of the opposite kind, being minimal in the interior and increasing towards the boundary. Finding such spatially inhomogeneous exact solutions to the covariant Toner-Tu model is definitely a crucial starting point to being able to understand how the phenomenology of active matter in flat space translates to its curved variants.
In addition to the steady state with spontaneous flow, flocks in flat space also have dissipative sound modes with a linear (but angle dependent) dispersion [42]. It is here that all three players : curvature, order and activity come together with dramatic consequences. The presence of curvature gaps the sound mode spectrum at long wavelengths leading to a band structure with nontrivial topology protected by the broken time-reversal and Galilean symmetry in the system. We demonstrate this by computing the Chern number for the bands and show that this is a generic feature of flocks on curved surfaces. The most interesting result the non-trivial band topology has is to localize "edge" modes of density or goldstone mode fluctuations along special geodesics on the surface, at which the gap in the sound spectrum vanishes. The rather novel feature here is that the system is not artificially engineered as a metamaterial with some underlying lattice structure [33,81], nor does it require any external forcing or fields of any kind [34]. The spontaneous flow is generated by activity breaking both time reversal and Galilean invariance simultaneously, while the curvature is responsible for the spectral gap in the ordered phase.
Topological excitations of the type described here are "protected" in the sense that they are robust against static perturbations and heterogeneities in the medium through which they propagate. While quantifying the limits of such topological protection in active systems will require numerical work and remains to be explored, our work demonstrating that such topologically protected propagating modes are a generic consequence of active flows on curved surfaces raises the question of whether nature may use this mechanism to guide and direct the robust transmission of intercellular physical forces in curved environments. It is therefore tempting to offer some speculation to the possible relevance of our findings to biology. In a number of developmental phenomena, from wound healing to morphogenesis and organ development, living cells migrate collectively, offering an intriguing realization of a polar active fluid. While a full understanding of the mechanisms that regulate collective cell migration is still out of reach, it is now widely recognized that the transmission of physical forces plays an important role, alongside biochemical signalling. For instance, propagating mechanical waves have been shown to mediate cooperative force transmission among epithelial cells in wound healing assays [82]. In many biological processes cell motion takes place on curved surfaces, as in cell renewal and repair in the highly folded intestine [83] and the shaping of the early limb bud in developing embryos [84]. The effect of curvature on the dynamics of epithelial cells is beginning to be explored in vitro by examining collective cell migration on cylindrical capillaries of varying radii [21]. While cylinders have zero gaussian curvature which would not yield topologically protected states, these experiments clearly demonstrate that curvature affects cell morphology and dynamics by enhancing cell speed and cell extrusion. A more direct application of the work described here would be to a polar version of the active nematic vesicles described in Ref. [10]. Here active vesicles were engineered by confining an active suspension of microtubule-kinesin bundles to the surface of a lipid vesicle. The interplay of activity and curvature yields a number of dynamical structures, including spontaneously oscillating defect textures and folding nematic bands, and ultimately activity-driven shape deformations of the vesicle. Our work may also be relevant to the physics of cell membranes that are activated through coupling to the polymerizing acto-myosin cortex, as modeled in recent work by Maitra et al. [31]. Finally, in the spirit of colloidal crystals on curved interfaces and reconstituted active systems, one might also envision synthetic experimental realizations of the topological sound modes investigated here, by depositing active Janus colloidal rods or active bio-filament motor complexes on the surface of a droplet or a vesicle shell.

VII. ACKNOWLEDGMENTS
We thank Brad Marston for inspiring this work through a beautiful talk delivered at KITP based on Ref. [39]. We also thank Rastko Sknepnek and Vladimir Juričić for useful discussions. This work was supported by the National Science Foundation at Syracuse University through awards DMR-1609208 (MCM & SS) and DGE-1068780 (MCM) and at KITP under Grant No. NSF PHY-1125915. All authors thank the Syracuse Soft Matter Program for support and the KITP for its hospitality during the completion of some of the work.
Appendix A: Steady state and linearization for λ2, λ3 = 0 Including the two additional λ 2 p∇ · p and λ 3 ∇|p| 2 nonlinearities, the equation for the polarization order parameter is modified to We don't include the viscous terms ν, ν ′ , consistent with our approximations in the main text. For the azimuthally symmetric ansatz (on the sphere and the catenoid), ∇ µ p µ ss = 0 identically and hence the λ 2 terms do not affect the steady state profile. On the contrary, the λ 3 ∇|p| 2 term acts as an additional polarization dependent contribution to the scalar pressure P ∼ v 1 ρ − λ 3 |p| 2 , which if large can lead to density and splay instabilities [49]. We shall disregard this and only work in the regime where λ 3 is not large enough to destabilize the entire system. In-cluding it, the steady state equations on the sphere (Eqs. (6), (7)) get modified to λ sin θ cos θ(p ϕ ss Once again, setting X(θ) = ρ ss (θ) − ρ c , we get the same equation as before, only now with a modified coefficient depending on λ 3 .
with the exponent now changed to η ′ = λa/(v 1 b − aλ 3 ). For η ′ > 0 (to have a physical density profile), we require v 1 b > aλ 3 , which is nothing but a condition to have stable pressure and a positive compressibility. Hence the effect of λ 3 = 0 is to only change the density profile through the η exponent, the functional form remaining the same. This is true even for the catenoid where the exponent is the same as on the sphere and given by η ′ = λa/(v 1 b − aλ 3 ). Note that this inhomogeneous profile does not exist for an equilbrium polar liquid crystal for which λ = 0 (but possibly λ 2 , λ 3 = 0 allowing for spontaneous splay [56]), leading to η ′ = 0. Hence the inhomogeneous steady state we obtain is only possible in an active system.
We use the same notation as we used in the main text (along withλ 2,3 = λ 2,3 p 0 ). At long times (t ≫ β −1 ), the polarization magnitude u is still a fast mode and to leading order, it gets slaved to the density fluctuations in the same fashion as before (u ≃ αδρ/β). Consistent with the level of approximation used in the main text (Eq. (31)), the long time and long wavelength dynamics is governed only by the two slow modes δρ and v, with a dynamical matrix D(q) of the same form as found in the absence of λ 2 and λ 3 .
It is easy to see now that the profile of the localized equatorial density mode on the sphere, which was ∝ sin µ ′ /v ′ 1 θ (for λ 2,3 = 0, µ/v 1 = η as given in Eq. (36)) is no longer the same as the steady state density profile of the flock (∼ sin η ′ θ), when λ 3 = 0.
A similar result also holds true for the catenoid. Hence we find that, even upon including additional convective non-linearities (which are lower order in gradients compared to the viscous terms), all of the qualitative properties of the steady state and the topologically protected modes remains the same, with the only modification being a more detailed dependence of the localization length on some of the material parameters in the system.

Appendix B: Polar flock on the catenoid
For an azimuthally symmetric ordered steady state on the catenoid, just as on the sphere, we neglect the viscous and elastic stresses, and use the ansatz : ρ = ρ ss (y), p y = 0 and p ϕ = p ϕ ss (y). Plugging this into the Eqs. (1) and (3), we find that the continuity equation is satisfied identically and Eq. (3) reduces to (for ν = 0) p ϕ ss a(ρ ss − ρ c ) − bR 2 cosh 2 y(p ϕ ss ) 2 = 0.
Setting X(y) = ρ ss (y) − ρ c , we solve the equations in the same fashion as before to get ∂ y X = η tanh yX, where η = λa/bv 1 (the same exponent as on the sphere). The steady state density profile is then ρ ss (y) = ρ c + (ρ min − ρ c ) cosh η y.
where ρ min is the minimum density of the flock attained on the neck of the catenoid (y = 0). Unlike the sphere, the catenoid is not a compact surface, so in reality one would have a finite sample with boundaries. The mean density ρ 0 is given by the spatial average of the steady state profile, where · denotes a spatial average over the entire surface. For a catenoid of height L (euclidean height in the zdirection when embedded in R 3 ) and radius of curvature Writing B η = 1/ cosh η y < 1, we obtain the density profile quoted in the main text (Eq. (40)) We expect the viscous and elastic stresses to be less important on a weakly curved surface close to the neck, in particular when the characteristic scale of curvature (∼ R) is much greater than the equilibrium correlation length (ξ ∼ ν/a(ρ 0 − ρ c )). Additionally the density and polarization (along with their gradients) grow larger as we go away from the neck. So close to the boundaries of a large sample, one would have to account for higher order nonlinearities along with the elastic stresses, which would then become important.
At long time (t ≫ β −1 ), the fast polarization magnitude u decays and is slaved to the density field u ≃ αδρ/β (to lowest order) and the slow dynamics at long wavelengths is dominated by only now with µ = αλ/β > 0. Hence at the same level of approximation used earlier for the sphere (neglecting viscous stresses and the parameter regime given in Eq. (31)), the long time dynamics of sound excitations in a polar flock on a curved surface is generically described by equations of the form given above, or consequently by the linear dynamical matrix D(q) (Eq. (28)), possibly upto some coordinate rescaling.
As the only modifications are in the definitions of the parameters, many of the predictions made in the case of the sphere apply here too. In particular the sound mode spectrum is still gapped at q = 0 for non-zero m and the bands have a non-trivial topology given by the Chern numbers C ± (see Eq. (33)). As m = 0 at the neck of the catenoid (y = 0), changing sign on either side, we have one topologically protected mode localized at the neck. A consequence of the negative curvature of the surface is that, in contrast to the sphere, m < 0 for y > 0. Due to this, the edge mode takes on a different form (in the comoving frame) |δΨ edge = ψ 0 e´y 0 m(y ′ )dy ′ +iqxx 0 1 .
Now the edge mode is a localized mode of transverse goldstone fluctuations with density fluctuations completely absent. Additionally, the edge mode spectrum is ω edge = 0 to lowest order in q x implying that the edge mode is stationary in the comoving (with speed ∼ λp 0 ) frame. This too connects the two bulk bands and is topologically protected. Using m(y) = −2 tanh y for the catenoid, this gives the profile of the localized mode, the lab frame version of which is quoted in the main text (Eq. (42)).