{'fraction_non_alphanumeric': 0.039859666920280666, 'fraction_numerical': 0.014265787471468425, 'mean_word_length': 4.363523010655181, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chirality-driven edge flow and non-Hermitian topology in active nematic cells

Many of the biological phenomena involve collective dynamics driven by self-propelled motion and nonequilibrium force (i.e., activity) that result in features unexpected from equilibrium physics. On the other hand, biological experiments utilizing molecular motors, bacteria, and mammalian cells have served as ideal setups to probe the effect of activity in materials and compare with theory. As has been established, however, biomolecules are chiral in nature, which can lead to the chiral patterning of cells and even to the left-right symmetry breaking in our body. The general mechanism of how the dynamics of bio-matters can couple with its own inherent chirality to produce macroscopic patterns is yet to be elucidated. Here we report that cultured neural progenitor cells (NPCs), which undergo self-propelled motion with nematic cell-to-cell interactions, exhibit large scale chiral patterns when flowing out from containers made by gel. Moreover, a robust chiral cell flow is produced along the boundary when the NPCs are cultured on substrates with edges. Perturbation by actomyosin inhibitors allowed control over the chirality, resulting in the switching of the direction of the chiral patterning and boundary flow. As predicted by a hydrodynamic theory analogous to the non-Hermitian Schrodinger equation, we find an edge-localized unidirectional mode in the Fourier spectrum of the cell density, which corresponds to the topological Kelvin wave. These results establish a novel mechanism of flow that emerges from a pool of bipolar cells, and demonstrate how topological concepts from condensed matter physics can naturally arise in chiral active systems and multi-cellular phenomena.

The mechanism of how chiral motions at the level of components can give rise to robust macroscopic patterning has recently gained much interest in nonequilibrium physics [21][22][23][24]. In biological systems, the chirality of cytoskeletons and molecular motors are typically responsible for the chiral patterns observed at the cellular [10-13] and multi-cellular levels [14], with some shown to be directly causing the left-right asymmetry in our body plans [15][16][17][18].
Recent works [20,[25][26][27][28][29][30] have elucidated how chirality in classical systems can lead to situations similar to topological insulators, a phenomenon heavily studied in quantum systems [31]. In these settings, waves are observed at the boundary or the interface of the system in a predictable manner from the celebrated bulk-edge correspondence [32] while the net flow is suppressed in the bulk. An intriguing example is the Kelvin wave in geophysics, which explains the robust equatorial wave driven by the chirality of the earth rotation and has recently been found to have a topological origin [19].
Whether such robust waves due to chirality and topology exists in biological systems is an interesting question. Indeed, spontaneous flow is found in a wide range of scales in biology, from the cortical flow [33], chiral flow in bacteria [34], up to the collective cell migration in vivo [15,35]. Yet, to our knowledge, there is still no example of topological waves playing a role in the biological context.
Here we use neural progenitor cells (NPCs), an active nematic system [6-8, [36][37][38][39], as a model to investigate how chirality of biomaterials can couple with activity to generate unidirectional flow. We first establish that the collective dynamics of NPCs are chiral by an assay involving cells flowing out from a container made by gel. Next, we observe how the NPCs under the confined geometry with circular and linear boundaries exhibit chiral edge currents. By numerical simulation and the analysis of the hydrodynamic theory, we find that the edge flow is a natural consequence of a chiral active nematic confined within a region with boundary. We further find that the fluctuation in the cell density at the edge shows a pattern of the topological Kelvin mode, resembling the situation in geophysics [19]. The results clarify the scenario of chiral cells producing topological edge flow, thus expanding the possibility of applying condensed matter concepts in understanding multi-cellular phenomena.

Chiral patterning in neural progenitors
We first prepared a dense cell suspension in Matrigel, a commonly used extracellular matrix which is mainly composed of the laminin protein. When the gel droplets were deposited on glass and incubated in culture medium (Fig. 1a), cells spontaneously migrated out from the gel within a day. Because of the high cell density, the cells were nematically aligned soon after exiting from the gel (Fig. 1b).
After two days of culturing, when the cells had spread out from the gel due to migration and cell divisions, the cell population showed a clear chiral pattern (Fig. 1c). The angle of alignment, quantified by applying the tensor method [40,41] to the phase contrast images, became tilted at the positions further away from the center of the gel (Fig. 1d, left). The angle of tilting depended on the clone of the NPCs used in the experiment (Extended Data Fig. 1a), but the direction of chi-rality was always the same; the cell alignment pattern tilted toward the right respective to the radial axis. We asked if perturbations on cytoskeleton dynamics can affect the chiral patterning, as has been tested in experiments using cell lines [14] and tissue dynamics in chicken development [15]. The high concentration of cytoskeleton inhibitors prevented the spreading of the cells by blocking cell migration and division, and also disrupted the coherence of the nematic patterning (see Methods and Extended Data Fig. 1b). Nevertheless, several inhibitors, when used at low concentration, modulated the chirality of the pattern while maintaining the nematic ordering (Fig. 1d,e,f, Extended Data Fig. 1c). Of the tested inhibitors, Jasplakinolide, an actin stabilizer, and Latrunculin A, an actin polymerization inhibitor, weakened the right-handedness of the cell migration pattern [14]. On the other hand, Blebbistatin, a myosin II specific inhibitor, enhanced the right-handedness. Microtubule related inhibitors had less effect on chirality, and merely decreased the coherence of nematic order at higher concentrations.
To see if the chirality exists at the single-cell level, we tracked the motion of NPCs at low cell density (Extended Data Fig. 1b). Although the motion of the cells is markedly different in the sparse condition, the cells showed a clear bias toward turning right (∼ 0.1 rad/hour, Extended Data Fig. 1d,e,f). However, we found no evidence of this single-cell-level chirality being affected by cytoskeleton inhibitors at least in the dosage range where cells were able to migrate. This indicates that the inhibitors are affecting the chirality of cell-tocell interactions, with possible mechanisms being the friction between rotating cell bodies as in the case of growth cone filopodia [13], or the asymmetric sliding between neighboring cells [42].

Chiral edge flow
We next asked if the chirality can be observed under other settings by employing a microcontact printing strategy to confine the cells [43] (see Methods). When we confined the cell migration area with a circular boundary by laminin coating and passivation (Fig. 2a), the cells showed unidirectional flow at the edge (Fig. 2b). By quantifying the average velocity by tracking the cells using the cell nucleus signal (H2B-mCherry), we found that the flow is strongly localized at the edge. This is remarkable since the individual cells do not have a net directionality; cells can migrate in either direction that respects the bipolar cell shape. Indeed, the cells flip their direction of motion stochastically once every 2-3 hours on average, which is why net flow is weak in the bulk except for the flow around topological defects [8]. To confirm that this edge flow is arising from the same chiral effect in the macroscopic patterning observed in the gel experiment we perturbed the cells by Jasplakinolide. As expected, the edge flow switched its direction (Fig. 2c,d).
To further quantify the chiral edge flow in an even simpler setting, we again used the laminin stamps to contain the cells in stripes with various widths, L, which allowed the cells to show near-perfect nematic order (Fig. 3a)  cells showed similar edge flow as in the case of circular shaped boundaries; on the right edge, the cells were moving toward the positive y-direction, whereas on the left side the cells were migrating toward the negative (Fig. 3a). A similar pattern of cell flow has been observed in retinal progenitor cells [38], although the chirality in NPCs was consistent and the edge flow was an order of magnitude faster, likely due to a different mechanism of chiral flow generation. The cell density was higher near the boundary compared to the center region of the stripe (Fig. 3b, Extended Data Fig. 2a), and the cell flow in the y-direction was also concentrated at the edge (Fig. 3c). Both of these features were observed even when cell division was blocked by Mitomycin C (Extended Data Fig. 2b), indicating that the flow and the inhomogeneous cell density is not due to the pattern of cell growth (Extended Data Fig. 2e). Perturbation by Blebbistatin and Jasplakinolide resulted in slightly enhanced and reversed edge flow, respectively (Fig. 3d,e, Extended Data Fig. 2c,d), as expected from the gel-based assay (Fig. 1d).
Numerical simulation and theory of chiral active nematics To investigate the mechanism of the chiral edge flow, we conducted numerical simulation of an agent-based model that describes a dense active nematic system in two-dimensions. In this model, the cells are represented by agents undergoing unidirectional motion with velocity v 0 , stochastic flipping with average flipping time τ, and rotational diffusion set by the constant D φ (Fig. 3f). The agents interact with each other through repulsion and nematic alignment (Fig. 3g), with additional chirality controlled by the parameter θ, which quantifies the rotation rate of the angle of the agents (Fig. 3f, see Methods for the detail of the numerical simulation). The pattern within the stripe confinement, where the boundary condition was set to make the agents align to the edge, resembled the features observed in experiment when θ and τ were both nonzero (Fig. 3g, Extended Data Fig. 3a,b). The density profile, the pattern of alignment, and finite velocity near the edge were also reproduced when we set θ = 0.2 rad/hour. As expected, this value of chirality is the same order but larger than the chirality estimated from the cell tracks under low density (Extended Data Fig. 1f).
From the experimental observations and numerical simulations, it is clear that the chiral symmetry breaking, θ 0, is key in the emergence of the edge flow. To seek a theoretical explanation of this phenomenon, we turned to the hydrodynamic theory of dense active nematics, where the cell density, cell flow, and nematic order undergo time evolution simultaneously (see Supplementary Information). Without chirality, the hydrodynamic equation has a steady-state solution describing a perfect nematic order, where the cell density is uniform (ρ = ρ 0 ), velocity field is zero (P x,y = 0), and the nematic order is parallel to the edge (Q yy = −Q xx = r > 0 and Q xy = 0, where Q is the nematic tensor and r is positive parameter that quantifies the extent of alignment) [44][45][46]. By adding small cell chirality, we obtain a linearized equation of time evolution, (1) Here, Ψ = (δρ/ρ 0 , δP x , δP y , δQ xx , δQ xy ) T describes the fluctuations of density (δρ), velocity field (δP x,y ) , and nematic tensor (δQ xx , δQ xy ) around the perfect nematic ordered state. M is a 5 × 5 matrix including spatial derivatives, and S = (0, 0, 0, 0, 2θr) is a constant global driving that arises due to θ (see Supplementary Information). Importantly, M also includes θ-dependent terms, which are analogous to the Coriolis force in geophysics [19] but with the origin being the chiral nature of the cell rather than the rotation of the earth. Within the linear regime, the edge flow can be explained as the combination of chirality-induced deformation of nematic pattern and the pattern-induced active flow [47]. First, the chirality of the cells induces a tilt in the steady-state nematic ordering (i.e., ∂ x δQ xx,xy 0). Then, the spatial non-uniformity of the nematic order leads to the active force field, which is captured in the steady-state equation: where η = v 0 /(4/τ + 2D φ ) (see Supplementary Information). sults of experiment and numerical simulation (Extended Data Fig. 4). We note, however, that the situation of the experiment is expected to be out of the linear regime for large L due to the relatively large chirality, θ. Nevertheless, we here study the linear equation expecting that the basic physics is preserved even when taking in the higher-order terms in the hydrodynamic equation (see Supplementary Information).

Spectrum of the linearized dynamics
Equation (1) can be considered as a non-Hermitian Schrödinger equation with an external field. We therefore calculated the complex spectrum of the non-Hermitian Hamiltonian, H = iM, under the stripe boundary condition similar to the experiment. We found that the edge modes that are most significant in the cell density fluctuation has a pattern similar to the topological Kelvin mode obtained in fluid mechanics [19,20] (Fig. 4a). By changing the extent of chirality, we further found that the amplitude of the edge localization becomes weaker as θ becomes close to zero, whereas the slope of the modes in the wavenumber-frequency plane does not change according to θ (Extended Data Fig. 5c,6c,7a). This is also consistent with the edge-localized topological Kelvin mode, which we can identify through a simplified linear equation that assumes decoupling between the velocity and the nematic fields (see Supplementary Information). We obtain, as complex eigenvalues E of H, for k 2 y > k 2 0 and zero otherwise, where k 0 = [8(1 + r)] −1/2 η −1 . ReE + (k y ) is plotted in Fig. 4a (dashed lines), which is in agreement with the spectrum of the full model. Therefore, we expect that the mode of the density fluctuation at the edge should follow ReE ± (k y ) ± √ (1 + r)/2v 0 k y for large k y , which is a θ-independent dispersion relation.
To compare with this result, we calculated the power spectrum of the cell density by using the fluorescent signal of H2B-mCherry (Fig. 4b). First, we found that there is a pattern of unidirectional mode observed only at the edges and not in the bulk region (Fig. 4c). Second, the spectrum calculated from experiments with cytoskeleton perturbations showed that the group velocity of this unidirectional mode, obtained as the fitted slope to the spectrum for large k y , is independent of the steady-state edge velocity and is rather correlated with the single cell velocity (v 0 ) of the NPCs under each condition (Fig. 4d,e). Lastly, this unidirectional mode weakened (Fig. 4f) and switched direction (Fig. 4g) upon the application of Jasplakinolide. These results are consistent with the interpretation that the unidirectional edge mode observed in the fluctuation of the cell density is the topological Kelvin wave. We also found that the spectrum of the particle density fluctuation in the agent-based model follows a similar pattern (Extended Data Fig. 7b,c).
The main difference between our non-Hermitian setup and the previously studied Hermitian system [19] is that there is no band gap observed within the realistic range of parameters. Nevertheless, we found that the band gap can be induced by assuming larger chirality in the same equation, where the bands can be found to have non-trivial topological Chern numbers (Extended Data Fig. 5, see Supplementary Information). The mode described by Eq. (3) existed irrespective of the energy gap, which matched with one of the topological edge modes in the case where the gap was open. The gapped Hamiltonian can be continuously deformed to a Hermitian Hamiltonian without closing the gap, and the property of the edge localized mode did not change along this path. The extent of edge localization of the unidirectional modes depended on θ, and the regime of the experiment and simulation (θ = 0.2 rad/hours) is where there is significant edge localization of the modes without the energy gap. Therefore, we find that the edge modes generated by the NPCs, the agent-based model simulation, and the non-Hermitian theory are all correspondents of the topological edge modes observed in Hermitian and gapped systems [19].
In the band structure analysis of the non-Hermitian model, we further noticed that many of the modes identified in the bulk band for the periodic boundary condition are localized at the edge under the open boundary condition (Extended Data Fig. 5e). This is a phenomenon called the non-Hermitian skin effect, which has also recently been shown to have a topological origin [48,49]. The direct consequence of this effect in chiral active matter is yet to be elucidated.

Discussion and conclusion
In this work, we have shown that the chiral nature of the NPCs can produce not only a chiral spiral pattern in the collective cell dynamics but also a robust unidirectional steady flow and the Kelvin wave at the edge in a confined culture setting. Results of numerical simulation and theory show that the dynamics at the edge, i.e., the geometry-induced steadystate active flow and the topological edge mode, are universal features of chiral active nematics.
Recent works have discovered that specific molecules such as myosins Ic, Id, and V are responsible for the chirality of cell motion [13,16,18], which implies that the change in the balance between these molecules and less chiral myosins such as myosin II can lead to changes in the extent of chirality presented at the cell level. Although the direction of cell chirality can be different across cell types and species, we expect to find similar edge effects in other real biosystems since the mechanism of flow and wave generation does not depend on the detail of the cell biology. It will be particularly interesting to find in vivo examples of topological modes, such as surface flow inside cells and cell flow in tubes and organs. Phys. Condens. Matter 32, 193001 (2020  Gel drop assay Cells were detached from dishes and spun down to make a pellet of 2.5 × 10 6 cells. The pellet was suspended in 50 µl icecold Matrigel which was stored in aliquots at -20 deg. The cell suspension in gel was plated onto glass base dishes using pipette tips which were stored in the refrigerator until use. Each drop consisted of 2 µl of suspension. After 3 minutes of incubation in 37 deg, the plate was filled with warm culture medium with additional Matrigel suspension (0.15%) and left overnight in the incubator. The next day, the media in the glass base dishes were changed to fresh media with or without the inhibitors. After two to three days of culturing, the glass base dishes with the gel drops were placed under the fluorescent microscope (Zeiss AxioObserver7). Fluorescence and phase contrast images typically of size 7.3 mm x 7.3 mm were obtained by the 10x lens using the tiling acquisition setup.
Image analysis of the chiral patterning Colour maps representing the tilt angle of the cell alignment were obtained by the tensor method [8, 40,41]. For an acquired phase contrast image I(r), where r = (x, y) denotes the two-dimensional position (in pixels), we calculated the differ-ential tensor: with After applying a spatial Gaussian filter to J to obtain where the size of the filter was 65 µm, we calculated the local principal angle of alignment θ(r) := 1 2 arctan and the coherence The coherence, which takes values between 0 and 1, quantifies the extent of anisotropy in the image. For the phase contrast images of the NPCs, the perfectly aligned situation corresponded to C ∼ 0.3, as shown in Extended Data Fig. 1b. For Fig. 1d, we masked out the regions of the gel and the regions where cells are absent by setting upper and lower thresholds in the intensity of the H2B-mCherry channel. We calculated the mean coherence and tilt angle by averaging the gradient tensor over all the pixels in each image using the same mask (Fig. 1d,e). All analyses were conducted using a custom Python code.
Low cell density experiment NPCs tend to become unhealthy under low cell density conditions, which can be rescued by applying conditioned media. We prepared conditioned media by taking the supernatants from the plates with sub-confluent cell density. Typically, we plated NPCs at a quarter of confluent density, changed medium the next day, and collected the supernatant the day after. The supernatant was spun down 800 g 5 minutes, filtered 0.22 µm, and aliquoted and frozen down in -20 deg.
For low cell density experiments, we plated 2 × 10 4 cells per glass base dish in 1.5 ml thawed and warmed conditioned media with 0.15% Matrigel and waited three hours until the cells adhered to the dish surface. We added 0.5 ml of conditioned media with or without cytoskeleton inhibitors and left the dish overnight in incubation. The next day, the glass base dishes with low cell density were put on the confocal microscope (LSM800, Zeiss) to take tiled images of the fluorescent channel for H2B-mCherry every 5 minutes using a 10x lens. For the tracking of cells, we used the TrackMate plugin included in the Fiji package [50].
From the obtained single-cell tracks, we calculated the average angle change as a function of time (Extended Data  Fig. 1f) using a custom Python code. Using the angle of displacement θ i (t) of the i-th cell at between time frames t and t + 1, we calculated the average angle change by where the average ... is taken over the different cells (i) and time (s).

Microcontact printing experiment
For the micropatterning experiment [43], we first designed photomasks and ordered printing on positive film sheets (Tokyo Lithmatic Corporation, Tokyo). Next, we prepared a substrate of negative photoresist (SU-8 3025, MicroChem) on silicon wafer by spin coating, applied the photomask after soft baking, and exposed to UV light using an LED spotlight. After developing and baking the photoresist pattern by the standard procedure, we made stamps by pouring polydimethylsiloxane (PDMS, SYLGARD 184 Silicone Elastomer Kit, Dow Corning 98-0898) onto the pattern and incubating at 80 deg for more than two hours. To make dishes with micropatterned laminin, we first cleaned the PDMS stamp surface using an oxygen plasma cleaner for two minutes, deposited 50 µl of cold 20% Matrigel solution (in PBS) on its surface, and incubated 30 mins at room temperature. We then removed the matrigel solution by an aspirator, and gently placed the stamp on the dish and gently pushed with a finger to remove air bubbles before applying a 5-6 g weight. We used µ-Dish (35 mm High, uncoated, Nippon Genetics ib81151, ibidi 81151) instead of glass base dishes since they had better passivation properties. After carefully removing the stamp, we applied 2 ml of 0.2 % Pluronic F-127 solution (Molecular Probes P6866) to the entire dish and incubated 60 mins at room temperature. We then washed the dish four times with 1 ml PBS, and lastly poured in the cell solution at 1/4 of confluency and left overnight in the incubator.
The next day, we washed the dish three times using culture media to remove the floating cells, and applied 2 ml of conditioned media with or without the inhibitors before the live imaging.

Calculation of cell flow and density
To calculate cell flow, we first generated the tracks of the individual cells by using the H2B-mCherry signal and the Track-Mate plugin included in the Fiji package [50]. We set a threshold in the detection of cell nuclei so that only the brightest 10-30 % of the cells in the culture were tracked, in order to reduce the probability of mis-tracking. Using the data of tracks, we calculated the local average velocity of flow by averaging in the azimuthal direction (Fig. 2c), y-direction (Fig. 3c, Extended Data Fig. 2), or within the boxes of size 30 µm (Figs. 2b,d,3a). We calculated the single cell velocity (v 0 ) for each condition from the mean-square displacement obtained from the tracks. For the cell density, we used either the counts of the cells in the tracks (Fig. 3b, Extended Data   Fig. 2a,b,c,d) or the sum of the H2B-mCherry signal by subtracting the background (Extended Data Fig. 2e). All analyses were conducted using a custom Python code.

Numerical simulation of the agent-based model
In the numerical simulations (Figs. 3f,g,h,i, Extended Data Fig. 7b,c), we used a model where the cells are represented by agents [51]; the i-th agent at time t is characterized by the twodimensional coordinates r i (t) = (x i (t), y i (t)) and the direction of motion (director) p i (t) = (p ix (t), p iy (t)). Here, p i (t) is a two-dimensional unit vector, |p i (t)| = 1. The time-evolution of r i (t) and p i (t) is defined as where v 0 ,p i , and F ji are the velocity of agents, mean director in the neighbored of the i-th agent, and exclusion force from the j-th agent to the i−th agent, respectively. R η,τ [O] represents the random noise and velocity flipping. By operating R η,τ , a unit vector O = (cos ϕ, sin ϕ) is mapped to R η,τ [O] = (cos(ϕ + η 1 π + η 2 π), sin(ϕ + η 1 π + η 2 π)), where η 1 ∈ [−h 0 , h 0 ] is a uniform random noise, and h 0 (≤ 1) is a parameter characterizing the randomness of the direction. η 2 is another stochastic variable which takes 0 or 1, with the probability of taking 1 given by ∆t/τ, where τ is the velocity flipping rate. The mean director (molecular field) is given bȳ where (p i × p j ) z = p ix p jy − p iy p jx ,ẑ × p j = (−p jy , p jx ) and Note that the summation includes the i-th agent. The first term in Eq. (14) represents the nematic alignment, while the second and third terms represent the chirality-driven alignments parameterized by θ. The chirality-driven alignment is introduced so that the i-th agent nematically interacts with the jth agent with a tilt. To derive the second and third terms in Eq. (14), we consider the rotation of the j-th agent around the z axis (p j → p j = p j − θ∆tẑ × p j ), and plug p j into the nematic alignment term. The autonomous rotation of agents is included in the last term in Eq. (14).
In addition to the chiral nematic alignment, we included the exclusion force F ji to avoid over-clustering of the alignments. We assume the short-ranged repulsive interaction: where x i j = x i − x j , and Θ(x) is the Heaviside step function. k and a are the strength and the length scale of the repulsive force, respectively. We assume that the system is finite along the x-direction (x ∈ [−L/2, L/2]), and is periodic in the y-direction. Near the edges, the cells align their shape with the boundary and can only move parallel to the edge. To mimic this boundary effect, we impose the following condition for the cells at position −L/2 < x ix < −L/2 + 0.50a or L/2 − 0.50a < x ix < L/2: Here, a, and L are the size of the agent (same as the length scale of the repulsive force), and length of the system along the x-direction, respectively. The numerical conditions we used for Figs. 3h,i,4h, Extended Data Figs. 3a,b,7b,c are ρ 0 = 4000 cells/mm 2 , v 0 =40 µm/hour, ∆t =0.02 hours, h 0 =0.06, R 0 =35 µm, k =0.12 µm −1 , a =12.5 µm, with τ=2.5 hours or ∞ (no flipping), θ=0.2 rad/hour or 0 rad/hour (zero chirality). We used a custom C++ code to run the simulation, and made the plots using a custom Python code.

Spectral analysis of image data and simulation
To obtain the power spectrum of the cell density, we used the timelapse image from the H2B-mCherry channel ρ(r, t), where r = (x, y) and t denote the two-dimensional position (in pixels) and time (frame number), respectively. We first took the stripe experiment data with L = 300 µm, and cropped out the 25 µm edge regions from the left and right, and the 100 µm width region from the center. The signal was normalized by subtracting the mean to minimize the spatio-temporal inhomogeneity in the images:ρ(r, t) = ρ (r, t) − t ρ (r, t)/T with ρ (r, t) = ρ(r, t) − r ρ(r, t)/N A , where T and N A are the number of time frames and the number of pixels in the sums, respectively. After taking the sum over x,ρ(y, t) = xρ (r, t), we calculated the discrete Fourier transform ofρ(y, t) to obtain the power spectrum P(k y , ω), and averaged them across multiple regions for each conditions (Figs. 4c,f,g). To obtain the power spectrum of the edge, we assumed 180 • rotational symmetry and flipped the y-axis for the left edge spectrum and added it to the spectrum of the right edge, before the averaging over regions to generate the plots in Fig. 4c,f,g.
To obtain the group velocity from the power spectrum, we first obtained Ω(k y ) = Argmax ω P(k y , ω), and then fit Ω(k y ) with a linear function, Ω(k y ) = v s k y at k y v 0 τ > 5 to obtain the slope v s .
For the Fourier analysis of the simulation, we took data corresponding to the L =300 µm setup, and first calculated the density of the particles along y at each time points at the edge (25 µm wide) and the bulk (100 µm wide) regions ( Fig. 4b) with the binning size of 25 µm in the y-direction. We calculated the power spectrum of these one-dimensional densities by Fourier transform to obtain Extended Data Fig. 7b,c. All analyses were conducted using a custom Python code. numerical calculations have been performed on cluster computers at RIKEN iTHEMS.

HYDRODYNAMIC MODEL AND ITS TOPOLOGY
We derive the hydrodynamic model of the chiral active nematic system using the Boltzmann-Ginzburg-Landau approach [44][45][46]. We start from the Boltzmann equation where f (r, φ, t) is the density of cells at position r = (x, y) and angle φ at time t. The second term in the left hand side describes the contribution where the cells keep moving with speed v 0 in the direction specified by the angle φ with e φ being the unit vector along the angle φ. The third term in the left hand side describes chiral motion of cells, which bends the motion of cells, with strength of the chirality characterized by the parameter θ. The first term in the right hand side accounts for the stochastic flipping of the direction of motion with timescale τ. The second term in the right hand side is the rotational diffusion term with strength D φ . The term I col is the collision integral, for which we only keep contributions that are essential in obtaining the nematic order, as we explain later.
To obtain an effective hydrodynamic equation, we consider the Fourier expansion: and derive an equation governing the first three terms, k = 0, 1, 2, of the Fourier coefficients. Note that since f (r, φ, t) is real, the Fourier coefficients satisfy f −k = f * k . Physically, the k = 0 component describes the density of particles, ρ = f 0 . The k = 1 component is related to the polar vector P = (P x , P y ) through ρP = (Re f 1 , Im f 1 ). The k = 2 component is related to the nematic tensor Q through ρQ xx = Re f 2 /2 and ρQ xy = Im f 2 /2. Expanding the Boltzmann equation in Fourier series and taking k = 0, 1, 2, we obtain where ∆ ≡ ∂ 2 x + ∂ 2 y is the Laplacian, and the diffusion coefficient comes from considering the next order term, f 3 , in the Boltzmann equation and solving for f 3 assuming ∂ t f 3 = 0 and ignoring f 4 and higher order terms. In the same way, we obtain a term involving ν 2 ∆ in the equation for f 2 . When defining l = 10 µm as the minimal length scale (corresponding to the cell size) and normalizing the equation using l, we find that ν 2 is overestimated and can become larger than θ. For the calculations we perform below, we choose the sign of ν 2 to be the same as θ, as expected from the Boltzmann equation, but the value of normalized ν 2 to be smaller than θ. In the equation for f 1 we also have a term proportional to ν 1 ∆. This term does not come from the Boltzmann equation, but we include this term since it is necessary to regularize the theory at small length scale. A similar situation has been analyzed in more detail in [27]. It is natural to take the sign of ν 1 to be the same as θ, by considering the direction of the pressure that a rotating object feels inside the active fluids. These terms, ν 1 and ν 2 , represent non-dissipative forces and torques originating from the chiral motion, where ν 1 is also known as the odd (Hall) viscosity in fluid mechanics.
The final term in the equation for f 2 is derived from the collision integral; when the chirality is absent, this term is responsible for producing the nematic order, For the parameters µ[ρ] and ξ in the term coming from the collision integral, we use microscopic estimates from the Boltzmann-Ginzburg-Landau approach [44,46]: (S6) With this background, we analyze the fluctuations of f 0 , f 1 , and f 2 around the nematically ordered state, which should be achieved under the stripe boundary condition (boundary at x = ±L/2). We assume that the nematic order is built along ydirection, so that in the absence of fluctuations, f 2 = −ρ 0 r < 0. We expand up to first order in fluctuations of density, polar vector, and nematic tensor as f 0 = ρ 0 + δρ, f 1 = ρ 0 (δP x + iδP y ), and f 2 = −ρ 0 r + δ f 2 = −ρ 0 r − rδρ + ρ 0 (δQ xx + iδQ xy ). Keeping up to the second order in these fluctuations and θ and r, we obtain This is the equation whose steady-state and topological properties we are going to analyze. In what follows, we take cell size 10 µm as unit of length, and one hour as the unit of time. This yields v 0 = 4 and ρ 0 = 0.4, which are the parameters adopted in the agent-based simulations and roughly corresponds to the experimental situation. Taking D φ = 0.19, which reproduces experimental results reasonably well as we see below, the strength of the nematic order in our system is r ≈ 0.61.

Steady-state profile
We first analytically explore the chiral edge flow as a steady-state solution of the linearized hydrodynamic equation (S7). We look for a solution which is time-independent and homogeneous along the y-direction, and also assume δP x = 0. Then, the linearlized hydrodynamic equations become: where we set m ≡ µ[ρ 0 ]−8D φ for notational simplicity. Note that Eq. (S9) represents the so-called active current (Eq. (1) in main text). We can find an analytical solution of the set of equations Eq. (S8)-(S11) if we neglect −νr∂ 2 x δρ in Eq. (S10). Therefore, for the sake of analytical derivation, we ignore the term −νr∂ 2 x δρ in Eq. (S10). We have checked that the numerically obtained band structure is largely unaltered by the term −νr∂ 2 x δρ. To solve Eqs. (S8)-(S11), we first write δP y in terms of δQ xy using Eq. (S9). Then, combining Eq. (S8) and Eq. (S10), we can eliminate δρ. Combining the remaining equations, we can finally write a closed differential equation for δQ xx : where , (S13) . (S15) A general solution of this differential equation can be written as where and C 1 , C 2 , C 3 , and C 4 are constants of integral to be determined from boundary conditions. For boundary conditions, we assume δQ xx = δQ xy = 0 at the boundaries x = ±L/2,. We further assume that the integral of the density fluctuation is zero, that is, the number of cells is conserved: With these boundary conditions, we can fully determine the constants of integral. The steady-state solution with these boundary conditions is where the coefficients C 1 and C 2 are and ρ c is a constant to be determined from the condition L/2 −L/2 dxδρ = 0. In Extended Data Fig. 4, we plot the steady-state for θ = 0.2, ν 1 = ν 2 = 0.02 with different L. We see that δP y , which is the average velocity of the steady-state in y-direction, as well as the density fluctuation δρ/ρ 0 is localized at both edges. The steady-state has velocity in positive y-direction at the right edge, and negative direction at the left edge, which is consistent with what we observe in experiment and numerical simulation.
From the steady-state results, we obtain δρ/ρ 0 ∼ L 2 θ 2 /mν, δP y ∼ Lθ 2 /v 0 , and δQ xx ∼ δQ xy ∼ L 2 θ/ν in the lowest order of θ. This corresponds to the fact that the linearization around the nematically ordered state with no flow is justified only when θ is smaller than the inverse of the time scales in the system. When fixing θ and increasing the width of the system L, the condition for the linearization becomes violated, and the nonlinear terms in the hydrodynamic equations should be taken into account. Nevertheless, we here continue using the linear equation for the sake of simplicity, and with the expectation that the basic physics will be unaltered by the nonlinear terms since the steady-state solutions are still qualitatively similar to the experimental and numerical observations. We note that the large L limit of the linear equation is still well-behaved since δρ/ρ 0 , δP y , δQ xx , and δQ xy all converge to O(1) values.

Edge-localized topological Kelvin mode
The steady-state obtained above is the time-independent solution of the linearized hydrodynamic equation. We now show that there are also normal modes of the equation without a source term, which is localized at the edges and propagates unidirectionally. These modes turn out to have a topological origin, and is essentially the Kelvin modes known in geophysics.
As before, we consider a configuration where the system is bounded in the x-direction but long in the y-direction. The momentum k y along y-direction is thus a good "quantum" number, and we can replace ∂ y by ik y in the equation. We look for the solution with δP x = 0, and also neglect, for simplicity, the mixing effects between the momentum and nematic variables. We later compare the full numerical solution of the model with the analytical solution obtained here and confirm that the analytical solution matches the full numerical solution well, validating the approximation we employ here.
The simplified linearized hydrodynamic equation then takes the following form: By solving the linearized equations, we obtain the dispersion relation: We find that the dispersion relation does not depend on θ, and becomes flat at k y ≈ 0, i.e., ReE(k y ) = 0 if |k y | is small. In this momentum range, the waves are overdamped. On the other hand, outside of the range where ReE(k y ) = 0, waves propagate with finite life-time.
Next we discuss the localization of the waves. Since we are interested in the propagating wave, we focus on the momentum range where |k y | is sufficiently large such that the dispersion relation becomes The amplitude of the waves satisfies the differential equation Assuming that ν 1 is small and keeping leading orders in ν 1 , we obtain where D 1 and D 2 are constants to be determined from boundary conditions. From this expression of δP y , we observe that if the sign of θ and ν 1 are the same, the up-going wave [described by E + (k y )] and the down-going wave [described by E − (k y )] are localized to the opposite sides of the edges, so that unidirectional propagating waves appears at the edges. As we see below, this is nothing but the topological edge modes predicted from the bulk-edgecorrespondence, and is what we observed in power spectrum obtained from experiments and numerical simulations (Fig. 4, Extended Data Figs. 5,6,7).
A noticeable feature we can observe here is that the dispersion relation itself does not depend on the strength of the chirality θ, but the localization length does.

Effective Hamiltonian and its topological properties
We now discuss the numerical band structure and topological properties of the linearized hydrodynamic equation. In order to make connection with quantum mechanics and the physics of topological insulators, we multiply both sides of the equation by an imaginary unit i. Defining the equation takes the following form: where withp i ≡ −i∂ i being the momentum operator andp 2 ≡p 2 x +p 2 y . The equation takes the same form as the Scrhödinger equation under an external source s, with Ψ serving as a wavefunction and H as the Hamiltonian. Note that the matrix H is non-Hermitian.

Case of gapped Hamiltonian
We first explore the band structure of the non-Hermitian Hamiltonian H when the strength of the chirality is large and an energy gap opens. We take the y-direction to be long so that the momentum along the y-direction is a good quantum number. In the x-direction, we may take either a periodic boundary condition or an open boundary condition with edges. In Extended Data We observe that in the presence of edges, there are modes crossing the gaps between the bands, which are the topological edge modes. In fact, the Chern number of the lowest two bands of our non-Hermitian model is six, which agrees with the six edge states present in Extended Data Fig. 5c according to the bulk-boundary correspondence. In Extended Data Fig. 5c, we also plot the analytical result for the edge-localized mode obtained in Eq. (S26) as the red lines. We see that one of the six topological edge modes agree well with the analytical line, confirming that this edge state is the topological Kelvin mode. Another mode, which is dispatched from Re(E)τ ≈ ±6 and extends far in momentum, corresponds to the Yanai mode also observed in geophysics [19]. There are four additional topological edge states that are only present when the nematic tensor is included in the equation.
To see how the modes are localized at the edges, we plotted the spectrum with the colours indicating the localization at the edges (Extended Data Fig. 5e,f). We see that all six topological edge modes with positive group velocity are localized at the right edge, whereas those with negative group are localized at the left edge. This observation is again consistent with the bulk-edge correspondence.
Topology of non-Hermitian systems has been recently attracting considerable interest [52][53][54]. In our system, we can continuously transform our non-Hermitian Hamiltonian into a Hermitian Hamiltonian without closing the gap in the real part of the complex energy spectrum. In the terminology of [54], our system is characterized by the symmetry class A with a line gap at ReE(k) ≈ 0, where E(k) is the complex energy as a function of the momentum k. The Hermitian Hamiltonian to which we can continuously deform our non-Hermitian Hamiltonian is In Extended Data Fig. 5g, we plot the real part of the energy spectrum as a function of a parameter which connects H and H Herm ; the left end of the figure corresponds to H and the right end to H Herm . We observe that the energy gap does not close when connecting H and H Herm . We note the essential difference between conventional topological lattice models and our system. Our system, unlike systems defined on a lattice, does not have a compact momentum space. Typically, topological invariant such as the Chern number should be defined on a compact parameter space. The Chern number can be calculated by integrating the Berry curvature over the compact parameter space, and the integral is guaranteed to be an integer. In our system, the momentum space is not compact, so the Chern number, calculated analogously to lattice systems, is not guaranteed to take an integer value. However, it has been shown that when odd viscosity terms (ν 1 and ν 2 ) are present, the Chern number is guaranteed to be an integer [27,55,56]. Indeed, when we calculate the Chern number of our Hamiltonian for the energy gap around the zero energy, we obtain a value which approaches 6 as we take into account larger areas in momentum space for integrating the Berry curvature (Extended Data Fig. 5h). We also note that when θ > 0 and ν 2 < 0, the Chern number of the lowest two bands are two, and there are two topological edge states. When both ν 1 and ν 2 are negative, there is no topological edge state and the Chern number is zero.
As plotted in Extended Data Fig. 4b, the velocity field in the y-direction is localized at the edges. The integrated velocity along the y-direction essentially picks up the k y = 0 component of the eigenstate. Our model, at k y = 0, does not have any topological edge state as one can observe from the band structure. Therefore, we consider that the edge localization of the steady-state itself is a distinct mechanism from the topological edge states analyzed above.
In fact, we notice that many of the modes in the bands in the bulk localize at the edge in the case of the open boundary condition (Extended Data Fig. 5e). This phenomena of the bulk modes localizing in the presence of the boundary is called the non-Hermitian skin effect, and has also been recently found to have a topological explanation [48,49].
Case where the gaps are closed By comparing the experiment and numerical results of the agent-based model, we estimate that the realistic value of chirality is around θ = 0.2. In Extended Data Fig. 6a, we plot the band structure for θ = 0.2 and ν 1 = ν 2 = 0.02. For ν, we used the value calculated from Eq. (S4).
We find that the bulk energy gap is closed, which is due to the large dissipation in the hydrodynamic equation. However, the mode described by Eq. (S26) still exists in the spectrum (Extended Data Fig. 6c), indicating that the Kelvin mode may be surviving even in the case where the gap is closed.
To confirm this, we plotted the spectrum with the colours indicating the localization at the edges and the darkness indicating the contribution to the density fluctuation for each mode (Fig. 4a in the main text). Here, the contribution to the density fluctuation was calculated by normalizing the wave function and looking at the δρ/ρ 0 component. We observe that, among all the modes which exist around zero energy, only the modes forming the large X in the figure survives after extracting modes with significant density fluctuation, which fits with the Kelvin mode [Eq. (S26)]. When lowering the chirality (θ = 0.02, Extended Data Fig. 7a), we found that the similar modes survived but in a non-localized manner. This is consistent with the experiment, where the slope of the dispersion of the power spectrum does not seem to be altered by the change of θ by the application of Jasplakinolide except for the flip in the direction (Fig. 4e in the main text). The existence of the wave modes for the case of small θ, and the change in the localization by increasing θ, were also observed in the spectral analysis of the agent-based model (Extended Data  Fig. 7b,c).   Extended Data Fig. 7. Band structure and localization. a, Energy spectrum of the effective non-Hermitian Hamiltonian calculated with the stripe boundary condition with θ = 0.02 rad/hour. Colour indicates the position of the localization of the modes, and the intensity of the colour indicates the contribution of the mode to the density fluctuation. b, Power spectrum of the agent density calculated for the edge and bulk regions from the numerical simulation with θ = 0.2 rad/hour. c, Power spectrum of the agent density calculated for the edge and bulk regions from the numerical simulation with θ = 0 rad/hour.