###### Abstract

We introduce a minimal model for a collection of self-propelled apolar active particles, also called as ‘active nematic’, on a two-dimensional substrate and study the order-disorder transition with the variation of density. The particles interact with their neighbours within the framework of the Lebwohl-Lasher model and move asymmetrically, along their orientation, to unoccupied nearest neighbour lattice sites. At a density lower than the equilibrium isotropic-nematic transition density, the active nematic shows a first order transition from the isotropic state to a banded state. The banded state extends over a range of density, and the scalar order parameter of the system shows a plateau like behaviour, similar to that of the magnetic systems. In the large density limit the active nematic shows a bistable behaviour between a homogeneous ordered state with global ordering and an inhomogeneous mixed state with local ordering. The study of the above phases with density variation is scant and gives significant insight of complex behaviours of many biological systems.

###### pacs:

87.10.Rt, 05.65.+b, 64.60.BdIntroduction :— Active systems are composed of self-propelled particles where each particle extracts energy from its surroundings and dissipates it through motion towards a direction determined by its orientation. These kind of systems are ubiquitous in nature, ranging from very small scale systems inside the cell to larger scales harada (); nedelec (); rauch (); benjacob (); animalgroups (); helbing (); feder (); kuusela31 (); hubbard (), vibrated granular media vnarayan (); kudrolli () etc., and have been studied extensively through experiments, theories and simulations sriramrmp (); tonertusr (); rev (). A collection of head-tail symmetric ‘apolar’ active particles with an average mutual parallel alignment is said to be in a ‘nematic’ state, whereas in an ‘isotropic’ state particles remain randomly oriented. An active system where fluid media do not play important role in emergence of ordered state, and thus the hydrodynamic interactions can be ignored, is called a ‘dry active system’ kemkemer (); vnarayan (); animalgroups (); serra (); schaller (); surrey ().

Active nature of particles introduces a nonequilibrium coupling between density and orientation field, as represented in terms of curvature coupling current in literature sradititoner (); shradhanjop (); sriramrmp (). Such coupling in active nematic induces unusual properties like large density fluctuation sradititoner (); chateprl2006 () and growth kinetics faster than as in usual conserved model shradhatrans2014 (). Recent studies of the active nematic found a defect-ordered nematic state aparnaredner (); shimanatcomm (); yeomans () as opposed to the equilibrium nematic for high particle densities. Recent experiment on amolyiod flibrils ncommam () also found a phase with coexisting aligned and disordered fibril domains, similar to the defect-ordered state obtained in simulations. But few investigations have been done on the behaviours of the active nematic in various density limits, especially at low densities. Here we introduce a minimal model for two-dimensional active nematic and compare various ordering phases of active and equilibrium nematic in different density limits. The ordering in the system is characterised in terms of a scalar order parameter which is the positive eigen value of nematic order parameter Q pgdgenne () in two-dimensions. In the low density limit both active and equilibrium systems are in the isotropic (I) state with particles randomly oriented throughout the whole system (see Fig. 1(b) - I), resulting in a small . The Phase diagram of the active nematic as a function of packing density (see Fig. 1(a)) shows a jump in close to , whereas in the equilibrium case goes continuously to larger values and an isotropic to nematic (I-N) transition occurs close to . In the equilibrium nematic (EN) state particles remain homogeneously oriented in the system (see Fig. 1(b) - EN). At the active system goes from the isotropic to a banded state (BS) where particles cluster and align in the perpendicular direction to the long axis of the band (see Fig. 1(b) - BS-1). With increasing density more number of particles participate in band formation (see Fig. 1(b) - BS-2) and follows a plateau over a range of density. In the large density limit active system shows bistability between a homogeneous ordered (HO) (see Fig. 1(b) - HO) and an inhomogeneous mixed (IM) or local ordered state (see Fig. 1(b) - IM). This IM state is very similar to defect-ordered nematic state in ref. aparnaredner (); shimanatcomm (); yeomans ().

Model :— We consider a two dimensional square lattice. At each vertex ‘’ we define an occupation variable , which can take values (occupied) or (unoccupied), and an orientation variable , which lies between and because of the apolar nature of the particles. Each particle interacts with its nearest neighbours through modified Lebwohl - Lasher Hamiltonian llasher ()

(1) |

where is the interaction strength between two neighbouring particles. This model is analogous to the diluted XY-model with nonmagnetic impurities dilutedxymodel (), where impurities and spins are analogous to vacancies and particles respectively in the present model.

Orientation evolves through Monte - Carlo (MC) updates mcbinder () following the Hamiltonian in Eq. 1. Unlike the diluted XY-model, particles also move on the lattice. Depending on the type of movement we define two kinds of models. (i) ‘Equilibrium model’ (EM) - a particle can diffuse to any unoccupied nearest-neighbouring site, and therefore satisfies the detailed balance condition. (ii) ‘Active model’ (AM) - a particle can move to only those unoccupied nearest-neighbouring sites which are in the direction that makes the least inclination with the particle orientation. Details of the model and particle movement are shown in Supplemental Material SM () section I. The asymmetric move of the active particles does not staisfy the detailed balance condition and arises in general because of the self-propelled nature of the particles in many biological kemkemer (); paxton () and granular systems vnarayan (). These moves produce an active curvature coupling current in coarse-grained hydrodynamic equations of motion shradhanjop (); sradititoner ().

Numerical details :— We consider a collection of particles with random orientation , homogeneously distributed on a square lattice () with periodic boundary condition. The packing density of the system is . We choose a particle randomly and move it to an unoccupied neighbouring site, followed by an orientation updation through MC. We use MC steps to evolve the system to its steady state and all the results have been obtained by averaging over next MC steps. Twenty four realizations have been used for better statistics.

We calculate the scalar order parameter

(2) |

which is small in the isotropic state and close to in the ordered state. First we calculate for EM as a function of inverse temperature for different densities. As shown in Supplemental Material SM () section II, the critical temperature is approximated as . Critical temperature decreases with the lowering of the packing density , similar trend is found in the study of diluted XY-model dilutedxymodel () for varying nonmagnetic site density. In rest of our calculations temperature is kept fixed at and packing density is varied from small values to complete filling .

Phase diagram :— At low densities, , the active system is in the isotropic state where the particles with random orientation remain homogeneously distributed throughout the system, and therefore holds vanishingly small values. The jump occurs in at . For particles cluster in, and both ordered state with high local density and disordered state with low local density coexist (see Fig. 1(b) - BS-1). Mean alignment inside the band is perpendicular to the long axis of the band. In the moderate density limits, band formation in more favourable than lane formation (mean alignment parallel to the long axis of the structure) because the number of particles that can have translational motion is much larger in the banded state, and therefore entropy favours band formation. Similar mechanism create the bending instability close to order-disorder transition in the active polar systems chate (); shradhapre (). The above transition from I state to BS occurs at density lower than the corresponding equilibrium I-N transition density (see Fig. 1(a)).

As we further increase , unlike the equilibrium system where increases monotonically with , the active system shows very small change in for a range of density. This plateau like appearance of with variation in is very similar to plateau phase in magnetization versus field curve of magnetic systems plateauphase (). If an energy gap exists between two consecutive magnetic states, a finite field is required for the magnetic system to go from the lower to the higher state. So until that finite field is applied, the increasing field keeps the system magnetization to be unchanged, and the system is called to be in the plateau phase. With increasing packing density in the plateau regime of the active nematic more particles participate in band formation (see Fig. 1(b) - BS-1 and BS-2). On further increment of density, close to equilibrium I-N transition , transverse fluctuations lead the system to a mixed state shradhanjop (); shimaprl2011 ().

In the large limit active system shows a bistable behaviour with two distinct steady states; first, a state where is large and real space configuration is ‘homogeneous ordered’ (HO), and the second, an ‘inhomogeneous mixed’ (IM) state where realizes some moderate values. In the HO state though the particle orientation is homogeneous, large density inhomogeneity exists in the system (see Fig. 1(b) - HO). IM state is a local ordered state with many aligned clusters of high particle density. The system consists of many such aligned clusters of high density separated from low density disordered regions and mean alignment in each cluster is in different directions (see Fig. 1(b) - IM). IM state is similar to the defect-ordered state recently found in the study of ref. shimanatcomm (); aparnaredner (); yeomans (), with large number of defects in high density active nematic.

Renormalised mean field study for small :— We also calculate the jump in the scalar order parameter and the shift in the transition density using the Renormalised mean field (RMF) method of the coupled coarse-grained hydrodynamic equations of motion for the number density and the order parameter for active nematic shradhanjop (); sradititoner ().

(3) |

and

(4) |

where, is the unit vector along the orientation and is the position of particle . We can obtain the number density by coarse-graining over small subvolume. Eqs. 3 and 4 can be obtained either from symmetry arguments as in ref. sradititoner () or from microscopic rule based model shradhanjop (). Density equation 3 is a continuity equation , because the total number of particles is conserved. Current , where the first term consists of two parts, an active curvature coupling current and anisotropic diffusion , which can also be present in the equilibrium model. The second term in density equation is an isotropic diffusion term. First two terms in the order parameter equation is the alignment term. We choose as a function of density which changes sign at some critical density . Third term is coupling to density and last term is diffusion in order parameter and written for equal elastic constant approximation for two-dimensional nematic.

A homogeneous steady state solution of Eqs. 3 and 4 gives a mean field transition from isotropic to nematic state at density where changes sign. Using RMF method we calculate the effective free energy close to order-disorder transition when is small. We consider density fluctuations and neglect order parameter fluctuations. The effective free energy

(5) |

where , where is a constant. , and . Both and are positive. A detail calculation for is shown in Supplemental Material SM () section III. The density flcutuations introduce a new cubic order term protortional to the activity strength , in the free energy . The presence of such term produces a jump at a lower density . This type of jump and shift in transition because of flcutuations are also called as fluctuation dominated first order phase transition in statistical mechanics coleman () and widely studied in many systems fdfopt (). The jump in and the shift in is proportional to the activity parameter and for we recover the equilibrium transition.

Two-point orientation correlation function :— To further characterise the system we also calculate the two-point orientation correlation at different packing densities, where signifies an average over many realisations. In Fig. 2 we plot vs. on log-log scale, for , , , and for active model and , , and for equilibrium model. For very small packing density , decays exponentially in the active systems. Therefore the system is in short-range-ordered (SRO) isotropic state. At , decays following the power law and therefore the system is in a quasi-long-range-ordered (QLRO) state. At high packing densities correlation functions confirm the bistability in the active systems. For (see Fig. 2(a)) shows power law decay in HO state, whereas in IM state decays abruptly after a distance . The abrupt change in at a certain distance indicates the presence of local ordered clusters. In contrast, the equilibrium systems show a transition from SRO isotropic state at low density to QLRO nematic state at , similar to Berezinskii - Kosterlitz - Thouless (BKT) transition Berezinskii (); KT () in the diluted XY-model dilutedxymodel ().

Orientation distribution :— We also compare the steady state properties of active and equilibrium models in the high density limit. First we calculate the steady state (static) orientation distribution from one snapshot of particle orientation. Both active HO and equilibrium nematic show a Gaussian distribution of orientation (see Fig. 3(a)). Hence in HO state orientation fluctuations of particles are of same kind as in equilibrium model and the system is in QLRO state. Distribution in the IM state is very broad and spanning the whole range of orientation. Hence the system has many local ordered regions with different orientations.

We also calculate the time averaged distribution of mean orientation of all the particles in active HO and equilibrium nematic states. is calculated from mean of all particle orientaions averaging over a long time (from to ) in the steady state. for HO is narrow in comparison to the broad distribution for equilibrium model (see Fig. 3(b)). Narrow distribution of implies that orientation autocorrelation in active system decay slowly as comapared to the corresponding equilibrium model which is in agreement with the slow decay of velocity autocorrelation in bacterial suspension wulib ().

Summary :— In this letter we have introduced a minimal model for active nematic and found three distinct phases with the variation in density. At low densities the active nematic is in disordered isotropic state with very small correlation between the particles. With increasing density active nematic undergoes a fluctuation induced first order phase transition from the isotropic to the banded state where large number of particles participate in band formation. Large density fluctuations in the active systems change the nature of the transition and shift the transition density to smaller value as compared to the equilibrium isotropic nematic transition. At large densities equilibrium nematic is in QLRO nematic state, whereas active nematic goes from the banded state to either the homogeneous ordered (high ) or the inhomogeneous mixed (moderate ) state. This inhomogeneous mixed state is similar to the phase with coexisting aligned and disordered fibril domains found in recent experiment ncommam (). Experiments on thin layer of agitated granular rods, elongated living cells, bacterial colony of apolar B. subtilis etc. at different densities can realize the different phases we found here. In the present model we have frozen the motion of the active particles in the transverse direction, i.e. the activity strength is kept large. It will be interesting to see the evolution of different phases with the particles having a small probability to move in transverse directions as well.

###### Acknowledgements.

S. M. acknowledges Thomas Niedermayer for useful discussions. S. M. and M. K. acknowledges financial support from the Department of Science and Technology, India, under INSPIRE award 2012 and Ramanujan Fellowship respectively.## References

- (1) Y. Harada, A. Nogushi, A. Kishino and T. Yanagida, Nature (London) 326, 805 (1987); M. Badoual, F. JÃ¼licher and J. Prost, Proc. Natl. Acad. Sci. U.S.A. 99, 6696 (2002).
- (2) F. J. NÃ©dÃ©lec, T. Surrey, A. C. Maggs and S. Leibler, Nature (London) 389, 305 (1997).
- (3) E. Rauch, M. Millonas and D. Chialvo, Phys. Lett. A 207, 185 (1995).
- (4) E. Ben-Jacob, I. Cohen, O. Shochet, A. Tenenbaum, A. CzirÃ³k and T. Vicsek, Phys. Rev. Lett. 75, 2899 (1995).
- (5) Animal groups in Three Dimensions, edited by J. K. Parrish and W. M. Hamner (Cambridge University Press, Cambridge, 1997).
- (6) D. Helbing, I. Farkas and T. Vicsek, Nature 407, 487 (2000); D. Helbing, I. J. Farkas, and T. Vicsek, Phys. Rev. Lett. 84, 1240 (2000).
- (7) T. Feder, Phys. Today 60(10), 28 (2007); C. Feare, The Starlings (Oxford University Press, Oxford, 1984).
- (8) E. Kuusela, J. M. Lahtinen and T. Ala-Nissila, Phys. Rev. Lett. 90, 094502 (2003).
- (9) S. Hubbard, P. Babak, S. Sigurdsson and K. Magnusson, Ecol. Modell. 174, 359 (2004).
- (10) V. Narayan, N. Menon and S. Ramaswamy, J. Stat. Mech. P01005 (2006); V. Narayan, S. Ramaswamy, and N. Menon, Science 317, 105 (2007).
- (11) D.L. Blair, T. Neicu and A. Kudrolli, Phys. Rev. E 67, 031303 (2003).
- (12) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143 (2013).
- (13) J. Toner, Y. Tu, and S. Ramaswamy, Ann. Phys. (Amsterdam) 318, 170 (2005).
- (14) S. Ramaswamy, Annu. Rev. Condens. Matter Phys. 1, 323 (2010) .
- (15) R. Kemkemer, D. Kling, D. Kaufmann and H. Gruler, Eur. Phys. J. E 1, 215, (2000)
- (16) X. Serra-Picamal, V. Conte, R. Vincent, E. Anon, D. T. Tambe, E. Bazellieres, J. P. Butler, J. J. Fredberg and X. Trepat, Nat. Phys. 8, 628 (2012).
- (17) V. Schaller, C. Weber, C. Semmrich, E. Frey and A. R. Bausch, Nature (London) 467, 73 (2010).
- (18) T. Surrey, F.âJ. NÃ©dÃ©lec, S. Leibler and E. Karsenti, Science 292, 1167 (2001).
- (19) S. Ramaswamy, R. A. Simha and J. Toner, Europhys. Lett. 62, 196 (2003).
- (20) E. Bertin, H. Chaté, F. Ginelli, S. Mishra, A. Peshkov and S. Ramaswamy, New J. of Phys. 15, 085032 (2013).
- (21) H. Chaté, F. Ginelli, and R. Montagne Phys. Rev. Lett. 96, 180602 (2006).
- (22) S. Mishra, S. Puri, and S. Ramaswamy, Phil. Trans. of the Royal Soc. A: Mathematical, Physical and Engineering Sciences 372, 20130364 (2014).
- (23) E. Putzig, G. S. Redner, A. Baskaran, A. Baskaran, arXiv:1506.03501 (2015).
- (24) Xia-qing Shi and Yu-qiang Ma, Nat. Comm. 4, 3013 (2013).
- (25) S. P. Thampi, R. Golestanian and J. M. Yeomans, Euro. Phys. Lett. 105, 18001 (2014).
- (26) S. Jordens, L. Isa, I. Usov and R. Mezzenga, Nat. Comm. 4, 1917 (2013).
- (27) P. G. de Gennes and J. Prost, The Physics of Liquid Crystals (Clarendon Press, Oxford, 1995).
- (28) âA. Lebwohl and G. Lasher, Phys. Rev. A 6, 426 (1972).
- (29) Y. M. Blanter, Y. E. Lozovik and A. Y. Morozov, Phys. Scr. 52, 237 (1995); S. A. Leonel, P. Z. Coura, A. R. Pereira, L. A. S. MÃ³l, and B. V. Costa, Phys. Rev. B 67, 104426 (2003).
- (30) D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, Second edition (Cambridge University Press, Cambridge, 2005).
- (31) Supplementary Material for (I) model figure, (II) estimate of critical temperature in equilibrium model and (III) renormalised mean field calculation of effective free energy for small .
- (32) W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. St. Angelo, Y. Cao, T. E. Mallouk, P. E. Lammert and V. H. Crespi, J. Am. Chem. Soc. 126, 13424 (2004).
- (33) H. Chaté, F. Ginelli, G. Grégoire and F. Raynaud, Phys. Rev. E 77, 046113 (2008).
- (34) S. Mishra, A. Baskaran and M.âC. Marchetti, Phys. Rev. E 81, 061916 (2010).
- (35) M. Takigawa and F. Mila, Introduction to Frustrated Magnetism: Materials, Experiments, Theory (Springer, Heidelberg, 2011), ch 10
- (36) Xia-qing Shi, Yu-qiang Ma, arXiv:1011.5408 (2010).
- (37) S. Coleman and E. Weinberg, Phys. Rev. D 7, 1888 (1973).
- (38) B. I. Halperin, T. C. Lubensky and S. K. Ma, Phys. Rev. Lett. 32, 292 (1974); J. H. Chen, T. C. Lubensky and D. R. Nelson, Phys. Rev. B 17, 4274 (1978).
- (39) V. L. Berezinskii, Sov. Phys. JETP 34(3), 610 (1972).
- (40) M. Kosterlitz and D. Thouless, J. Phys. C 6, 1181 (1973).
- (41) Xiao-Lun Wu and A. Libchaber, Phys. Rev. Lett., 84, 3017 (2000).

Supplemental Material

## Appendix A I. Model figure

## Appendix B II. Estimate of critical temprature in equilibrium model

## Appendix C III. Renormalised mean field (RMF) study of active nematic for small scalar order parameter

In this section we will write an effective renormalised mean field free energy
for scalar order parameter for small . We keep the density fluctuations
and ignore the order parameter fluctuations in the coupled
hydrodynamic equations of motion for active nematic.
Density fluctuation produces a cubic order term in the effective
free energy for scalar order parameter and such term produces a jump
in at a new transition density lower than
equilibrium I-N transition point .
Shift in transition density and jump is directly
proportional to the activity
parameter and we recover equilibrium limit for zero .

We write coupled hydrodynamic equations of motion for density and
order parameter
where nematic order parameter pgdgenne ()

(6) |

being the coarse grained angle at position and time . Density equation

(7) |

and order parameter equation

(8) |

Density Eq. 7 is a continuity equation , where has two parts, active and diffusive. Details of these two currents are given in the main text. is the activity parameter, present because of self-propelled nature of the particles, is the coupling of density in equation. and are the diffusion coefficients in density and order parameter equations respectively, and are the alignment terms and ingeneral depends on the model parameters. For metric distance interacting models njopshradha () is a function of density and changes sign at some critical density. We choose and . Steady state solution of homogeneous Eq. 7 is , we add small perturbation to mean density . In the staedy state density fluctuation can be obtained from Eq. 7,

(9) |

where and and keep the lowest order terms in and

(10) |

and

(11) |

Here we assume nematic is aligned along one direction and there is variation only along the perpendicular direction. Hence we can choose either of equations 10 or 11. Two constants and are the value of density where nematic order parameter is zero.

We use Eq. 10 and substitute the solution for density in equation for and obtain an effective equation for .

(12) |

We neglect all the derivative terms and keep only polynomial in , i.e. we neglect higher order fluctuations. We can do taylor expansion of about mean density , , where . Substitute the expression for from Eq. 10 hence

(13) |

We can write an effective free energy for

(14) |

hence

(15) |

Therefore

(16) |

where , and and is a conatant. Hence fluctuation in density introduces a cubic order term in effective free energy . Effective free energy in Eq. 16 is similar to Landau free energy with cubic order term chaiklub (). We calculate jump and new critical density from coexistence condition for free energy. Steady state solutions of order parameter ( for isotropic and for ordered state) are given by

(17) |

Non-zero is given by . Coexistence condition implies

(18) |

hence we get the solution

(19) |

and

(20) |

Hence the jump at new critical point is . Since and hence , the new critical density

(21) |

is shifted to lower density in comparison to equilibrium . Eq. 21 gives the expression for new transition density as given in main text. Hence using renormalised mean field theory we find a jump at a lower density as compared to equilibrium I-N transition density.

## References

- (42) P. G. de Gennes and J. Prost, The Physics of Liquid Crystals (Clarendon Press, Oxford, 1995).
- (43) E. Bertin, H. Chaté, F. Ginelli, S. Mishra, A. Peshkov and S. Ramaswamy, New J. of Phys. 15, 085032 (2013).
- (44) P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics (Cambridge University Press, Cambridge, 2000).