Isotropic-nematic transition of self-propelled rods in three dimensions
Matthias C. Bott Felix Winterhalter Matthieu Marechal Abhinav Sharma Joseph M. Brader René Wittmann
cond-mat.soft cond-mat.stat-mech






Soft Matter Theory, University of Fribourg, CH-1700 Fribourg, Switzerland
matthiaschristian.bott@unifr.ch


Institut für Theoretische Physik, Universität Erlangen-Nürnberg, 91058 Erlangen, Germany


Institut für Theoretische Physik, Universität Erlangen-Nürnberg, 91058 Erlangen, Germany


Leibniz-Institut für Polymerforschung Dresden, 01069 Dresden, Germany


Soft Matter Theory, University of Fribourg, CH-1700 Fribourg, Switzerland


Soft Matter Theory, University of Fribourg, CH-1700 Fribourg, Switzerland




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. 


Isotropic-nematic transition of self-propelled rods in three dimensions
    R. Wittmann
    December 30, 2023
=======================================================================





§ INTRODUCTION


Collective nonequilibrium behavior in suspensions of active Brownian particles (ABPs) is the subject of much 
current research interest <cit.>. 
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 <cit.>, light 
activated colloids <cit.> and colloids with artificial flagella <cit.>. 
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 <cit.> to describe the early-stage dynamics of motility-induced phase separation <cit.>,
identifying an effective equilibrium regime <cit.>,
defining effective interaction potentials <cit.>
or employing linear-response theory <cit.>.

More fundamentally, a better understanding of active 
pressure <cit.> or chemical potential <cit.> is required
to provide a solid framework for active thermodynamics <cit.>.
 Recently, also the question of how activity influences the well-studied phase transitions in a passive system of soft disks has been addressed in detail <cit.>.


While spherical ABPs are ideal for exploring basic concepts, 
suspensions of anisotropic ABPs are perhaps more relevant, 
as these better represent the generic type of particles encountered in nature 
<cit.>. 
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 (infinite)
activity <cit.> in two dimensions (2D). 
Simulations of large 2D systems (with rotational diffusion) <cit.> 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 <cit.>. At higher densities, a laning phase is found,
which does have nematic order on the range of the simulation domain, but is not homogeneous <cit.>.
Experiments on a fluidized monolayer of rods 
have identified giant number fluctuations in such states <cit.>. 
 For experiments on very long (and thus non-Brownian) bacteria in quasi-2D, a nematic phase with long-range order was reported <cit.>.
 
The (enhanced) nematic ordering of a biologically inspired 2D nematic model has been studied in 
simulation <cit.>. In extensions of the Vicsek model to incorporate local nematic ordering (rather than the polar ordering of the original Vicsek model), the region of stability of homogeneous active nematic phase (with giant
number fluctuations) in 2D is determined <cit.>.
It is unsure whether this phase has 
long-range or quasi-long-range order in these agent-based simulations; it even seems to depend on the details of the model (compare Refs. Ginelli2010 and Ngo2014).
 Returning to the SPRs interaction model, the lack of observations of a homogeneous nematic phase in previous simulations (in 2D) means that the IN phase boundary 
for overdamped SPRs has not been addressed explicitly. 
 To avoid confusion of terminology we emphasize that in the literature
`active nematics' usually address anisotropic particles driven randomly 
back and forth along their axis <cit.>. 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 <cit.>, the density at the IN phase transition is insensitive to activity. 
In contrast, more general collision-based models <cit.> 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 <cit.>. 
In general, the existence of a (nonpolar) nematic phase and its phase boundary 
remains an open problem.



In this paper we address the activity dependence of the IN transition of SPRs close to equilibrium using 
overdamped Brownian dynamics simulations. We have chosen to mostly work in three dimensions (3D), because (i) we expected a greater stability of the homogeneous nematic phase
in 3D than in 2D (as for the phase-separated state of ABPs <cit.>) (ii) giant number fluctuations are predicted to be reduced in 3D compared to 2D for active nematics <cit.>
(it should be noted that the validity of the linearized theory from which these predictions stem is debated for 2D systems <cit.>)
and (iii) the equilibrium nematic phase has long-ranged order in 3D (as opposed to quasi-long-range order in 2D), 
which makes it easier to describe using particle-resolved computer simulations and other theoretical approaches—leading to less finite-size effects 
in the former—than in 2D.
We will further consider systems of lower aspect ratio, well away from the Onsager limit. 

Although we will focus on the repulsive Gay-Berne (rGB) model,
which is a classic model of thermotropic liquid crystals <cit.>, 
we have also performed simulations of lyotropic hard 
spherocylinders (HSC) <cit.>
to ensure that our findings are not model-specific.
For both of the considered models we observe a shift of the IN transition towards higher 
densities as the system is driven out of equilibrium by turning on the activity, see Fig. <ref>. 
The IN transition line is independent of system size, in contrast to the 
transition between isotropic and polar-cluster states found at higher activities. 

A preliminary survey of simulations of the rGB model in 2D (see the supplementary information (SI) <cit.>)
indicates that nematic order is more rapidly disrupted by turning on activity than in 3D, 
but finite size effects preclude a definite statement.
 This might be the reason why such an active IN transition has not been reported earlier.

The remainder of this paper is arranged as follows.
 In Sec. <ref> we present the two different models we considered
and describe the simulation method employed for each of them.
In Sec. <ref> we present the simulation results of the two models and identify the nematic phase at finite activity.
Finally, in Sec. <ref> we discuss our findings and provide an outlook.






§ SIMULATION METHODS





 §.§ Gay-Berne model


The regular Gay-Berne potential <cit.> between a pair of particles 
is of the form of a Lennard-Jones potential, whose depth and range depend on the inter-particle separation and the particle orientations:


    ϕ_ gb(û_1, û_2, 𝐫)    = 
       4ϵ(û_1, û_2, 𝐫̂)
     
     {( σ_0/r - σ(û_1, û_2, 𝐫̂) + σ_0)^12
       - ( σ_0/r - σ(û_1, û_2, 𝐫̂) + σ_0)^6}.


Here the unit vectors û_1 and û_2 specify the orientation of the interacting particles 1 and 2 and 𝐫, r, 𝐫̂ their center to center-vector, -distance and -direction.
The attraction depth ϵ(û_1, û_2, 𝐫̂) and the range σ(û_1, û_2, 𝐫̂) of the particle interaction are dependent on the orientation. 
The legthscale is σ_0, which in our case is the width of the particle (see also appendix <ref> for detailed description).

Here we consider a purely soft-repulsive WCA-like version <cit.>, the rGB model


    ϕ_ rgb   (û_1,û_2,𝐫) 
       =ϕ_ gb(û_1,û_2,𝐫) + ϵ(û_1,û_2,𝐫)    r<r_ min
    
        0                          r ≥ r_ min,

obtained by shifting 
and truncating the Gay-Berne potential.

We study particles with a length-to-width ratio which is roughly given by κ^-1=3 (see the sketch in Fig. <ref>).
The packing fraction is defined as η = N V_E/V, with V being the volume of the simulation box and V_E = 4πσ_0^3/(3κ) the volume of the ellipsoidal particle.
 
The position and orientation vectors evolve in time according to the coupled Langevin equations in the overdamped limit,


    d/dt𝐫_i(t) =    γ^-1 𝐅_i(t)  + v_0 û_i(t) + ξ_i(t), 
    d/dtû_i(t) =    α^-1(𝐓_i(t) ×û_i ) + (η_i(t) ×û_i(t)).


where γ=k_B T/D_T and α=k_B T/D_R are friction coefficients determined by the 
translational and rotational diffusion coefficients D_T and D_R, which we define as D_T/σ_0^2 = 3 D_R.
The force and torque acting on particle i are related to the total potential energy, which we take to be 
a sum of pair potentials, U_N=∑_i<jϕ(û_i, û_j, 𝐫_ij), 
according to 
𝐅_i(t) = -∇_i U_N and 𝐓_i(t) = - û_i ×∂ U_N/∂û_i. 
The stochastic vectors ξ_i(t) and η_i(t) are Gaussian distributed with 
zero mean and have the time correlations 
⟨ξ_i(t)⊗ξ_j(t') ⟩ = 2 D_T1δ_ijδ(t-t')
and
⟨η_i(t) ⊗η_j(t') ⟩ = 2 D_R1δ_ijδ(t-t'), where ⊗ refers to the dyadic product of two vectors. 
The active component of the dynamics enters via the second term in (<ref>), where v_0 is a constant self-propulsion velocity.
We define the dimensionless time t^*= t D_T/σ_0^2, and activity as v_0^* = v_0 σ_0 /D_T.



 §.§ Hard spherocylinders (HSC)


We also consider HSC interacting via a nearly hard-core 
potential with aspect ratio l=5 (well-studied in equilibrium) (see also the sketches in Fig. <ref>).

The aspect ratio l=L/σ_0 of a HSC is given by the ratio of the cylinder length L and the diameter σ_0 of the capping hemispheres.
Its volume is thus given by V_S=π(l/4+1/6) σ_0^3 and the packing fraction follows as η = N V_S/V.
For the system of active HSC, 
we use an existing simulation framework, the  part of waLBerla <cit.>. The  part is a massively parallel framework for molecular dynamics (MD) 
and a similar technique, the discrete element method; we do not use the Lattice Boltzmann technique (for which the waLBerla framework is better known) in this work.
We implemented a friction and noise term in this MD framework (while keeping the inertia term), which means that the Langevin equation we are using is not fully overdamped:

    d𝐩_i/dt   = -Ξ_i ·𝐯_i +𝐅_S,i +𝐅_R,i
    d𝐋_i/dt   = -γ_r ω_i +𝐓_S,i +𝐓_R,i ,

where 𝐩_i=m𝐯_i and 𝐋_i=I_i ·ω_i are the momentum and angular momentum of a particle with mass m and inertia tensor I_i. 
Here 𝐯_i and ω_i denote the translational and angular velocity, respectively.
The first term on the right hand side of Eq. (<ref>) and (<ref>) accounts for friction due to the viscous dissipation. The translational friction tensor Ξ_i depends on the 
translational friction coefficients γ_ and γ_ for motion parallel and perpendicular to the symmetry axis û_i of particle i: 


    Ξ_i=γ_û_i⊗û_i+γ_(𝕀-û_i⊗û_i) ,


(where 𝕀 denotes the identity matrix). For reasons of symmetry, only angular velocities ω_i perpendicular to the symmetry axis û_i of the HSC are considered, therefore,
the rotational friction coefficient γ_r for rotation of the particle axis suffices to describe the viscous torque. 

The subscript S in Eq. (<ref>) and (<ref>) indicates the systematic contributions to the force 𝐅_i and torque 𝐓_i, respectively.
The first contribution to the systematic force 𝐅_S,i are the particle interactions:
The particles interact only when they intersect. Overlaps are resolved by applying a fully elastic linear spring force model to any contact points. The restitution force F_rest in the direction of the contact
normal 𝐧 acting at the contact point is given by
F_rest = k δ
with k the stiffness of the potential and δ the penetration depth. We set the stiffness k to a high value: β kσ_HSC^2= 4· 10^4,
such that more than 99% of the collisions at the higher densities have a penetration δ/σ_HSC<0.02. 
The second contribution to 𝐅_S,i models the self-propulsion: 𝐅_SP,i= γ_ v_0 û_i.

The random contributions 𝐅_R,i and 𝐓_R,i, originating from collisions with solvent molecules as mentioned above, have a Gaussian probability distribution. 
The corresponding correlation functions are related to the viscous friction according to the fluctuation-dissipation theorem for particles i and j:



    ⟨𝐅_R,i(t)⟩   =⟨𝐓_R,i(t)⟩=0⟨𝐅_R,i(t)⊗𝐅_R,j(t')⟩   =2  T Ξ_i δ_ijδ(t-t')
    ⟨𝐓_R,i(t)⊗𝐓_R,j(t')⟩   =2  T γ_r(𝕀-û_i⊗û_i) δ_ijδ(t-t')
    ⟨𝐅_R,i(t)⊗𝐓_R,j(t')⟩   =0


with Boltzmann's constant  and the temperature T. δ(t-t') represents the Dirac delta distribution, which in the case of discrete time steps of size dt is replaced by δ_tt'/dt. 
Since the angular velocity ω_i is kept perpendicular to the symmetry axis of the HSC, only random torques 𝐓_R,i normal to this axis are applied.

We define the mass m and the moment of inertia I_i of the particle such that the relaxation time 
τ_I = m/γ̅ with 1/γ̅ = 1/γ_ + 2/γ_
for the (linear) momentum is 100 times smaller 
than the Brownian time scale σ_0^2/D_T
and effects of the inertia are thus expected to be small. For the HSC system, the translational diffusion constant is given by D_T=k_B
T/γ̅.






 §.§ Parameters



We analyze the orientational behavior of the system by measuring the time averages S and P of the nematic order 
parameter S(t) <cit.> and the polar order parameter P(t), respectively, defined at each instant of time t as


    S(t) = 1/N∑_i=1^N 3 (û_i ·𝐧̂)^2 - 1/2 ,   P(t) = 1/N| ∑_i=1^N û_i ·𝐧̂ | ,


where û_i is the instantaneous orientation vector of particle i and 𝐧̂ 
is the nematic director, see appendix <ref>.

Both quantities take values between 0 and 1, with the extreme values indicating full disorder or perfect order, respectively.  
To distinguish between the different states we choose the 
threshold values S_t=0.35 for the onset of nematic and P_t=0.55 for polar order.

For the rGB ellipsoids, different choices for these thresholds would lead to a slight shift of the 
transition lines in the phase diagram, but not affect our main conclusions. For the HSC, the IN 
transition is more strongly discontinuous leading to a large jump in the order parameter, so the location of the IN
transition in the phase diagram is not affected by small changes of the threshold S_t. In the HSC system,
we never observe global polar order (due to the system size and the relatively low activities), so the value of P_t is irrelevant.


In the rGB model, we simulate N=500 particles for various activities and densities.
The size of the simulation box is determined by the packing fraction which ranges from η = 0.1 to η = 0.59 in our computations.
This corresponds to a side length of the simulation box ranging between b ≈ 19.9 σ_0 to b ≈ 11.0 σ_0.

Our simulations yield a passive IN transition at η = 0.53, which is in agreement with 
previous studies <cit.>.
To rule out finite-size effects, we repeated some simulation runs for N=1000 particles.

To show that the active IN transition is not specific to the rGB model we compare the results to larger-scale simulations of the active HSC model using N≃ 40000 particles.
The IN transition in the passive HSC system lies at η=0.415, 
in good agreement with previous work <cit.>.
For a finite activity we have also performed some simulations for N=20000 particles and observed no significant changes in the results compared to N=40000.


To make a proper comparison of both systems 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 <cit.> and define the dimensionless

    Pe'=√(D_a/D_p) = v_0/√(6 D_T D_R)=v_0^*√(D_T)/σ_0√(6  D_R) .

Furthermore, we rescale the density by its value at the equilibrium IN transition for each system. 







§ NUMERICAL RESULTS




For the rGB model, we mapped out an exemplary full finite-size phase diagram,
shown in Fig. <ref>, which 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 S>S_t and P>P_t.
Its occurrence here is a known artifact of a finite system <cit.>, 
which we detail below. The purpose of showing it here is to indicate the onset of large polar cluster formation (a `large' cluster contains a few hundred particles).

Outside of this region we can expect that the simulation results for the rGB model indicating an isotropic phase are trustworthy.
The large finite-size effects in this polar state are not to be confused with the distinct polar fluctuations observed at relatively low activity near the IN phase boundary, which arise due to a combination of
the finite-size effects and the enhanced tendency of the rods to align parallel within the active nematic phase,
indicate a true phase transition.


Our main result is that we observe a nonequilibrium nematic phase with S>S_t but P<P_t,
whose boundary bends to the right with increasing v_0, suggesting that introducing a moderate amount of activity can suppress 
orientational ordering. 

For both the rGB and HSC systems the phenomenology depicted in Fig. <ref> is consistent; the 
IN transition line moves to higher densities as the activity is increased. 
We explicitly verified that the location of this active IN transition in each model is independent of the system size,
as indicated in Fig. <ref>.



Despite the similarity of the rescaled rGB and HSC phase boundaries in Fig. <ref> there are 
quantitative differences presumably related to both the interparticle interactions and the aspect ratio. 
Most notably, we observe in Fig. <ref> a different microstructure and
the rGB system begins to exhibit local polar order as the IN phase boundary is approached.
 In addition, the equilibrium phase diagram differs in the two models in that a crystal is found for ellipsoids at higher densities, while the 
phase diagram for HSC with this aspect ratio features also
a smectic phase (that is, a phase with fluid-like layers in which the particles are nematically ordered with the director normal to the layers).
We see remnants of the latter in the active HSC system in the form of smectic clusters (not shown). These clusters usually span the system. Larger system sizes
are required to characterize these system-spanning smectic clusters, which we leave for future work. 















To understand the differences between the active nematic phase in the two models, let us first analyze the rGB system in some more detail.
In Fig. <ref>a we show the
time-averaged global order parameters from Eq. (<ref>) at packing fractions η
close to the IN phase boundary in the active rGB system.

Following a path of state points by increasing v_0 at a fixed packing fraction 
leads to a decrease in the nematic order parameter S, which 
eventually falls below our chosen threshold S_t.  

Beyond this transition point, we classify the state as isotropic
and conclude that the activity destabilizes the nematic phase. 

The transition packing fraction shifts to higher values at higher activity.
We also observe a reduced slope of S(η) at higher activity, 
which is why different threshold values S_t would result in a slightly different phase boundary.


As illustrated by the behavior of P 
in the inset of Fig. <ref>a,
the emergence of nematic order is not associated with persistent global polar ordering.

   
To make a clear statement about the behavior of the active system, it is important to discuss the role of fluctuations. 
In the global isotropic phase it is well known <cit.> 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. <ref> for the rGB ellipsoids thus corresponds to a single cluster spanning the whole system.
On increasing the system size, the associated “phase boundary” shifts to higher v_0 for a given η, 
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 <cit.>. 

In the states which we characterize as nematic the polar fluctuations are much more prominent 
than one would expect for an isotropic phase with the same parameters.
 In fact, even in the actual isotropic phase found at the same density but higher activity, the fluctuations are significantly weaker.
In Fig. <ref>b we show the time evolution of
both order parameters associated with the nematic snapshot in Fig. <ref>b. 
The pronounced temporal fluctuations near the transition result in slightly larger errors of the time-averaged values compared to the bulk phases
and also rationalize the decrease of the slope of S(η) at higher activity, observed in Fig. <ref>a.
Moving deeper into the nematic phase, the nematic order parameter can be determined quite accurately. 



In other words, we suspect that the fluctuations discussed above for the rGB ellipsoids are related to an enhancement 
of unphysical, finite-size induced self-interactions due to the persistent motion of the 
aligned rods in the nematic phase. 
However, the following considerations support our claim that the 
IN phase boundary depicted in Fig. <ref> is generic.

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 distinct 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, as indicated in Fig. <ref>b, 
whereas the average nematic order parameter is robust, i.e., the IN phase boundary in Fig. <ref> does not change.
We even found indications that the transition, i.e., the change of the nematic order parameter in Fig. <ref>a, becomes sharper in a larger system.

 For more details on the finite-size effects and fluctuations in the rGB simulations see the SI <cit.>.






The above discussion is corroborated by our simulations of the HSC system, where we do not observe significant fluctuations of the global order parameters in the nematic phase,
which is similar to the equilibrium nematic order parameter, even relatively close to the phase boundary.
In particular, the IN transition is always rather sharp, as illustrated in Fig. <ref>. This is, at least partially, due to the much larger system
size of the HSC system. The error bars, which are a measure of the standard deviation (not the standard error), in the nematic order parameter increases drastically
in the isotropic phase at higher propulsion speed, especially Pe'=1.5 and near the transition. Large error bars are an indication of large fluctuations, such as those found in the rGB system.
However, for the HSC system, the fluctuations are much less pronounced than in the smaller rGB system.
 It is known that the IN transition is first order in equilibrium, although the coexistence region is very small <cit.>. Since the region of bistability cannot completely disappear if an
infinitesimal propulsion speed is imposed, there must be a
(small) region of bistability at nonzero propulsion speed (at least for small v_0^*).
The jump in the order parameter, the magnitude of which is only weakly affected by the self-propulsion (see Fig. <ref>),
indicates that the IN transition remains discontinuous, but the expected region of bistability is smaller than our density resolution for all v_0^*
(including the equilibrium system).

Due to the different aspect ratios of the two types of particles and the resulting difference in friction,
the swimming speed at which the transition starts to shift towards higher densities is reduced.
In the isotropic phase near the IN transition, the nematic order parameter fluctuates strongly as a function of time for the larger swimming speeds.
This explains the raggedness of the curves
in Fig. <ref> for larger velocities (especially Pe'=1.5).


When investigating by eye the snapshots of the isotropic phase in the HSC system, we made the following observations:
In the isotropic phase, as for the rGB model, polar clusters were found that increase in size
when increasing activity or density (however, we made sure that the polar clusters of HSC never span the system). 
In contrast, the typical nematic configurations of the HSC do not show strong local polar ordering,
even at swimming speeds where the isotropic phase clearly exhibits large polar clusters, compare Figs. <ref>c and <ref>d.







§ CONCLUSIONS



In conclusion, we identified in 3D and for small aspect ratios a homogeneous nematic phase, close to equilibrium, 
which can be clearly distinguished from the isotropic phase, even in a relatively small system. 
 By homogeneous, we mean that there are no appreciable inhomogeneities in the local density.
This nonequilibrium nematic phase is gradually destabilized by activity
and we observe no evidence for giant number fluctuations (but we cannot exclude the possibility entirely). 
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 <cit.> for 2D is thus not universal.

 The reason for our new observations could be related to the shortness of the considered particles.
In systems of long rods, especially in 2D,
head-to-side collisions dominate and will rotate the particles 
towards either a parallel or anti-parallel orientation, 
 which can be used as an argument in favor of enhancing the nematic order rather than destroying it. 
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 <cit.>.
 The finding that the active nematic phase in 3D systems should be more stable than in 2D,
is still reasonable since in the latter case the rods have less directions in which they can escape upon a collision, so they would cluster more readily.



The investigation of the effect of low-level activity on established equilibrium states
should, of course, also be carried out with hydrodynamic interactions taken into account.
One fundamental question to be addressed is then whether there are major differences between 
 this more realistic model and our overdamped simulations. 
If momentum is conserved in such a more realistic model, is it possible that the nematic phase becomes unstable with respect to 
inhomogeneous flows with large wave length, as predicted by (linear) hydrodynamic theory <cit.>.
In our system, these large-wavelength instabilities are suppressed <cit.> by the friction and noise terms in
the equations of motion; as a result, momentum is not conserved. Similarly, walls in an experimental system also act to violate momentum conservation <cit.>.
It will be interesting to see to what extent these effects suffice to recover the behavior found in this work.
In any case, our work will provide an important benchmark to understand the role of the ignored hydrodynamic interactions in the near-equilibrium regime.
We thus hope that our work will motivate experiments on nematic phases of SPRs in near-equilibrium.
Such active liquid crystals could, for example, be constructed by rendering active a system of synthetic colloidal rods <cit.>.
We are thus confident that the problem of active perturbations of equilibrium phases, 
 pioneered by our simulations (and similar efforts for other systems <cit.>) 
is not only of pure theoretical interest. 






An open task for our overdamped simulations is to provide a
more fundamental quantitative understanding of the nature of the observed nonequilibrium IN transition
and the active nematic phase in particular.

Although the transition appears to remain of first order (as in equilibrium),
which we suspect from the sharp increase of the nematic order parameter in the HSC system,
a more careful analysis is required for a definite statement.
Since we have established a clear connection to the IN transition in equilibrium,
we do not believe that underlying mechanism driving this transition
is comparable to a liquid-vapor-like motility-induced phase separation <cit.>.
The latter (not to be confused with polar clustering) could rather be observed within the isotropic region of the phase diagram, i.e.,
at lower density, higher activity and, perhaps, only for shorter rods.

In this sense, and to obtain clarity on the collision argument, it will also be of interest to study the influence of both the aspect ratio and interparticle interactions on the active IN phase boundary. 

To properly characterize the active nematic phase and the transition region
a detailed analysis of different pair correlation functions and the orientational distribution will be presented in future work.

Along these lines, we will also explore in detail the high-density region in the HSC system to conclusively argue about the existence of an active smectic phase.











The most important open task is, however, on the theoretical side.
It would be desirable to have a first-principles theoretical 
approach to confirm our surprising predictions of the activity dependence of the IN phase boundary, even if this is limited to 
low activities, close to equilibrium. 
One obvious possibility would be to develop a linear-response theory <cit.> for an anisotropic and active system.
While the phase behavior of spherical ABPs can be explained solely by effective attractions <cit.> 
and that of active nematic rods by an effective (longer) aspect ratio <cit.>,
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 <cit.>.
This observation suggests an intuitive mapping 
to describe the IN transition in qualitative agreement with our simulations, which is yet to be quantified.
Another promising and possibly computationally efficient approach would be an implementation 
within dynamical density functional theory <cit.> for anisotropic and active systems <cit.>,
which recently has been generalized also to microswimmers in a hydrodynamic medium <cit.>.


In conclusion, there is much opportunity for further experimental, theoretical and numerical studies of the active nematic phase of SPRs.  
Beyond the bulk system, these should also address the Frank elastic behavior, the response to (time-dependent) external fields and inhomogeneous 
systems in the presence of confining walls.



§ ACKNOWLEDGEMENTS



The authors want to thank Sebastian Kapfer for a careful reading of the manuscript and valuable comments and
Hartmut Löwen for helpful discussions.
R. Wittmann and J. M. Brader acknowledge funding provided by the Swiss National Science Foundation,
M. Bott the support by the Swiss National Science Foundation through the National Center of Competence in Research Bio-Inspired Materials
 and F. Winterhalter acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the Research Unit “Geometry and Physics of Spatial Random Systems” (GPSRS) under grant number ME1361/11.
Finally, the authors gratefully acknowledge the compute resources and support provided by the Erlangen Regional Computing Center (RRZE).









§ NUMERICAL DETAILS



In this appendix we give a detailed description of the Gay-Berne (GB) model, present how to calculate the forces and torque
and describe how to extract the required order parameters from the numerical data.




 §.§ The Gay-Berne model



The Gay-Berne interaction potential for anisotropic particles is given by


    ϕ_ gb(û_1, û_2, 𝐫)    = 
       4ϵ(û_1, û_2, 𝐫̂)
     
     {( σ_0/r - σ(û_1, û_2, 𝐫̂) + σ_0)^12
       - ( σ_0/r - σ(û_1, û_2, 𝐫̂) + σ_0)^6} ,


With the unit vectors û_1 and û_2 specifying the orientation of the interacting particles 1 and 2 and 𝐫, r, 𝐫̂ their center to center-vector, -distance and -direction.
The attraction depth ϵ(û_1, û_2, 𝐫̂) and the range σ(û_1, û_2, 𝐫̂) of the particle interaction are dependent on the orientation. 

The shape of Gay-Berne particles is defined through an anisotropy parameter


    χ = 1/κ^2-1/1/κ^2+1,     κ = σ_s/σ_e.


Here σ_e is the “length” of the particle defined by the end-to-end interaction and σ_s the “width” of the particle defined though the side-to-side interaction.
For infinitely long cigar shaped particles (σ_e →∞) the anisotropy parameter χ→ 1, in contrast for infinitely thin oblate-like particles (σ_s → 0) we have χ→ -1.

The orientation dependent interaction range is given by


    σ(û_1, û_2, 𝐫̂)    = σ_0 [ 1 - 1/2χ{(𝐫̂·û_1 + 𝐫̂·û_2)^2/1 + χ(û_1 ·û_2)
       + (𝐫̂·û_1 - 𝐫̂·û_2)^2/1 - χ(û_1 ·û_2)}]^-1/2.


The well depth is defined as


    ϵ(û_1, û_2, 𝐫̂) = ϵ_0  ϵ(û_1, û_2)^ν ϵ'(û_1, û_2, 𝐫̂)^μ,



    ϵ(û_1, û_2) = { 1 - χ^2  (û_1, û_2)^2 }^-1/2,



    ϵ'(û_1, û_2, 𝐫̂)    = 1 - χ'/2{(𝐫̂·û_1 + 𝐫̂·û_2)^2/1 + χ'(û_1 ·û_2)
       + (𝐫̂·û_1 - 𝐫̂·û_2)^2/1 - χ'(û_1 ·û_2)},


with the parameter χ' describing the anisotropy in the well depth:


    χ' = 1 - κ'^1/μ/1 + κ'^1/μ,     κ' = ϵ_s/ϵ_e.


In the limit of spherical particles i.e. σ_e = σ_s, ϵ_e = ϵ_s one finds χ = 0 and χ' = 0 and therefore the Gay-Berne potential becomes a regular Lennard-Jones interaction.

In our study we employed a soft repulsive WCA-like version of the Gay-Berne interaction, which is obtained by shifting and truncating the Gay-Berne potential:


    ϕ_ wca   (û_1,û_2,𝐫) 
       =ϕ_ gb(û_1,û_2,𝐫) + ϵ(û_1,û_2,𝐫)    r<r_ min
    
        0                            r ≥ r_ min,


where r_ min(û_1,û_2,𝐫̂) = 2^1/6σ_0 
+ σ(û_1,û_2,𝐫̂) - σ_0 is the 
minimum of the Gay-Berne interaction.



 §.§ Calculation of Force and Torque


The expression to calculate the force and torque we denote for the regular Gay-Berne interaction (<ref>). 
For this we introduce the scaled variable


    R = r - σ(û_1, û_2, 𝐫) + σ_0 /σ_0


and the interaction potential can be written as


    ϕ_ gb(û_1, û_2, 𝐫) = 4ϵ(û_1, û_2, 𝐫̂) {(1/R)^12 - (1/R)^6}.


Furthermore we define a function g(X) for the relative orientation of two interacting particles


    g(X) = 1 - X/2{(𝐫̂·û_1 + 𝐫̂·û_2)^2/1 + X (û_1 ·û_2) + (𝐫̂·û_1 - 𝐫̂·û_2)^2/1 - X (û_1 ·û_2)}.


Hence we have


    σ(û_1, û_2, 𝐫) = σ_0   g(χ)^-1/2


and

    ϵ'(û_1, û_2, 𝐫̂)  = g(χ').


When changing the distance vector 𝐫 between two particles the interparticle vector 𝐫̂ changes as well, which can be made explicit in g(X)


    g(X) = 1 - X/2r^2{(𝐫·û_1 + 𝐫·û_2)^2/1 + X (û_1 ·û_2) + (𝐫·û_1 - 𝐫·û_2)^2/1 - X (û_1 ·û_2)}.




  §.§.§ Force


The force is given by 𝐅 = -∇ϕ_ gb. We here denote the expression for the force in x-direction:


    -F_x  = ∂ϕ_ gb/∂ x   = 4 ∂ϵ/∂ x{(1/R)^12 - (1/R)^6}
       + 4ϵ{(6/R)^7 - (12/R)^13}∂ R/∂ x,


with

    ∂ϵ/∂ x = ϵ_0 ϵ(û_1, û_2)^νμϵ'(û_1, û_2, 𝐫̂)^μ-1∂ g(χ')/∂ x,



    ∂ R/∂ x = 1/σ_0(x/r - ∂σ/∂ x),


and


    ∂σ/∂ x = -1/2σ_0 g(χ)^-3/2∂ g(χ)/∂ x.


Finally the derivative of the orientation function g(X) is given by


    ∂ g(X)/∂ x   = x X/r^4{(𝐫·û_1 + 𝐫·û_2)^2/1 + X (û_1 ·û_2) + (𝐫·û_1 - 𝐫·û_2)^2/1 - X (û_1 ·û_2)}
       - X/r^2{(𝐫·û_1 + 𝐫·û_2)/1 + X (û_1 ·û_2) (û_1^x + û_2^x) 
       + (𝐫·û_1 - 𝐫·û_2)/1 - X (û_1 ·û_2) (û_1^x - û_2^x) }.


The force in y- and z-direction can be calculated equivalently.



  §.§.§ Torque


Due to the angular dependence of the Gay-Berne potential particles experience torque.
So far we only determined the center to center force.
We can calculate the torque from an equivalent force 𝐄 acting on a point at unit distance from the center of the particles.
This equivalent force can be calculated from the derivative of the potential with respect to the unit vector


    𝐄 = -  (
    		[ ∂ϕ_ gb/ ∂û_1^x; ∂ϕ_ gb/ ∂û_1^y; ∂ϕ_ gb/ ∂û_1^z ] )


Again we denote the derivatives of ϕ_ gb with respect to û_1^x, but in other directions and for particle 2 one obtains equivalent results.


    -E_x = ∂ϕ_ gb/∂û_1^x   = 4 ∂ϵ/∂û_1^x{(1/R)^12 - (1/R)^6}
       + 4 ϵ{6/R^7 - 12/R^13}∂ R/∂û_1^x.


Where


    ∂ϵ(û_1, û_2, 𝐫̂)/∂û_1^x   =   ϵ_0 νϵ(û_1, û_2)^ν-1∂ϵ(û_1, û_2)/∂û_1^xϵ'(û_1, û_2, 𝐫̂)^μ
       +   ϵ_0 νϵ(û_1, û_2)^νμϵ'(û_1, û_2, 𝐫̂)^μ-1∂ g(χ')/∂û_1^x,



    ∂ϵ(û_1, û_2)/∂û_1^x = ϵ(û_1, û_2)^3 χ^2 û_1 ·û_2 û_2^x,



    ∂ R/∂û_1^x = 1/2( σ(û_1, û_2, 𝐫)/σ_0)^3 ∂ g(χ)/∂û_1^x.


And the derivative of the orientation function is given by


    ∂ g(X)/∂û_1^x = 
       -X/2 [ r̂_x {2(𝐫̂·û_1 + 𝐫̂·û_2)/1 + X (û_1 ·û_2) + 2(𝐫̂·û_1 - 𝐫̂·û_2)/1 - X (û_1 ·û_2)}
       + X û_2^x {(𝐫̂·û_1 - 𝐫̂·û_2)^2/(1 - X (û_1 ·û_2))^2 - (𝐫̂·û_1 + 𝐫̂·û_2)^2/(1 + X (û_1 ·û_2))^2}].


The second term in the first bracket changes sign when taking the derivative with respect to orientation of particle 2.
Finally we obtain the torque by the cross product of 𝐄 and the orientation vector û_1,


    T = û_1 ×𝐄.




 §.§ Calculation of order parameters
 


In this appendix we provide information on how to extract the order parameters from the simulation data.
The orientational behavior of an ensemble of N anisotropic particles can be analyzed using an order parameter S, which is defined as


    S = ∑_i=1^N 3 cos^2β_i - 1/2N,


where β_i is the angle between the orientation-vector of particle i and the nematic director (unit vector indicating the mean orientation of the particles).
The order parameter can take values between 0 and 1, where S=0 indicates that the system is in a fully isotropic state with random orientation and S=1 means perfect alignment of the particles.
However in simulations the nematic director is not known a priori. Following reference <cit.> we consider a tensorial order parameter 


    Q_αβ =  1/N∑_i=1^N 3/2û_iαû_iβ - 1/2δ_αβ,     α,β = x,y,z.


This 2nd-rank tensor can immediately be computed from the single particle orientations û. 
It has three eigenvalues, of which the largest is the order parameter and the corresponding eigenvector the nematic director 𝐧̂.

In systems with polar order, the nematic order parameter is also nonzero and the nematic director lies along the direction of polar order.
Thus, we can reuse the nematic director to calculate the polar order parameter, when we define it as


    P = 1/N| ∑_i=1^N û_i ·𝐧̂ |.


Again the polar order parameter ranges between 0 and 1, indicating no polar alignment and perfect polar ordering, respectively.



