Isotropic-nematic transition of self-propelled rods in three dimensions

Isotropic-nematic transition of self-propelled rods in three dimensions


Using overdamped Brownian dynamics simulations we investigate the isotropic–nematic (IN) transition of self-propelled rods in three spatial dimensions. For two well-known model systems (Gay-Berne potential and hard spherocylinders) we find that turning on activity moves to higher densities the phase boundary separating an isotropic phase from a (nonpolar) nematic phase. This active IN phase boundary is distinct from the boundary between isotropic and polar-cluster states previously reported in two-dimensional simulation studies and, unlike the latter, is not sensitive to the system size. We thus identify a generic feature of anisotropic active particles in three dimensions.

Collective nonequilibrium behavior in suspensions of active Brownian particles (ABPs) is the subject of much current research interest Romanczuk et al. (2012). Not only do these systems exhibit novel dynamics and phase behavior, they are also relevant for understanding self-organization phenomena in nature. Much of the interest in ABPs has been driven by the introduction of new experimental model systems, such as catalytic Janus particles Palacci et al. (2010); Erbe et al. (2008); Howse et al. (2007), light activated colloids Palacci et al. (2013) and colloids with artificial flagella Dreyfus et al. (2005). Additionally, studies of minimal spherical active models have triggered a whole new branch of fundamental research in nonequilibrium statistical mechanics. The striking similarities to an equilibrium system have been exploited by developing a Cahn-Hilliard-like mechanism Speck et al. (2014) to describe the early-stage dynamics of motility-induced phase separation Tailleur and Cates (2008); Buttinoni et al. (2013); Cates and Tailleur (2015), identifying an effective equilibrium regime Maggi et al. (2015); Farage et al. (2015); Fodor et al. (2016); Wittmann et al. (2017), defining effective interaction potentials Cates and Tailleur (2015); Schwarz-Linek et al. (2012); Wittmann and Brader (2016); Wittmann et al. (2017) or employing linear-response theory Sharma and Brader (2016). More fundamentally, a better understanding of active pressure Takatori et al. (2014); Solon et al. (2015) or chemical potential Dijkstra et al. (2016) is required to provide a solid framework for active thermodynamics Speck (2016); Solon et al. (2016).

While spherical ABPs are ideal for exploring basic concepts, suspensions of anisotropic active Brownian particles are perhaps more relevant, as these better represent the generic type of particles encountered in nature Elgeti et al. (2015); Marchetti et al. (2013). Self-propelled rods (SPRs), the anisotropic analog of ABPs, for which the self propulsion along the long axis of the particle breaks the up-down symmetry, exhibit a rich dynamical phase behavior at high activity Wensink and Löwen (2012). Simulations of large two–dimensional systems (with rotational diffusion) Yang et al. (2010); Abkenar et al. (2013); Weitz et al. (2015) reveal that at densities below the passive isotropic–nematic (IN) transition, the initially isotropic state begins to destabilize due to the emergence of moving polar clusters, which grow in size upon increasing activity but do not form a global phase Weitz et al. (2015). Experiments on a fluidized monolayer of rods have identified giant number fluctuations in such states Narayan et al. (2007). The (enhanced) nematic ordering of a biologically inspired two-dimensional nematic model has been studied in simulation Kraikivski et al. (2006); however, the IN phase boundary for overdamped SPRs has not been addressed explicitly. To avoid confusion of terminology we note that in the literature ‘active nematics’ usually address anisotropic particles driven randomly back and forth along their axis Marchetti et al. (2013). In this sense the term ‘nematic’ refers to the particle symmetry. The present work concerns SPRs and the terms ‘polar’ and ‘nematic’ will be reserved to describe collective states.

The theoretical understanding of the phase behavior of SPRs is a difficult problem. According to an early mean-field approach Baskaran and Marchetti (2008a), the density at the IN phase transition is insensitive to activity. In contrast, more general collision-based models Baskaran and Marchetti (2008b); Baskaran and Cristina Marchetti (2010) predict that the transition density decreases with increasing activity. For overdamped (Langevin) dynamics, the current numerical evidence suggesting that activity might stabilize nematic order of SPRs only arises from the observation that, as the density increases, the destabilization of the isotropic phase with respect to polar fluctuations occurs at lower activities Yang et al. (2010); Weitz et al. (2015). In general, the existence of a (nonpolar) nematic phase and its phase boundary remain an open problem.

Figure 1: Isotropic–nematic (IN) transition of self-propelled rods (SPRs). The solid line and the filled region corresponds to soft-repulsive Gay-Berne (GB) ellipsoids and the dashed line and symbols to hard spherocylinders (HSC). The particles, sketched on the right indicating their (approximate) length-to-width ratio, are propelled with constant velocity in the direction of their orientation (red arrows). In both systems, the state points are given by the modified Péclet number (swim speed rescaled by a particle-geometry dependent factor) and the packing fraction , relative to the IN transition packing fraction in equilibrium. For large we also find smectic clusters for HSC.

In this letter we address the activity dependence of the IN transition close to equilibrium in three dimensions using overdamped Brownian dynamics simulations. Although we will focus here on the repulsive Gay-Berne (GB) model, which is a classic model of thermotropic liquid crystals Gay and Berne (1981); Adams et al. (1987); Luckhurst et al. (1990), we have also performed simulations of lyotropic hard spherocylinders (HSC) Bolhuis and Frenkel (1997); Wittmann et al. (2016) to ensure that our findings are not model-specific. We will concentrate on systems of lower aspect ratio, well away from the Onsager limit. For both of the considered models we observe in Fig. 1 a shift of the IN transition towards higher densities as the system is driven out of equilibrium by turning on the activity. The IN transition line is independent of system size, in contrast to the transition between isotropic and polar-cluster states found at higher activities.

The reason that the active IN transition has not been previously reported is that all comparable simulation studies Yang et al. (2010); Abkenar et al. (2013); Weitz et al. (2015) of active liquid crystals have focused on (pseudo) two-dimensional systems at relatively high activity, for which only isotropic and (system-size dependent) polar-cluster states are observed. It is known that certain collective phenomena, such as motility-induced phase separation Stenhammar et al. (2014), occur at lower activity in 2D than in 3D, because particles can pass each other more easily in higher dimensions. It could thus be anticipated that nematic order in lower dimensional systems will be correspondingly more fragile with respect to activity. This expectation is validated by our finite-size simulations of the GB model in 2D (see the supplementary information (SI) sup ()), which show that nematic order is rapidly disrupted by turning on activity. Moreover, the scaling arguments of Ramaswamy et al. Ramaswamy et al. (2003) predict that global nematic order cannot withstand a finite activity: the relative root-mean-square fluctuations in the particle number in a subvolume scale as , independent of the propulsion speed Ramaswamy et al. (2003); Marchetti et al. (2013). For this does not vanish in the thermodynamic limit and there is thus no homogeneous nematic phase. To study a true active IN transition three-dimensional systems must be considered.

In our simulations, we consider a system of interacting, anisotropic particles for which the position and orientation of particle are specified by the vector and unit vector , respectively (see the sketch in Fig. 1). These vectors evolve in time according to the coupled Langevin equations in the overdamped limit,


where and are friction coefficients determined by the translational and rotational diffusion coefficients and . The force and torque acting on particle are related to the total potential energy, which we take to be a sum of pair potentials, , according to and . The stochastic vectors and are Gaussian distributed with zero mean and have the time correlations and . The active component of the dynamics enters via the second term in (1), where is a constant self-propulsion velocity. We define dimensionless time, , and activity, , where is a typical length-scale of the studied interaction potentials (see SI sup ()).

The GB potential Gay and Berne (1981); Adams et al. (1987); Luckhurst et al. (1990) between a pair of particles is of the form of a Lennard-Jones potential, whose depth and range depend on the separation and mutual orientation. Here we consider a purely soft-repulsive WCA-like version Rull (1995), obtained by shifting and truncating the GB potential, where the length-to-width ratio is roughly given by . We also consider spherocylinders interacting via a nearly hard-core potential with aspect ratio (well-studied in equilibrium). The massively parallel framework we employ allows us to efficiently gather statistics by considering very large systems. A detailed description of both models and sets of simulations is given in the SI sup (), see also the sketches in Fig. 1.

Figure 2: Phase diagram of the three-dimensional GB model in the packing fraction–activity plane. The stability of the isotropic (I), nematic (N) and local polar (P) states is determined by analyzing global order parameters. Typical snapshots are shown for three distinct state points, as indicated by triangles. The IN boundary (see Fig. 1 for a closeup) is insensitive to changes in the system size and separates true nonequilibrium phases. In contrast, the boundary separating I and P states is system-size dependent; the region P merely indicates where system-spanning polar clusters occur.

We analyze the orientational behavior of the system by measuring the time averages and of the nematic order parameter Eppenga and Frenkel (1984) and the polar order parameter , respectively, defined at each instant of time


where is the orientation vector of particle and is the nematic director. Both quantities take values between and , with the extreme values indicating full disorder or perfect order, respectively. To distinguish between the different states we choose the threshold values for the onset of nematic and for polar order.

In the GB model, we simulate particles and define the packing fraction , with being the volume of the simulation box and the volume of one ellipsoidal particle (see also SI sup ()). Our simulations yield a passive IN transition at , in agreement with previous studies Rull (1995). The resulting phase diagram, shown in Fig. 2, reveals three distinct states in the density-activity plane: isotropic, nematic and polar. We characterize the polar state, in which the majority of particles are driven in the same direction, by and . Its occurrence here is a known artifact of a finite system Yang et al. (2010); Abkenar et al. (2013); Weitz et al. (2015), which we detail below.

Our main result is that we observe a nonequilibrium nematic phase with but , whose boundary bends to the right with increasing , suggesting that introducing a moderate amount of activity can suppress orientational ordering. Simulations at selected state points with larger numbers of particles (N=1000) verify that the location of this active IN transition is independent of the system size (see SI sup ()). To check that the active IN transition is not specific to the GB model we compare the results to large-scale simulations of the active HSC model, where , using particles. For both systems the phenomenology depicted in Fig. 1 is consistent; the IN transition line moves to higher densities as the activity is increased.

Figure 3: Simulation snapshots of the isotropic (I) and nematic (N) phases for the active GB system with and (a) and (b) and for the active HSC system at packing fractions just (c) below and (d) above the transition at (see Fig. 1)


In order to make a proper comparison we define a dimensionless swimming speed that does not contain arbitrary length- and time-scales. We thus consider the (square root of) the active part of the single-particle diffusivity relative to the passive part,  ten Hagen et al. (2011). Furthermore, we rescale the density by its value at the equilibrium IN transition for each system. Despite the similarity of the rescaled GB and HSC phase boundaries in Fig. 1 there are quantitative differences arising presumably from both the interparticle interactions and the aspect ratio. Most notably, we observe in Fig. 3 a different microstructure. The nematic phase in the HSC system does not differ much from the equilibrium nematic, even quite close to the phase boundary, which results in a sharp IN transition also for a finite activity (see (SI) for more details sup ()). In contrast, the GB system begins to exhibit local polar order as the IN phase boundary is approached, which we now analyze in more detail.

Figure 4: Nematic and polar order parameter used to determine the phase behavior in the active GB simulation. The state points indicated by a square and a circle correspond to the snapshots in Figs. 3a and b, respectively. (a) Time average of (the inset shows ) for different velocities as a function of packing fraction . The broken line indicates our threshold for determining the IN phase boundary. (b) Time evolution of both S and P for and .

In Fig. 4a we show the time-averaged global order parameters from Eq. (3) at packing fractions close to the IN phase boundary in the active GB system. The packing fraction at which decreases below the chosen threshold, , shifts to higher values as the activity is increased. This is accompanied by a reduced slope of and significant orientational fluctuations in the active system. It is thus increasingly difficult to unambiguously resolve the phase boundary with our method. As illustrated by the behavior of in the inset of Fig. 4a, the emergence of nematic order is not associated with global polar ordering. Following a path of state points by increasing at a fixed packing fraction leads to a decrease in the nematic order parameter , which eventually falls below our chosen threshold and we classify the state as isotropic, i.e., the activity destabilizes the nematic phase.

To better understand the behavior of the active system, it is important to discuss the role of fluctuations. In the global isotropic phase it is well known Yang et al. (2010); Abkenar et al. (2013); Weitz et al. (2015) that there emerge local polar clusters with a critical size, which increases upon increasing the activity or the density. The local polar state depicted in Fig. 2 for the GB ellipsoids thus corresponds to a single cluster spanning the whole system. On increasing the system size, the associated “phase boundary” shifts to higher for a given (see SI), which consistently verifies that the polar state in our finite-size simulation does not represent a true nonequilibrium phase with global order in an infinite system Yang et al. (2010); Abkenar et al. (2013); Weitz et al. (2015). In the states which we characterize as active nematic the polar fluctuations are much more prominent than one would expect if the state points were isotropic. In Fig. 4b we show the temporal fluctuations of both order parameters associated with the nematic snapshot in Fig. 3b. These rationalize the smaller slope of the time averaged value of at higher activity, as observed in Fig. 4a.

It might be suspected that the fluctuations discussed above for the GB ellipsoids are related to an enhancement of unphysical, finite-size induced self-interactions due to the persistent motion of the aligned rods. However, the following considerations support our claim that the IN phase boundary depicted in Fig. 1 is generic (see (SI) sup ()). Firstly, we stress that the observed long-time behavior is independent of the (either polar or isotropic) initial conditions. Secondly, upon further increasing the activity, the nematic phase eventually turns into a well-defined isotropic phase with significantly fewer fluctuations, which points to a well-defined phase transition even if the fluctuations in the nematic phase partially arise from finite-size effects. Finally, both the lifetime and the magnitude of the described fluctuations decrease with increasing system size, whereas the nematic order parameter is robust. This point is corroborated by our large-scale simulations of the HSC system, where we do not observe significant fluctuations of the global order parameters and typical nematic configurations do not show strong local polar ordering, even at swimming speeds where the isotropic phase does clearly exhibit large polar clusters, compare Figs. 3c and d.

In conclusion, we identified in three dimensions and for small aspect ratios a homogeneous nematic phase, which can be clearly distinguished from the isotropic phase, even in a relatively small system. This nonequilibrium nematic order is gradually destabilized by activity. Our finding is not sensitive to the precise particle shape or the details of the interaction, provided the rods are short. The activity-induced stabilization of the nematic phase predicted by mean-field theory Baskaran and Marchetti (2008b) is thus not universal.

Macroscopically, the nematic phase is entirely different in two dimensions, where giant number fluctuations occur and no homogeneous phase exists Ramaswamy et al. (2003). Nevertheless, the microscopic mechanisms, leading to the stabilization of nematic order inside the high-density regions should also be at work in the homogeneous phase, described in the present work. In systems of long rods, especially in two dimensions, head-to-side collisions dominate and will rotate the particles towards either a parallel or anti-parallel orientation, thus enhancing the nematic order. As the aspect ratio is reduced, head-on collisions become increasingly frequent, generating disorder and destabilizing the nematic phase. Moreover, as the aspect ratio is reduced the passive IN transition moves to higher densities, which further increases the relative importance of head-on collisions. The consequences of dense clustering and correlations beyond the mean-field level have not been taken into account in previous theoretical studies Baskaran and Marchetti (2008b).

Hydrodynamic theory predicts that the nematic phase is unstable in an infinite system of active particles with respect to inhomogeneous flows with large wave length Aditi Simha and Ramaswamy (2002); Marchetti et al. (2013). However, the simulation technique employed by us, Brownian dynamics, implies a fixed background with which momentum is exchanged through the noise and friction terms in the equations of motion. In this case large wavelength instabilities will be suppressed Aditi Simha and Ramaswamy (2002); Marchetti et al. (2013). In experiments, the stick-boundary conditions at the walls might also act to similarly suppress these instabilities at the low activities we consider, if the system is not too large. This suppression, which could be a topic of future research, should lead to a system that is qualitatively similar to the one we consider here.

Regarding the bulk system, it will be of interest to explore in detail the influence of both interparticle interactions and aspect ratio on the active IN phase boundary. On the theoretical side, it would be desirable to have a first-principles theoretical approach to describe the activity dependence of the IN phase boundary, even if this is limited to low activities, close to equilibrium. While the phase behavior of spherical ABPs can be explained solely by effective attractions Cates and Tailleur (2015); Schwarz-Linek et al. (2012); Wittmann and Brader (2016) and that of active nematic rods by an effective (longer) aspect ratio Kraikivski et al. (2006), an appropriate effective potential for SPRs should account for their characteristic broken up-down symmetry. The most simplistic passive model system with this property consists of hard pear-shaped objects, for which it has been detailed recently that the nematic phase destabilizes with increasing deviation from ellipsoidal shape Schönhöfer et al. (2017). This observation suggests an intuitive mapping to describe the IN transition in qualitative agreement with our simulations, which is yet to be quantified. Finally, there is much opportunity for further experimental or numerical studies of the bulk elasticity, the response to (time-dependent) external fields and inhomogeneous systems of SPRs in the active nematic phase.


  1. P. Romanczuk, M. Bär, W. Ebeling, B. Lindner,  and L. Schimansky-Geier, The European Physical Journal Special Topics 202, 1 (2012).
  2. J. Palacci, C. Cottin-Bizonne, C. Ybert,  and L. Bocquet, Physical Review Letters 105 (2010), 10.1103/physrevlett.105.088304.
  3. A. Erbe, M. Zientara, L. Baraban, C. Kreidler,  and P. Leiderer, Journal of Physics: Condensed Matter 20, 404215 (2008).
  4. J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh,  and R. Golestanian, Physical Review Letters 99 (2007), 10.1103/physrevlett.99.048102.
  5. J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine,  and P. M. Chaikin, Science 339, 936–940 (2013).
  6. R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone,  and J. Bibette, Nature 437, 862–865 (2005).
  7. T. Speck, J. Bialké, A. M. Menzel,  and H. Löwen, Physical Review Letters 112 (2014), 10.1103/physrevlett.112.218304.
  8. J. Tailleur and M. E. Cates, Physical Review Letters 100 (2008), 10.1103/physrevlett.100.218103.
  9. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger,  and T. Speck, Physical Review Letters 110 (2013), 10.1103/physrevlett.110.238301.
  10. M. E. Cates and J. Tailleur, Annual Review of Condensed Matter Physics 6, 219–244 (2015).
  11. C. Maggi, U. M. B. Marconi, N. Gnan,  and R. Di Leonardo, Scientific Reports 5 (2015), 10.1038/srep10742.
  12. T. F. F. Farage, P. Krinninger,  and J. M. Brader, Physical Review E 91 (2015), 10.1103/physreve.91.042310.
  13. É. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco,  and F. van Wijland, Physical Review Letters 117 (2016), 10.1103/physrevlett.117.038103.
  14. R. Wittmann, C. Maggi, A. Sharma, A. Scacchi, J. M. Brader,  and U. M. B. Marconi, Journal of Statistical Mechanics: Theory and Experiment 2017, 113207 (2017).
  15. J. Schwarz-Linek, C. Valeriani, A. Cacciuto, M. E. Cates, D. Marenduzzo, A. N. Morozov,  and W. C. K. Poon, Proceedings of the National Academy of Sciences 109, 4052–4057 (2012).
  16. R. Wittmann and J. M. Brader, EPL (Europhysics Letters) 114, 68004 (2016).
  17. A. Sharma and J. M. Brader, The Journal of Chemical Physics 145, 161101 (2016).
  18. S. Takatori, W. Yan,  and J. Brady, Physical Review Letters 113 (2014), 10.1103/physrevlett.113.028103.
  19. A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. E. Cates,  and J. Tailleur, Physical Review Letters 114 (2015), 10.1103/physrevlett.114.198301.
  20. M. Dijkstra, S. Paliwal, J. Rodenburg,  and R. van Roij,  (2016), 1609.02773 .
  21. T. Speck, EPL (Europhysics Letters) 114, 30006 (2016).
  22. A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri,  and J. Tailleur,  (2016), 1609.03483 .
  23. J. Elgeti, R. G. Winkler,  and G. Gompper, Reports on Progress in Physics 78, 056601 (2015).
  24. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao,  and R. A. Simha, Reviews of Modern Physics 85, 1143–1189 (2013).
  25. H. H. Wensink and H. Löwen, Journal of Physics: Condensed Matter 24, 464130 (2012).
  26. Y. Yang, V. Marceau,  and G. Gompper, Physical Review E 82 (2010), 10.1103/physreve.82.031904.
  27. M. Abkenar, K. Marx, T. Auth,  and G. Gompper, Physical Review E 88 (2013), 10.1103/physreve.88.062314.
  28. S. Weitz, A. Deutsch,  and F. Peruani, Physical Review E 92 (2015), 10.1103/physreve.92.012322.
  29. V. Narayan, S. Ramaswamy,  and N. Menon, Science 317, 105–108 (2007).
  30. P. Kraikivski, R. Lipowsky,  and J. Kierfeld, Physical Review Letters 96 (2006), 10.1103/physrevlett.96.258103.
  31. A. Baskaran and M. C. Marchetti, Physical Review E 77, 011920 (2008a).
  32. A. Baskaran and M. C. Marchetti, Physical Review Letters 101 (2008b), 10.1103/physrevlett.101.268101.
  33. A. Baskaran and M. Cristina Marchetti, Journal of Statistical Mechanics: Theory and Experiment 2010, P04019 (2010).
  34. J. G. Gay and B. J. Berne, The Journal of Chemical Physics 74, 3316 (1981).
  35. D. Adams, G. Luckhurst,  and R. Phippen, Molecular Physics 61, 1575 (1987).
  36. G. R. Luckhurst, R. A. Stephens,  and R. W. Phippen, Liquid Crystals 8, 451 (1990).
  37. P. Bolhuis and D. Frenkel, The Journal of chemical physics 106, 666 (1997).
  38. R. Wittmann, M. Marechal,  and K. Mecke, Journal of Physics: Condensed Matter 28, 244003 (2016).
  39. J. Stenhammar, D. Marenduzzo, R. J. Allen,  and M. E. Cates, Soft Matter 10, 1489–1499 (2014).
  40. “See supplementary material at (link will be supplied by the editor) supplemental for technical details about the simulations and preliminary results in two dimensions.” .
  41. S. Ramaswamy, R. A. Simha,  and J. Toner, EPL (Europhysics Letters) 62, 196 (2003).
  42. L. F. Rull, Physica A: Statistical Mechanics and its Applications 220, 113 (1995).
  43. R. Eppenga and D. Frenkel, Molecular Physics 52, 1303 (1984).
  44. B. ten Hagen, S. van Teeffelen,  and H. Löwen, Journal of Physics: Condensed Matter 23, 194119 (2011).
  45. R. Aditi Simha and S. Ramaswamy, Phys. Rev. Lett. 89, 058101 (2002).
  46. P. W. A. Schönhöfer, L. J. Ellison, M. Marechal, D. J. Cleaver,  and G. E. Schröder-Turk, Interface Focus 7 (2017), 10.1098/rsfs.2016.0161.
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Comments 0
The feedback must be of minumum 40 characters
Add comment

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question