Collective dynamics and pair-distribution function of active Brownian ellipsoids

The gist The system exhibits a rich state behavior that includes cluster formation, local polar order, polar flocks, and disordered states via computer simulations.

State Diagram Investigation

The study obtains the full state diagram of active Brownian ellipsoids depending on the P´eclet number and packing density via computer simulations. The system is found to exhibit a rich state behavior that includes cluster formation, local polar order, polar flocks, and disordered states. The state diagram of active ellipsoids is more complex than that of standard active Brownian particles (ABPs) since the anisotropic particle shapes allow for torques. Typically, motility-induced phase separation (MIPS) is suppressed in systems of active ellipsoids compared to active spheres due to the torque interactions. When undergoing MIPS, spherical particles collide and hinder each other’s movement. The torque resulting from particle-particle interaction can suppress aggregation at an existing cluster’s surface for nonspherical particles such as rods and ellipsoids moving along their long axes. Even minor anisotropies, like a length-to-diameter ratio of an elliptic particle such as κ = ael/bel = 1.0424 (with the length ael, the diameter bel, and the length-to-diameter ratio κ) can cause MIPS aggregates to dissolve. However, torque can also enhance cluster formation. This effect was reported by Zhang et al. [50], who used spherical particles with an anisotropic interaction where repulsion is stronger at the rear of the particles than at the front. While the torque resulting from shape interaction inhibits MIPS, it can lead to the formation of other phases. Increasing the length-to-diameter ratio further can lead to the emergence of local nematic order while the system is in global disorder. For even higher ratios (such as κ = 1.96), highly polarized domains with local nematic order emerge. Finally, extremely high anisotropies allow for the formation of a polar band with local smectic order. Jayaram et al. [16] found that MIPS clusters become polar domains and polar bands for increasing anisotropies in the case of infinite P´eclet numbers. Due to their movement, active particles collide; therefore, particle interaction hinders their movement. If the particles are elongated, they can slide past each other more easily and do not hinder each other’s movement that much, which reduces the trapping effect that causes phase separation. Rather than a persistent large single cluster as in the case of active spheres, one therefore observes a continuous appearance and disappearance of small clusters. This happens even for small anisotropies such as κ = 1.15 and in the absence of diffusion (infinite P´eclet number). Further increasing the anisotropy, another state transition occurs, and a single, highly polarized cluster emerges.

Equations of Motion and Dimensionless Variables

The general equations of motion for uniaxial ABPs in two spatial dimensions that can be influenced by interactions, external fields, and orientation-dependent propulsion read r˙ i = DTkBTFAubi + Fint,i(rj,ubj) + Fext(ri) + ξT, (1) and φ˙ i = DRkBTMint,i(rj,ubj) + ξR (2). As ellipsoids are anisotropic, their translational Brownian motion is also anisotropic. The translational Brownian noise ξT is implemented by a random force with zero mean. The variances of its components k and l are given by ⟨ξT,k(t), ξT,l(t′)⟩ = 2DT,kl(φ(t))δ(t − t′) where k, l ∈ x, y and δ(t− t′) is the Dirac delta function. The rotational Brownian noise is a random torque with zero mean, which variance is given by ⟨ξR(t), ξR(t′)⟩ = 2DRδ(t − t′). In this work, we focus on ellipsoids of revolution that have the same volume as a sphere with diameter σ0. The length of their axis of revolution, which is also referred to as the polar axis, is ael = 41/3 σ0. The axes perpendicular to the polar axis, which will be referred to as the equatorial axes, have a length of bel = 2−1/3 σ0. This means that the ellipsoids have a length-to-diameter ratio κ = 2 and the same volume as a sphere of the same diameter. Brownian particles undergo diffusive motion. While translational diffusion is typically isotropic for spheres, it is anisotropic for ellipsoids. The diffusive motion depends on the diffusion constant, which is a scalar for spheres and a matrix for anisotropic particles. For uniaxial particles, such as ellipsoids, there are two different translational diffusion coefficients depending on the two hydrodynamic friction coefficients γT,∥ and γT,⊥. The friction coefficients correspond to the friction for movement parallel and perpendicular to the particles’ polar and equatorial axes, respectively. Using the Einstein-Smoluchowski equation [55–57], these two friction coefficients lead to two diffusion coefficients. The diffusion matrix DT reads DT = DT,∥(ub ⊗ ub) + DT,⊥(I − ub ⊗ ub), (3) where I is the two-dimensional identity matrix, DT,∥ = kBT /γT,∥ the parallel diffusion coefficient, DT,⊥ = kBT /γT,⊥ the perpendicular diffusion coefficient.

Lennard-Jones Units and Dimensionless Variables

In the simulations, we use Lennard-Jones units. Distance and energy are measured as multiples of the distance σ0 and the energy ϵ (appearing in the interaction potential, see Appendix A) respectively. Friction coefficients are multiples of the friction coefficient γT,sp of a sphere with diameter σ0. Thus, the Lennard-Jones time is given by τLJ = σ 20γT,sp/ϵ = σ 20kBT / (ϵDT,sp) where DT,sp is the diffusion coefficient of a sphere that has the same volume as the ellipsoid. All parameters are measured in this unit system. For example, the active propulsion force FA is given by FA = 24ϵ/σ0. As in previous work [28, 44], we use two dimensionless variables to parametrize the state of the system, namely the P´eclet number for ellipsoids Pe and the global packing density Φ0. The P´eclet number is (following Ref. [22]) defined as Pe = v0aelDR, (5) with the propulsion speed v0 = FADT,∥/(kBT) of the ellipsoids, which defines Pe via the ratio of persistence length v0/DR and particle size ael. In the limit of spherical particle shapes, Eq. (5) corresponds to the standard definition of the P´eclet number for active spheres, except for a factor 3. The P´eclet number can be changed by tuning either the temperature T or the propulsion force FA. We chose to change the P´eclet number by changing T as we would need to cope with a different effective size of the particles due to stronger propulsion when tuning the propulsion force. The global packing density is defned as Φ0 = NpartApart/Asim with the number of particles in the system Npart, the area of a particle Apart = aelbelπ, and the area of the simulation domain Asim.

Simulation Details and State Classification

Since the influence of the particle shape on the collective dynamics has been studied thoroughly in previous work [21–23], we focus on the impact of Pe and Φ0 for a fixed length-to-diameter ratio κ = 2. The state diagram is of interest by itself, but is also required to find the parameter combinations for which the particle distribution stays homogeneous (which will later be used for obtaining the pair-distribution function). Thus, we performed computer simulations of active ellipsoids by numerically solving Eqs. (1) and (2) for different P´eclet numbers Pe and global packing densities Φ0. Then, we classified the observed states. The simulations were performed with a modified version of the molecular dynamics simulation package LAMMPS [58]. We modified the software package such that overdamped dynamics could be simulated. For integrating the equations of motion, we used the EulerMaruyama method. The simulations start with random initial positions in a quadratic simulation domain with a side length of 128σ0 and periodic boundary conditions.

1. The improved system can accurately predict the full state diagram of active Brownian ellipsoids depending on the P´eclet number and packing density via computer simulations, allowing for classification into states such as cluster states, local polar order, polar flocks, and disordered states.

2. The system can obtain a quantitative description of collective behavior by numerically investigating the pair-distribution function and obtaining an analytical representation that can be used in further theoretical work.

3. The improved system can perform quantitative analysis by calculating the reduced interaction energy per particle Eint/(ϵN), the average local polarization divided by density ⟨∥P ∥⟩/Φ0, and the average nematic order divided by density ⟨S⟩/Φ0 to distinguish between states.

4. The system can accurately model inter-particle interactions using a modified version of the molecular dynamics simulation package LAMMPS [58] with an overdamped dynamics implementation via the Euler-Maruyama method to simulate active ellipsoids.

5. The system can generate a predictive analytical approximation of the pair-distribution function, gapprox, by fitting its Fourier coefficients (up to third order) using empirical fit functions derived from simulation data and interpolated across parameter space via h(Pe, Φ0) = X2mX3n0 Pe m 2 Φ n0 um,n.