Order-disorder transition in active nematic: A lattice model study

We introduce a lattice model for active nematic composed of self-propelled apolar particles, study its different ordering states in the density-temperature parameter space, and compare with the corresponding equilibrium model. The active particles interact with their neighbours within the framework of the Lebwohl-Lasher model, and move anisotropically along their orientation to an unoccupied nearest neighbour lattice site. An interplay of the activity, thermal fluctuations and density gives rise distinct states in the system. For a fixed temperature, the active nematic shows a disordered isotropic state, a locally ordered inhomogeneous mixed state, and bistability between the inhomogeneous mixed and a homogeneous globally ordered state in different density regime. In the low temperature regime, the isotropic to the inhomogeneous mixed state transition occurs with a jump in the order parameter at a density less than the corresponding equilibrium disorder-order transition density. Our analytical calculations justify the shift in the transition density and the jump in the order parameter. We construct the phase diagram of the active nematic in the density-temperature plane.


INTRODUCTION
Self-propelled particles compose an interesting type of the active systems [1][2][3][4][5] where each particle extracts energy from its surroundings and dissipates it through motion and collision.Their examples range from very small intracellular scale to larger scales [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20].Also many artificially designed systems, e.g., vibrated granular media [21][22][23][24], active polar disks [25], active colloids [26][27][28][29] imitate the physics of the active systems.If n is the average alignment direction of a collection of such active particles, and the system remains invariant under the transformation n → −n, it is called 'active nematic'.Activity introduces many interesting properties which are absent in their thermal equilibrium counterparts.One of such interesting features is the presence of large density fluctuation in the ordered active nematic [30].Density is a key control parameter in various experiments and numerical simulations.Earlier studies on equilibrium nematic for a fixed temperature show the isotropic to nematic transition at some critical density [31].However, the effect of the density fluctuation in the active nematic is not well understood.
Most of the previous studies of the active systems are done either by using the coarse-grained hydrodynamic equations of motion [32] or microscopic rule based numerical simulation of agent based point particles [33] or Brownian dynamic simulation [34].Here we introduce a lattice model for a two-dimensional active nematic, explore various states of the system in the density-temperature plane, and compare it with the corresponding equilibrium model.In general, lattice model itself is interesting for development of simplified theories, and provides insight into complex systems.Our model is analogous to the previous lattice model of polar active spins [35,36]; but we include volume exclusion to avoid multiple occupancy on single site.Such volume exclusion limits the motion of particles towards an occupied neighbouring site, and introduces new features, e.g., typical pattern formation [20,37], density induced motility [38] in the system.
We construct a phase diagram for the active nematic in the density-temperature plane, as shown in Fig. 1(a).There we observe -(i) disordered isotropic (I) state in low density regime, (ii) locally ordered inhomogeneous mixed (IM) state in intermediate density regime, and (iii) bistability between the IM and a homogeneous globally ordered (HO) state in high density regime.In contrast to the continuous isotropic to nematic (I-N) transition in the equilibrium system, the I to IM state transition in the active nematic in the low temperature regime occurs with a jump in order parameter, as shown in Fig. 2(a).This transition occurs at a density lower than the equilibrium critical value, and the system forms clear bands (BS) in this regime.We finally justify the jump in the order parameter and the shift in the transition density by analytical study of the coarse-grained hydrodynamic equations written for the active model.

MODEL
We consider a collection of apolar particles on a two dimensional square lattice, as shown in schematic diagram Fig. 3(a).Occupation number 'n i ' of the i th lattice site can take values 1 (occupied) or 0 (unoccupied).Orientation θ i of apolar particle at the i th site can take any value between 0 and π.The model follows two sequential processes at every step; first, a particle moves to a nearest neighbouring site with some probability, and then orientation of the particle is updated based on its nematic interaction with its nearest neighbours.We define two kinds of models on the basis of particle movement: (i) 'Equilibrium model' (EM) -particle moves with equal probability 1/4 to any of the four neighbouring sites (Fig. 3(b)), (ii) 'Active model' (AM) -in this model particle movement occurs in two steps.First, it chooses a direction along which it is more inclined.As shown in Fig. 3(c,d), it chooses the direction of movement along BD if π/4 < θ ≤ 3π/4 and along AC otherwise.In the second step, it moves to a randomly selected site between the two nearest neighbouring sites along the chosen direction.For example, if BD is selected as the direction of movement, then the particle moves to randomly selected site B or D in the second step.In both the models, we consider volume exclusion, i.e., particle movement is allowed only if the selected site is unoccupied.
In both the models, the particles also interact with their nearest neighbours.The interaction depends on the relative orientation of the particles and is represented by a modified Lebwohl-Lasher Hamiltonian [39] where is the interaction strength between two neighbouring particles.The interaction in equation (1) governs the orientation update of the particle.We employ Metropolis Monte-Carlo (MC) algorithm [40] for orientation update of the particle after the movement trial.In both the models, an order parameter defining the global alignment of the system does not remain conserved during the MC orientation update described above.In actual granular or biological systems where mutual alignment emerges because of steric repulsion, orientation of particles need not to follow a conservation law.An order parameter defined by coarse-graining the orientation in our present model is a class of non-conserved order parameter: Model A as described by Hohenberg and Halperin [41].
Both the models EM and AM comprise of two different physical aspects -motion of the particles and nematic interaction amongst the nearest neighbours.If the particles are not allowed to move, the models reduce to an apolar analogue of the diluted XY-model with nonmagnetic impurities [42], where impurities and spins are analogous to vacancies and particles, respectively.However, unlike the diluted XY-model, particles in these models are dynamic.In the EM, the particle diffuses to neighbouring sites, whereas it moves anisotropically in the AM.The anisotropic movement of the active particles arises in general because of the self-propelled nature of the particles in many biological [43] and granular systems [21,22].This move produces an active curvature coupling current in the coarse-grained hydrodynamic equations of motion [30,32].The AM does not satisfy the detailed balance principle [40], because of the orientation update after the anisotropic movement.The coupling of the particle movement with the orientation update in our active model is analogous to the active Ising spin model introduced by Solon and Tailleur [35,36], where the probabilistic flip of the spins is an equilibrium process, whereas the out-of-equilibrium aspect of the model is attributed to the anisotropic movement probability of the spins.However, their orientation update algorithm [35,36] is similar to kinetic Monte-Carlo, whereas we use Metropolis Monte-Carlo algorithm to update particle orientation.

NUMERICAL STUDY
We consider a collection of N particles with random orientation θ ∈ [0, π] homogeneously distributed on a L × L lattice (L = 256, 400, 512) with periodic boundary.Packing density of the system is defined as C = N/(L × L).We choose a particle randomly, move it to a neighbouring site obeying exclusion, and then update its orientation using Metropolis algorithm.In each iteration, we repeat the same process for N number of times, and we use 1.5 × 10 6 iterations to achieve the steady state of the system.We obtain the steady state results by averaging the observables over next 1.5 × 10 6 iterations and use more than twenty realisations for better statistics.
The ordering in the system is characterised by a scalar order parameter defined as It is proportional to the positive eigenvalue of the nematic order parameter Q [31].It takes the minimum value 0 in the disordered state and the maximum value 1 in the complete ordered state.First we study the EM as a function of inverse temperature β = 1/k B T for different packing densities.As shown in Fig. A1, the system shows disordered isotropic to nematic state (I-N) transition with decreasing temperature.In contrast to the first order I-N transition in the equilibrium Lebwohl-Lasher model in three dimensions [39,44], we find continuous transition for the EM defined in two dimensions.The observed nature of transition supports the study by Mondal and Roy [45].Similar to the diluted XY-model [42], the critical inverse temperature β c (C) increases with density in the EM.

Phase diagram
We construct phase diagram for both the equilibrium model and the active model on the density-temperature plane.
As shown in Fig. 1(a), two distinct states appear in the EM -(i) an equilibrium isotropic (EI) state on the left side of the red boundary and (ii) an equilibrium nematic (EN) state on the right side of the red boundary.In the EI state, particles remain disordered and homogeneously distributed throughout the system.Consequently, the scalar order parameter S 0 in this state.With increasing density or decreasing temperature the particles get mutually ordered and form the EN state (S > 0).As shown in Fig. 2(a), for a fixed temperature the scalar order parameter increases continuously with increasing density, and the system enters into the nematic state.Both the particle orientation and the coarse-grained density remain homogeneous in the EN state, as shown in the real space snapshot Fig. 1(b).Similar to the EM, the active system remains in a homogeneous disordered isotropic (I) state in the high temperature and/or low packing density regime (cyan coloured regime in the phase diagram Fig. 1(a)).With increasing density or decreasing temperature, beyond the I state, the active system enters into an inhomogeneous mixed (IM) state (golden regime in the phase diagram Fig. 1(a)), where locally ordered high-density domains coexist with disordered low-density regions.In the low temperature regime (β ∈ [1.9, 2.2]), the I to IM state transition with increasing C occurs with a jump in the scalar order parameter S, as shown in Fig. 2(a).In the very beginning of the IM state, as indicated by cross symbols in Fig. 1(a), we find a banded state (BS) in the low temperature regime, where particles cluster and align themselves within a strip to form band.However, out of the strip the system remains disordered with low local density, as shown in the real space snapshot Fig. 1(b).On further increment of the packing density C, bands formed in different directions start mixing leaving the system with many locally ordered high density patches separated by low density disordered regions.Typical real space snapshots for the orientation and the coarse grained density in the IM state are shown in Fig. 1(b).The jump in the S − C curve reduces with increasing temperature, and no bands appear in the high temperature (β < 1.9) regime.shows that the I to BS transition occurs in the low temperature regime with a jump in S at a density lower than the corresponding equilibrium I-N transition density C IN .These bands appear because of the large activity strength.A linear stability analysis, as detailed later in this paper, shows that the large activity strength induces an instability in the disordered isotropic state.This instability goes away for small activity strength or at high temperature.We also do a renormalised mean field calculation of an effective free energy written for the active nematic.The calculation predicts a jump in the scalar order parameter and shows a shift in the disordered (S = 0) to ordered (S = 0) state transition density.Both the jump in S and the shift in the transition density reduce with the activity strength or increasing temperature.The I to BS transition is a first order transition.The shift in the disorder-order transition point is a common feature of the active systems.For large activity and low temperature, if the system density is above a certain value but less than C IN , the large density fluctuation present in these systems causes local alignment with local density higher than C IN .Large density fluctuation is an intrinsic feature of the active systems, and as shown in Fig. 4, we also observe the same in the ordered active states in our model.Due to activity these locally ordered regions move anisotropically and combine with nearby region with similar local ordering.So larger ordered region forms at mean density lower than the equilibrium I-N transition density.Therefore, we find a disordered to ordered state transition at a lower density.For large activity strength, I-BS transition occurs with the jump in scalar order parameter.In our numerical study we calculate the probability P (S) of the scalar order parameter averaging over many iterations and realisations near the I-BS transition point.Figure 2(b) shows P (S) has two peaks, which further supports the first order I-BS transition for large activity strength.
In the high density regime (red coloured regime in the phase diagram Fig. 1(a)), the AM shows bistability, i.e., it can be either in the locally ordered IM state or in a homogeneous globally ordered (HO) state.As shown in Fig. 2(a), the S − C curve for fixed temperature bifurcates in the high density regime; the lower branch corresponds to the earlier discussed IM state, whereas the higher branch indicates the existence of the globally ordered state.Figure 1(b) shows that the system possesses less density inhomogeneity in the HO state compared to the IM state.A finite size scaling of both the HO and the IM state, as shown in the Fig. 2(c), shows that the active nematic possesses non-zero finite order in both these states.Order parameter time series shown in Fig. 2(d) confirms the bistability of the system in the high density regime.Bistability is not generally seen in other agent based numerical simulations of point particles [33]; it appears because of finite filling constraint of the model.This feature can be suppressed if we allow more than one particle to sit together.In the complete filling limit C = 1.0, the AM is equivalent to the EM, and it shows the globally ordered HO state only.

Two-point orientation correlation
We further characterise various states on the basis of the two-point orientation correlation in the different states of the equilibrium and the active nematic.It is defined as g 2 (r) =< i n i n i+r cos [2 (θ i − θ i+r )] / i n 2 i > where r represents interparticle distance, and < .> signifies an average over many realisations.Figure 5(a, b) show g 2 (r) versus r plots on log-log scale for the AM and the EM, respectively, for a fixed inverse temperature β = 2.0.In the AM, g 2 (r) decays exponentially at low packing density C < 0.37, i.e., in the isotropic state.Therefore, the active isotropic is a short-range-ordered (SRO) state.In the BS at C = 0.38, g 2 (r) decays following a power law.Therefore, the system is in a quasi-long-range-ordered (QLRO) state.Ordering increases with density.At high packing density, correlation function confirms the bistability in the active system.At C = 0.82, g 2 (r) shows power law decay in the HO state, whereas in the IM state g 2 (r) decays abruptly after a distance r.The abrupt change in g 2 (r) at a certain distance indicates the presence of locally ordered clusters in the IM state.In contrast, the equilibrium system shows a transition from SRO (exponential decay) isotropic state at low density C < ∼ 0.48 to QLRO (power law decay) nematic state at high density C > ∼ 0.50.

Orientation distribution and autocorrelation of the mean orientation
We compare the steady state properties of the active and the equilibrium models in the high density limit.First we calculate the steady state (static) orientation distribution P (θ) from a snapshot of particle orientation θ.As shown in Fig. 6(a), both the active HO and the equilibrium nematic show Gaussian distribution of orientation.Peak position of P (θ) for both the EN and the HO state can appear at any point between 0 and π because of the continuous broken rotational symmetry of the Hamiltonian shown in equation (1).Data shown in Fig. 6(a) is for one realisation only, and for other realisations also the distribution P (θ) remains Gaussian with peak at other θ values.Therefore, orientation fluctuation of the particles in the active HO state is same as in the equilibrium nematic state.The distribution P (θ) in the IM state is very broad and spans over the whole range of orientation.Therefore, the system possess no global ordering in the IM state.
We also calculate the time averaged distribution P ( θ) of mean orientation of all the particles in the active HO and the equilibrium nematic states.The mean orientation θ(t) of all particles is calculated for each iteration time t in the steady state.The distribution P ( θ) of the mean orientation is obtained from these θ(t) data.This distribution is a measure of the fluctuation in the global orientation of the particles in the steady state.As shown in Fig. 6(b), P ( θ) in the active HO state is narrow in comparison to the broad distribution in the EN state.We also calculate the autocorrelation of the mean orientation Cθ(t) =< 1 t t τ =1 cos 2 θ (t 0 ) − θ (t 0 + τ ) > in the steady state.As shown in Fig. 6(c), Cθ(t) decreases with time in the EN state, but remains unchanged in the active HO state.Both these results imply that the fluctuation in the global orientation direction θ in the active HO state is small compared to the EN state.We do not calculate the mean orientation θ in the active IM state, because the system possesses no global ordering in this state.

PHENOMENOLOGICAL APPROACH TO UNDERSTAND LOW DENSITY STATES OF THE ACTIVE MODEL
In this section we write the hydrodynamic equations of motion for the active model and characterise the low density states of the system.The equations of motion for the slow variables of the system, i.e., the number density ℘(r, t) = l δ(r − R l (t)) and the order parameter ) are as follows [30,32]: and Here R l (t) represents position of the particle l, and m l = (cos θ l , sin θ l ) is the unit vector along the orientation θ l .The total number of particles being a conserved quantity of the system, equation (3) represents a continuity equation The first term of J i consists of two parts: an anisotropic diffusion current J p1 ∝ Q ij ∇ i ℘ and an active curvature coupling current J a ∝ a 0 ℘∇ j Q ij where a 0 is the activity strength of the system.The second term represents an isotropic diffusion J p2 ∝ ∇℘.The α terms in equation ( 4) represent mean field alignment in the system.We choose α 1 (℘) = ( ℘ ℘ IN − 1) as a function of density that changes sign at some critical density ℘ IN .The β term represents coupling with density.The last term represents diffusion in order parameter that is written under equal elastic constant approximation for two-dimensional nematic.The steady state solution ℘(r, t) = ℘ 0 and Π(r, t) = Π 0 , where Π 0 = α1(℘0) α2 , of equations ( 3) and ( 4) represents a homogeneous ordered state for α 1 (℘ 0 ) > 0 at ℘ 0 > ℘ IN , and a disordered isotropic state for α 1 (℘ 0 ) < 0 at ℘ 0 < ℘ IN .
We study the linear stability of the disordered isotropic state (Π 0 = 0) by examining the dynamics of spatially inhomogeneous fluctuations δ℘(r, t) = ℘(r, t) − ℘ 0 , δΠ 11 = Π 11 (r, t), and δΠ 12 = Π 12 (r, t).We obtain the linearised coupled equations of motion for small fluctuations as Using Fourier transformation we get linear set of equations in the Fourier space as where M is the coefficient matrix as obtained from equations ( 5), (6), and (7) after the transformation.We solve equation ( 9) for the hydrodynamic modes λ.We choose q x = q y = q √ 2 since both the directions are equivalent.Therefore, we obtain For small wave-vector q, we can find an unstable mode For small D ℘ and large actvitity a 0 this mode becomes unstable for q < q c , where provided ∆D = D Π −D ℘ is positive and a 0 β > 8D ℘ ∆D.Therefore, the unstable mode λ + causes the I -BS transition for small diffusivity, i.e., at low temperature, and for large activity strength a 0 .We also calculate the jump in the scalar order parameter S and the shift in the transition density from equations (3) and (4).A homogeneous steady state solution of these equations gives a mean field transition from the isotropic to the nematic state at density ℘ IN where α 1 (℘) changes sign.Using renormalised mean field (RMF) method we calculate an effective free energy F ef f (S) close to the order-disorder transition where S is small.We consider density fluctuations δ℘ and neglect order parameter fluctuations.The effective free energy is where b 2 = α 1 (℘) + α 1 c, where c is a constant.[46] and are widely studied for many systems [47,48].The jump in S and the shift in the transition density are proportional to the activity strength a 0 , and for a 0 = 0 we recover the equilibrium transition.

DISCUSSION
In our present work we have introduced a minimal lattice model for the active nematic and study different ordering states in the density-temperature plane.A brief summary of the results is as follows.In the low density regime, the system is in the disordered isotropic (I) state with short range orientation correlation amongst the particles.In the low temperature regime, large density fluctuation in the active system induces a first order transition from the isotropic to the banded state with a jump in the scalar order parameter at a density lower than the equilibrium isotropic-nematic (I-N) transition density.The linear stability analysis of the isotropic state shows an instability for large activity strength in the low temperature regime.Such instability governs the band formation at density below the equilibrium I-N transition density.As we further increase density, bands vanish and locally ordered patches appear in the inhomogeneous mixed (IM) state.Renormalised mean field calculation confirms the jump in the scalar order parameter and the shift in the transition density.With increasing temperature the shift in the transition density and the jump in scalar order parameter decreases, and no bands appear in the system.The IM state is a state with coexisting aligned and disordered domains, similar to the coexisting or defect-ordered states found in Ref. [33,34,[49][50][51][52][53].
In the high density regime, the active nematic shows switching between the IM (low S) and the homogeneous ordered (HO, high S) states, i.e., the system shows bistability.In the complete filling limit and with excluded volume assumption the active model reduces to the equilibrium model.Therefore, the active model tends to show a homogeneous nematic state in the high density regime.However, large activity strength makes the HO state unstable and leads the system to the IM state.This instability in the HO state is similar to the earlier studies in Ref. [34,54].Ngo et al. [33] considered a two dimensional off-lattice model for the active nematic without the exclusion constraint.In the low and moderate density regime, they show a homogeneous disordered phase and an inhomogeneous chaotic phase, which are similar to the isotropic and the IM states, respectively.Similar to their study, the spanning area of the IM state (golden regime in the phase diagram Fig. 1(a)) along the density axis decreases with the increasing temperature.In the high density limit, they note a homogeneous quasi-ordered phase only, which similar to the HO state in our study.However, we show the bistability between the HO and the IM state in this density limit.
In conclusion, our lattice model for the active nematic is a simple one to design and execute numerically, and easy to compare with the corresponding equilibrium model.It shows new features like the BS in the low temperature regime and the bistability in the high density regime, as well as some of the early characterised states, e.g., the IM state.It also shows many basic features of the active nematic like large number fluctuation, long-time decay of orientation correlation, transition from SRO isotropic to QLRO nematic state.The shift in the transition density due to activity strength compared to the equilibrium model can be tested in experimental systems where activity can be tuned.We expect the emergence of the bistability in the high density regime in a two dimensional experimental system composed of apolar particles with finite dimension and high activity strength.It would be interesting to study the model without volume exclusion.In this study, particle orientation has continuous symmetry of O(2).Therefore, the equilibrium limit of our model is an apolar analogue of the two-dimensional XY-model.One can also study the model with discrete orientation symmetry as in Ref. [20,35,36] and compare the results with the corresponding equilibrium model.
Here we assume the system is aligned along one direction, and the variation in orientation is only along the perpendicular direction.Therefore, we can choose either of equations (B3) or (B4).Two constants c 1 and c 2 are the fluctuations in density when the nematic order parameter is zero.Now from the equation for Π ij , we obtain an effective equation for S as We neglect all the derivative terms and retain only the polynomials in S, i.e., we neglect higher order fluctuations.The Taylor expansion of α 1 (℘) about the mean density ℘ 0 gives α 1 (℘) = α 1 (℘ 0 + δ℘) = α 1 (℘ 0 ) + α 1 δ℘ where α 1 = ∂α1 ∂℘ | ℘0 .This gives We can write an effective free energy F ef f (S) so that Substituting the expression for δ℘ from equation (B4), we obtain Therefore, where and b 4 = 1 2 ℘ 2 0 α 2 .Since the free energy is a state function, we have assumed the integration constant to be zero.Therefore, the fluctuation in the density introduces a cubic order term in the effective free energy F ef f (S).Effective free energy in equation (B9) is similar to the Landau free energy with a new cubic order term [44].Now we calculate the jump ∆S and the new critical density from the coexistence condition for free energy.Steady state solutions of order parameter (S = 0 for isotropic and S = 0 for ordered state) are given by is shifted to a lower density in comparison to the equilibrium transition density ℘ IN .Equation (B14) gives the expression for new transition density as given in the main text.Therefore, using renormalised mean field theory we find a jump ∆S at a lower density as compared to the equilibrium I-N transition density.

Figure 1 .
Figure 1.Phase diagram.(a) Phase diagram for both the equilibrium and the active nematic in the density -inverse temperature plane.The equilibrium system remains in the isotropic (EI) state in the low density regime (on the left of the solid line) and in the nematic (EN) state in the high density regime (on the right of the solid line).The active nematic goes from the disordered isotropic (I) state to the locally ordered inhomogeneous mixed (IM) state with increasing density or decreasing temperature.The I -IM transition occurs with the appearance of clear bands (BS) in the low temperature regime.In the high density regime the active nematic shows bistability between the IM and the homogeneous globally ordered (HO) state.(b) Upper panel shows particle inclination towards the horizontal direction.Colour bar ranging from zero to one indicates vertical to horizontal orientation, respectively.BS is the banded state configuration shown for (β , C) = (2.0,0.38).IM, HO and EN state configurations are shown for (β , C) = (2.0,0.78).Lower panel shows the coarse-grained density in the respective states.

Figure 2 .
Figure 2. Disorder -order transition.(a) Scalar order parameter versus packing density plot for 400 × 400 system size at β = 2.0.The equilibrium system (E, solid line) shows continuous isotropic to nematic state transition with increasing density.The active system goes from the isotropic (I) state to the locally ordered inhomogeneous mixed (IM, ×) state.In the high density regime the system shows bistability between the IM state and the homogeneous globally ordered (HO, •) state.(b) The I to IM transition at low temperature occurs with a jump in S where the particles form bands (BS).Distribution of the scalar order parameter near the I -BS transition at (β , C) = (2.0,0.37) shows two peaks.(c) Finite size scaling of S for both the HO and the IM state at (β , C) = (2.0,0.76).(d) Order parameter time series show that the active system flips in between the HO and the IM state in the bistable regime.Two time series are shown for two different parameter values in the high density regime.

Figure 3 .
Figure 3. Model figure.(a) Two dimensional square lattice with occupied (n = 1) or unoccupied (n = 0) sites.Filled circles indicate the occupied sites.Inclinations of the rods towards the horizontal direction show respective particle orientations θ ∈ [0, π].(b) Equilibrium move: particle can move to any of the four neighbouring sites with equal probability 1/4.(c, d) Active move: particle can move to either of its two neighbouring sites with probability 1/2, if unoccupied, in the direction it is more inclined to, i.e., along BD in (c), and AC in (d).

Figure 2 (
Figure2(a) shows that the I to BS transition occurs in the low temperature regime with a jump in S at a density lower than the corresponding equilibrium I-N transition density C IN .These bands appear because of the large activity strength.A linear stability analysis, as detailed later in this paper, shows that the large activity strength induces an instability in the disordered isotropic state.This instability goes away for small activity strength or at high temperature.We also do a renormalised mean field calculation of an effective free energy written for the active nematic.The calculation predicts a jump in the scalar order parameter and shows a shift in the disordered (S = 0) to ordered (S = 0) state transition density.Both the jump in S and the shift in the transition density reduce with the activity strength or increasing temperature.The I to BS transition is a first order transition.The shift in the disorder-order transition point is a common feature of the active systems.For large activity and low temperature, if the system density is above a certain value but less than C IN , the large density fluctuation present in these systems causes local alignment with local density higher than C IN .Large density fluctuation is an intrinsic feature of the active systems, and as shown in Fig.4, we also observe the same in the ordered active states in our model.Due to activity these locally ordered regions move anisotropically and combine with nearby region with similar local ordering.So larger ordered region forms at mean density lower than the equilibrium I-N transition density.Therefore, we find a disordered to ordered state transition at a lower density.For large activity strength, I-BS transition occurs with the jump in scalar order parameter.In our numerical study we calculate the probability P (S) of the scalar order parameter averaging over many iterations and realisations near the I-BS transition point.Figure2(b)shows P (S) has two peaks, which further supports the first order I-BS transition for large activity strength.In the high density regime (red coloured regime in the phase diagram Fig.1(a)), the AM shows bistability, i.e., it can be either in the locally ordered IM state or in a homogeneous globally ordered (HO) state.As shown in Fig.2(a), the S − C curve for fixed temperature bifurcates in the high density regime; the lower branch corresponds to the earlier discussed IM state, whereas the higher branch indicates the existence of the globally ordered state.Figure1(b) shows that the system possesses less density inhomogeneity in the HO state compared to the IM state.A finite size scaling of both the HO and the IM state, as shown in the Fig.2(c), shows that the active nematic possesses non-zero finite order in both these states.Order parameter time series shown in Fig.2(d) confirms the bistability of the system in the high density regime.Bistability is not generally seen in other agent based numerical simulations of point particles[33]; it appears because of finite filling constraint of the model.This feature can be suppressed if we allow more than one particle to sit together.In the complete filling limit C = 1.0, the AM is equivalent to the EM, and it shows the globally ordered HO state only.

Figure 5 .
Figure 5. Two-point orientation correlation shown for β = 2.0 on log-log scale.(a) Active model: g2(r) decays exponentially at low density ( , ) and algebraically at high density ( , ).In the bistable regime at high density ( ), g2(r) decays algebraically in the HO state and abruptly in the IM state.(b) Equilibrium model: g2(r) decays exponentially at low density ( , ) and algebraically at high density ( , +, ).Continuous lines are the respective fits, fitted for more than one decade.

Figure 6 .
Figure 6.Steady state characteristics of high density states.(a) Orientation distribution P (θ) of particles calculated from one snapshot in the steady state.P (θ) fits with Gaussian distribution (continuous lines) for both the HO and the EN states.The IM state shows broad distribution of θ.(b) Distribution of the mean orientation P ( θ) calculated from θ of each snapshot in the steady state.P ( θ) is broad for the EN state in comparison to the HO state.(c) Steady state autocorrelation Cθ(t) of the mean orientation of the system.All plots are shown for (β , C) = (2.0,0.80).

α 1 = 1 2D℘ , and b 4 = 1 2 ℘ 2 0 α 2 . 2 3
∂α 1 /∂℘| ℘0 , b 3 = a0℘0α Both b 3 and b 4 are positive.A detailed calculation for F ef f is shown in Appendix B. The density fluctuations introduce a new cubic order term in the free energy F ef f (S) that is proportional to the activity strength a 0 .The presence of such term produces a jump ∆S = S c = 2b3 3b4 at a density ℘ c = ℘ IN (1 − 2b 9b4 ) < ℘ IN .Fluctuation in density produces a jump in order parameter and shifts the critical density.Such type of fluctuation induced transitions are called fluctuation dominated first order phase transitions in statistical mechanics

δF ef f δS = −b 2 − b 3 S
+ b 4 S 2 S = 0 (B10) Non-zero S is given by −b 2 − b 3 S c + b 4 S 2 c = 0. Coexistence condition implies F ef f (S c ) = − b jump at the new critical point is ∆S = 2b3 3b4 .Since b 4 > 0 and hence b c 2 < 0, the new critical density