Stabilization of active matter by flow-vortex lattices and defect ordering

Active systems, from bacterial suspensions to cellular monolayers, are continuously driven out of equilibrium by local injection of energy from their constituent elements and exhibit turbulent-like and chaotic patterns. Here we demonstrate both theoretically and through numerical simulations, that the crossover between wet active systems, whose behaviour is dominated by hydrodynamics, and dry active matter where any flow is screened, can be achieved by using friction as a control parameter. Moreover, we discover unexpected vortex ordering at this wet–dry crossover. We show that the self organization of vortices into lattices is accompanied by the spatial ordering of topological defects leading to active crystal-like structures. The emergence of vortex lattices, which leads to the positional ordering of topological defects, suggests potential applications in the design and control of active materials.

W hether it is biological matter such as cytoskeletal networks, cellular colonies and suspensions of bacteria or synthetic systems such as Janus catalysts and vibrating granular rods, continuous injection of energy from the constituent elements leads to exotic behaviour such as collective motion [1][2][3] , pattern formation [4][5][6][7] , topological defects [8][9][10] and active turbulence [11][12][13] . Active systems where hydrodynamic interactions are key, such as suspensions of swimming bacteria, are referred to as wet. The energy input to such systems is dissipated by the viscosity of the fluid. However, in many active materials there are alternative forms of dissipation, such as wall friction between the active particles and a substrate 14 . In the limit where the frictional damping dominates, it screens out any velocity fields that are generated by the activity. For instance, the Vicsek model for migrating animal herds 15 , and assemblies of vibrated granular rods 16,17 lack long-range hydrodynamic interactions and are classified as dry active materials. In general, a given active material will fall between these limits with the relative magnitudes of frictional and dissipative damping controlling its position on the wet-dry spectrum. Examples of experimental two-dimensional active nematics where friction could be varied to control the crossover between the wet and dry limits include cells moving and dividing on substrates 18,19 , microswimmers confined between parallel walls 11 and suspensions of microtubules and molecular motors moving at an interface between two fluids 6 .
Here we unify two different classes of active matter by using friction as a control parameter to interpolate between wet and dry active matter, and at the wet-dry crossover we discover an unexpected regime, where otherwise disordered vortices selforganize into lattices interleaved with ordered arrays of topological defects. The vortex lattices and their corresponding network of ordered defects arise from the competition between friction and viscous dissipation and stablize the active system. Our results contribute to understanding the physics of matter operating out of equilibrium, with its potential in the design of active micro-and nano-machines.

Results
Wet-dry crossover. The behaviour of dry and wet active nematics is summarized pictorially at the right-and left-hand sides of Fig. 1. In the absence of momentum conservation, dry active matter is commonly described in terms of the concentration of active components and their orientational order [20][21][22] . The existence of any curvature in the orientational order induces an instability that leads to the formation of bands of concentration 21,22 . The concentration bands are, however, always unstable and eventually break into filaments that, in turn, coalesce and form new bands. The process is repeated and a chaotic regime is established 23 .
Introducing hydrodynamics leads to extra complexity due to the coupling of the density and orientational order of active particles to the fluid flow, and distinct dynamical features are manifest in the models of wet active nematics. In particular, there is a hydrodynamic instability of nematic regions that leads to the formation of walls, lines of high distortions in the director field, where elastic energy is stored. The walls are continuously broken up by the creation and annihilation of topological defects, a process that gives rise once again to unstable nematic regions. The formation and removal of the walls is maintained by the active forcing, and a state termed active turbulence is established, characterized by a chaotic velocity field with regions of high vorticity 9,24 .
Theory. We build on the nematohydrodynamic equations of liquid crystals to describe a natural route from wet to dry active nematics. To this end, the continuum description of the dynamics of passive liquid crystals is modified to account for the active stresses generated by constituent elements 25 . This continuum approach, which allows for coarse graining over the microscopic details, has proven successful in reproducing several experimental observations including the flow behaviour and defect dynamics observed in experiments on microtuble bundles driven by motor proteins 6,9,24 , the spatial organization of bacterial cultures in confined environments 26 , tissue growth 27,28 and fluidization 29 , and the flow fields associated with cell division in cellular monolayers 30 .
The variables needed to describe the hydrodynamics of a wet active nematic are: f the relative concentration of active and passive particles, r the total density, u the velocity vector, and Q ¼ 3q 2 nn À I=3 ð Þ , which is the nematic tensor, with q the magnitude of the orientational order, n the director, and I the identity tensor. The four coupled continuum equations describing the time evolution of these quantities are 31 where the mobility and the rotational diffusivity are denoted by G f and G Q , respectively. The co-rotation term, accounts for the response of the orientational order to the flow gradients characterized by the strain rate tensor E ij ¼ (q i u j þ q j u i )/2 and vorticity tensor O ij ¼ (q j u i À q i u j )/2. The alignment parameter l, which describes the responses of particles to the strain and vorticity takes different values for different particle shapes 32 with l40, lo0 and l ¼ 0 for rod-like, disk-like and spherical particles, respectively. The chemical potential m ¼ dF =df and molecular field represent the relaxation of the concentration and the orientational order to equilibrium, where is the free energy functional and a, A, C, K f , K Q , K c are material constants.
The terms on the left-hand side of equation (4) are the usual terms in the Navier-Stokes equation describing the inertia of the fluid. These are small at low Reynolds number, the limit relevant to active particles on the micron-scale or smaller. Contributions to the stress tensor, P ij , are the viscous dissipation Å viscous ij ¼ 2ZE ij , where Z is viscosity, the elastic stresses that generate backflow and encode mechanical forces on the fluid due to the relaxational motion of the active entities where P denotes the pressure, and, of primary relevance to our current argument, the active stress Å active ij ¼ À zQ ij , which is proportional to the orientational order parameter and distinguishes active from passive nematics 25 . The strength of the activity is set by z. The final term in equation (4) accounts for the friction, for example, with a substrate, with g, the associated friction coefficient. Details of this model in the context of passive and active liquid crystals can be found in refs 31,33-39. In a wet active nematic, the momentum propagation is more effectively suppressed as the friction coefficient is increased. In the limit of sufficiently high friction, the viscous and elastic stresses become small and the friction forces balance with the forces generated by the active term. In this regime the velocity can be approximated in terms of the nematic tensor as In writing equation (8), we are neglecting the fluid inertia, which is dominated by the viscosity for bacterial suspensions and cytoskeletal filaments. We also neglect the elastic terms in the stress tensor that generate backflow. Non-dimensionalizing the momentum equation (see Methods section for more details) shows that the viscous stress dominates the backflow by a factor controlled by the Ericksen number Er ¼ ZUL/K, where U, L are a characteristic velocity and length. Er is B10 2 in our simulations.
We chose parameters that are in a range that has been previously used to reproduce experimental results measuring the velocity correlations of microtuble bundles driven by molecular motors 6,24 and the flow fields of dividing Madine Darby Canine Kidney cells 30 . Moreover, the scaling suggested by equation (8), that velocity increases proportionately with activity, has been observed in both simulations 24 and experiments 6 .
Considering equation (8) as an equality and substituting in the equations (1) and (2) for the concentration and orientational order evolution, we reproduce the standard equations for dry active matter, The second term on the right-hand side of equation (9) is identical to that commonly introduced to represent the current due to curvature in the director field, which drives ordering in dry active nematics. Thus, we demonstrate that the equations of dry active nematics can be generically reproduced from the nematohydrodynamic equations, showing that the dry limit arises naturally when friction dominates viscosity in wet active materials. We use þ y to denote that flow-driven terms that are non-linear in f and Q also appear in equation (9) and (10). These are listed explicitly in the Methods section and their appearance can be expected from symmetry arguments. In addition, we have also used a kinetic approach 21 to show that our analytical derivation of the dry active nematic equations from equations of wet active nematics, matches the microscopic derivation at the kinetic equation level. We show that even higher order terms predicted by our derivation can be expected from the expansion of the kinetic equations to higher order (Methods section).
Simulation. To demonstrate the wet to dry crossover, we present numerical solutions of the active nematic equations. We first illustrate the effect of increasing the friction on the evolution of the concentration and order parameter in a strongly ordered nematic (Fig. 2). Figure 2a shows the time evolution of the concentration field (top row) and director field (bottom row) for an active nematic with an intermediate friction coefficient. As for zero-friction the nematic state is hydrodynamically unstable and walls form in the director field. However, both wall creation, and subsequent destruction by topological defects, are significantly slowed by the friction, which leads to a reduction in the root-mean-square (RMS) velocity by a factor B10 (Fig. 2c). Surprisingly, the slower dynamics is accompanied by an increase in the number of topological defects (Fig. 2d), a consequence of the larger number of walls at higher friction. Similar trends have been recently observed in active nematics without any concentration variation 40 . However, as described earlier, concentration variation is necessary to approach the dry limit where activity manifests in concentration phase separation. Linked to the wall formation is the emergence of coincident concentration bands, with a higher number of active particles at the walls, where the magnitude of the nematic ordering is reduced. The concentration ordering is driven by an advective flux of active particles towards the walls, down the gradient in Q.
In the steady state this is balanced by diffusion of the active particles from high to low concentrations. Figure 2c,d show that the active system behaves very differently for higher values of the friction. The RMS velocity is very low and topological defects are not formed. We note that although, topological defects do not appear in the dry limit, such defects have been observed in shaken rod experiments and in a simulation of active rods. As pointed out in ref. 41 this occurs, in their model at least, because the collisions include a relative rotation of the colliding particles as well as nematic alignment. Figure 2b indicates that concentration bands still appear in the dry limit but, in the absence of an advective flux, the way in which they are formed must be very different. The relevant mechanism is an instability driven by curvature in the nematic order, described by the final term in equation (9). The coupling between the concentration and the nematic order is established by the molecular potential resulting in a strong (weak) ordering at high (low) concentrations. Note the striking difference between the wet and dry limits: in the wet limit, concentration bands are initially formed perpendicular to the director, and the concentration and the nematic order are anti-correlated. In the dry limit, the bands form parallel to the director, and the concentration and nematic fields are correlated. The variation of the RMS velocity suggests the existence of different regimes (Fig. 2c, inset). From equation (8), the magnitude of the velocity can be estimated as u RMS Ezq/gL q , where L q is the characteristic length scale of the nematic order variation. When the frictional damping is strong, L q is set by the distance between the walls, which is approximately proportional to the screening length L sc $ ffiffiffiffiffiffiffiffiffiffi Z=gr p and, therefore, the velocity is expected to vary as u RMS Bg À 1/2 (Fig. 2c, inset: green line). At even higher frictions, there is not enough energy for hydrodynamic wall formation and the active fluxinduced instabilities of dry matter dominate. Therefore, L q does not depend on the screening length and u RMS Bg À 1 (Fig. 2c, inset: blue line).

Flow-vortex lattices and ordered topological defects
We next focus on the dynamical effects of friction at a temperature corresponding to the (passive) isotropic-nematic transition (Fig. 3). In a passive nematic there would be no orientational order. However, order can be stabilized as a result of extensional flows generated by the active forcing for zg40, that is for an extensile, rod-like system (z40, l40) or a contractile, disc-like system (zo0, lo0) 42 .
Considering first extensile rods, for the smallest value of friction considered active turbulence is established. Flowvortices move around chaotically, decay and re-form ( Fig. 3a and Supplementary Movies 1 and 2). As the friction coefficient is increased there is a striking change in behaviour. As is apparent in Fig. 3b, the velocity vortices form ordered, stationary arrays ( Supplementary Movies 3 and 4). Figure 3c indicates how the vortices become narrower at higher friction.
The vortex lattice is interleaved with an ordered network of topological defects (Fig. 3e,f). The defect lattice is directly linked to the velocity pattern: þ 1/2 defects are generated between pairs of counter-rotating vortices due to the change in the sign of the vorticity 24 , while À 1/2 defects are created between two pairs of co-rotating vortices due to distinct flow-induced reorientation of nematic directors (Fig. 3g). The shear flow in the middle of a vortex turns a nematogen towards a 45°orientation relative to the vortex line, while extensional and compressional flows at the edges of a vortex induce vertical and horizontal reorientations, respectively. Similar behaviour is seen for a system of contractile, disk-like particles. In the absence of friction, the activity-induced ordering results in the generation of turbulent-like flow patterns (Fig. 4a,b). As the friction coefficient is increased, vortices arrange to form a lattice (Fig. 4c-e). The long-range order of the vortex lattice is evident from measurements of the vorticityvorticity correlation functions (Fig. 4f), which are defined as , where x is the vorticity. A similar ordered array of vortices has recently been reported in experiments on dividing endothelial cells 19 where the vortex structure was associated with the division-induced flow field.
In the absence of friction, the vorticity length scale is set by the activity and the elastic constant 13,37 . However, for both extensile rod-like and contractile disk-like particles, on increasing friction, the vorticity length scale drops. When it becomes comparable to the screening length vortices are no longer disturbed by the flow and the vortex lattice is established (Fig. 5a).
The positional ordering of defects for extensile and contractile systems can be compared by calculating the structure factor S k ¼ P N d n¼1 e 2pi kÁxn ð Þ , where x n denotes the position of the nth defect 43 . At low friction, positions of topological defects are uncorrelated (Figs 3d and 5b,c; low friction), while positional ordering is observed at high friction values (Figs 3e,f and 5b,c; high friction). The mechanism of the ordering is different in extensile and contractile systems as the primary mode of instability in contractile systems is a splay distortion of the director field, while in extensile systems, bend distortions dominate 44 . This manifests itself in different structure factors of defect ordering for extensile rod-like and contractile disk-like particles, which correspond to rectangular and centred rectangular Bravais crystal structures 43 , respectively (Fig. 5b,c).

Discussion
In addition to the vortex lattice surrounding epithelial cells, a similar positional ordering of defects has been observed in active polar systems in the form of aster arrays [45][46][47] . Moreover, the vorticity pattern we report at high friction resembles selforganized vortices of sperm cells on a surface 48 and periodic vortex arrays reported for a motility assay of microtubules with short-range attractions 7 . Recently, orientational, but not positional, ordering of topological defects has been observed in experiments on microtubule bundles driven by molecular motors 49 , but the mechanism for this is not yet understood.
Friction is present in many active systems and can be tuned, providing mechanistic insights into pattern formation in active materials. For example, in the experiments reported in ref. 6, the microtubule bundles and kinesin molecular motors move on an oil-water interface. The active layer experiences friction due to the momentum transfer to the surrounding fluid and thus the friction can be varied by changing the viscosity difference between the upper and lower fluid layers 50 . In addition, experimental measurements of the flow fields of endothelial cells 19 could be extended by exploring the flow patterns on substrates with different coatings to alter the frictional properties. Indeed recent studies have demonstrated the important role of friction in altering the self-propulsion mechanism of a single cell in confined migration 51 and have shown that the cortical friction can stabilize actin patterns in epithelial tubes 52 . The latter study shows that in the absence of orientational order, the combined effects of friction, actin treadmilling and myosin contractility control the formation of the actin ring on the cell cortex 52 .
The ability to control flow vorticity and defect structure will be important in the design and operation of biological and biomimetic materials. Scientists are just starting to try to create active machines that mimic nature on the microscopic scale 10 and active matter is being considered by those traditionally interested in device physics and novel materials 53 . The formation of the vortex lattice leads to the positional ordering of topological defects that could provide an early step towards the design and control of active materials producing well defined velocity fields. For example, one might envisage that such a vortex lattice has the potential to drive an array of microscopic gears.

Methods
Non-dimensionalization. To characterize the relative importance of different stress contributions, we turn to the non-dimensional form of the momentum equation (equation (4)) introducing the characterisitic length and characterisitc velocity as L, U, respectively, we have x* ¼ xL, u* ¼ uU and p* ¼ p/(rU 2 ), where * denotes dimensionless variables. Equation (4) It shuold be noted that 1/Re where Re ¼ rUL/Z, is the Reynolds number, is the common prefactor for all the stress terms. For the parameters that we use in our simulations, both active stress (zL/ZU) and frictional stress (rgL 2 /Z) are of the same order for large enough friction (4O(0.01)), and much larger than unity. The viscous stresses are on the order of 1/Re and the elastic terms are on the order of 1/(ReEr), where Er ¼ ZUL/K is the Ericksen number. The Ericksen number characterizes the ratio of viscous to elastic forces and is on the order of O(10 2 ) in our simulations. Thus both viscous and elastic stresses are not important for large values of friction coefficient.
The high friction limit. When the energy generated by the active particles is primarily dissipated as friction we may write gu % Àzr Á Q, which gives an expression for the velocity field in terms of distortions in the nematic director field, u % Àðz=gÞr Á Q. Considering this as an equality and substituting into the equation (1) for the concentration evolution gives In this equation, the terms in G f represent diffusive dynamics arising from free energy gradients. The second term on the right-hand side is the current of active particles due to self-generated flow. This is usually introduced as a phenomenological, curvature-driven current in models of dry active nematics 3,21,23 . The third term is a non-linear term arising again from advection of concentration but rarely used in models of dry active systems. Similarly, substituting for the velocity field in the equation (2) for the nematic order parameter we obtain ARTICLE Here, the terms written on the first line are generally used to describe the dynamics of dry active nematics. Terms on the other lines, that are proportional to z/g, are generated by the presence of active flow and flow gradients, just like the curvaturedriven term that appeared in equation (6) in the main text. These contributions are not usually considered in dry active nematics 3,21 , because they represent nonlinear contributions which are assumed to be sub-dominant in defining the dynamics. We now follow the approach of Bertin et al. 21 to show explicitly that such flow coupling terms can arise in general descriptions of dry systems.
Microscopic derivation of the flow-driven terms. We briefly review the microscopic approach used by Bertin et al. 21  however there is no net motion along the orientation vector of each particle, thus imposing nematic symmetry. Considering the motion of the particles, taking into account their angular diffusion and interparticle collisions, Bertin et al. showed that the evolution of the Fourier transform components of the single particle probability distributionf k x; y; t ð Þfollow the master equations whereP k is the Fourier transform of the noise distribution for angular displacements and JP k ; k; q Also, r ¼ @ x þ i@ y and r Ã is the complex conjugate of r.
To close this infinite hierarchy of equations, a small parameter E is introduced that characterizes the distance to the isotropic-nematic transition. Assuming that f k $ jE j k ; @ t $ E 2 ; @ x $ E, the first two modes (k ¼ 0, 1) are given by where the equations are truncated O(e 4 ) and the k ¼ 2 mode is used to close the equations. Coefficients A 0 ; C 0 ; D 0 and E 0 represent material parameters, that may also depend on r. Identifying the zeroth mode as the concentration (f 0 ¼ f) and first mode as the concentration weighted nematic field (f 1 ¼ 2 fQ xx þ ifQ xy À Á ), Bertin et al. 21 obtained the following equations OðE 3 Þ: These are the familiar equations describing dry active nematics, see eg equations (6,7) in the main text or the first lines of equations (12,13) (A minor difference may be attributed to the approaches adopted: as is usual in kinetic descriptions the equations describe the evolution of fQ rather than Q, while continuum approaches are generally not written in terms of concentration weighted fields.) The two additional terms that appear O E 4 ð Þ in the analysis are D 0 r Ã r Ãf 2 1 þ E 0f 1 rrf 1 , the last two terms of equation (16). These terms, that are not used in Bertin et al. 21 have the role of advecting and co-rotating Q with the flow field. They give rise to all the non-linear terms in the last three lines of equation (13) except those proportional to the alignment parameter l. Thus we establish that the flow-driven terms in equation (13) can be present in the kinetic description of dry active systems, but appear only at higher orders.
It is interesting to note that kinetic theory treats the case l ¼ 0, that is, the role of the symmetric part of the velocity gradient tensor is not taken into consideration. Further physics will need to be incorporated in the kinetic approach to capture the flow aligning behaviour of active nematics (terms in the equations of motion that depend on l). Similarly modifying the interactions between the particles may also yield additional coupling and the higher order diffusive terms that appear on the second line of equation (12).
Numerics. The governing equations of active nematohydrodynamics (equations (1-4)) describe a coupled system of partial differential equations. We employ a hybrid lattice Boltzmann technique to solve them. In this method, the equations describing the evolution of density and momentum (equations (3,4)) are solved in tandem using a discretized version of the Boltzmann equation where the first and second moments of the particle distribution function give the density and momentum respectively. The Bhatnagar-Gross-Krook approximation with a single relaxation time is used in the collision operator. Equations (1,2) which respectively describe the evolution of the concentration and orientational order parameter are solved using the method of lines. All spatial differentials are discretized using the second order central difference scheme and time integration is performed using an Euler method. The time step for the method of lines is chosen as 1/10th of that for the lattice Boltzmann updates. The flow field used to evolve the order parameters is updated after every lattice Boltzmann time step while the active and passive stress are determined using the updated order parameter fields in every lattice Boltzmann step, hence ensuring the coupling between the equations in the algorithm. As usual in the lattice Boltzmann technique, discrete space and time steps are chosen as unity. Details can be found in refs 24,31. Simulations are performed on a 100 Â 100 lattice. Unless otherwise stated, the parameters used are G f ¼ 1.0, G Q ¼ 0.1, a ¼ 0, A ¼ À 1.0, C ¼ 0.6, K f ¼ 0.1, K Q ¼ 0.01, K c ¼ 0.0, and Z ¼ 0.05, in lattice units. We use z ¼ 0.008, l ¼ 0.7 for extensile, rod-like particles and z ¼ À 0.02, l ¼ À 0.3 for contractile, disk-like active systems.
It is noteworthy that material parameters in active systems are not known experimentally, and the research in this direction is in progress. In our simulations, we have used the usual lattice Boltzmann units 31 , and therefore an experimentalist may suitably dimensionalize the system of interest by choosing a length, time and mass scale 24,54 .