{'fraction_non_alphanumeric': 0.03373015873015873, 'fraction_numerical': 0.013741916519694298, 'mean_word_length': 4.687003760969494, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Crystallization in a dense suspension of self-propelled particles

Using Brownian dynamics computer simulations we show that a two-dimensional suspension of self-propelled ("active") colloidal particles crystallizes at sufficiently high densities. Compared to the equilibrium freezing of passive particles the freezing density is both significantly shifted and depends on the structural or dynamical criterion employed. In non-equilibrium the transition is accompanied by pronounced structural heterogeneities. This leads to a transition region between liquid and solid in which the suspension is globally ordered but unordered liquid-like"bubbles"still persist.

Using Brownian dynamics computer simulations we show that a two-dimensional suspension of self-propelled ("active") colloidal particles crystallizes at sufficiently high densities. Compared to the equilibrium freezing of passive particles the freezing density is both significantly shifted and depends on the structural or dynamical criterion employed. In non-equilibrium the transition is accompanied by pronounced structural heterogeneities. This leads to a transition region between liquid and solid in which the suspension is globally ordered but unordered liquid-like "bubbles" still persist. Recently, the collective dynamics of self-propelled ("active") particles has become a topic of intense research [1,2] resulting in a wealth of new non-equilibrium phenomena like swarming [3,4], clustering [5][6][7] and active swirling [8]. These phenomena have been observed both in dense bacterial solutions [9] and in artificial microswimmers [10]. Excellent model systems for selfpropelled particles are colloidal suspensions, where the motility of colloidal particles can be achieved and steered by magnetic beads acting as artificial flagella [11], by catalytic reactions at Janus-particles [12], or by laser-heated metal-capped particles [13].
The purpose of the present Letter is to show that selfmotile interacting colloidal particles in two dimensions still freeze into a crystalline lattice displaying long-ranged orientational order despite the fact that energy is injected incessantly. We explore the nature of this nonequilibrium transition by Brownian dynamics computer simulations of a Yukawa model of self-propelled particles. We use a minimal model without explicit alignment of particle orientations. In equilibrium, i.e. in the absence of self-propagation, freezing and melting of colloidal suspensions is well understood [14]. But even for passive particles it is known that freezing is seriously affected and changed under non-equilibrium conditions, e.g. in a timeoscillatory external force field [15] or in shear flow [16]. Recently it has also been shown that active matter can reach steady states with frozen fluctuations [17].
For self-propelled particles we find that the freezing transition is largely shifted relative to its equilibrium location. This shift cannot be explained by a simple scaling using the concept of an effective temperature [18]; quite in contrast to sedimentation profiles of suspensions [19] or the long-time diffusion of single propelled particles [20]. Rather, the transition points based on different criteria for melting and freezing, which agree in equilibrium, diverge. In particular, the dynamical Lindemann-like melting [21,22] and freezing criteria [23,24] define a transition region between liquid and solid characterized by inhomogeneities of the orientational order parameter.
We study a suspension of N self-propelled particles moving in two dimensions and immersed in a solvent.
Even though the particles are driven we assume that the solvent remains in equilibrium at the well-defined temperature T . The overdamped motion of the ith particle is described througḣ The noise ξ i models the stochastic interactions with the solvent molecules. It has zero mean and correlations , where 1 is the identity matrix. Throughout the paper we employ dimensionless quantities and measure energy in units of k B T , length in units of ρ −1/2 , and time in units of (ρD 0 ) −1 . Here, ρ is the number density and D 0 is the bare diffusion coefficient. Particles interact pairwise through the repulsive Yukawa potential with screening length λ −1 and dimensionless coupling parameter Γ ≡ V 0 √ ρ/k B T , where V 0 is the bare potential strength. The total potential energy then becomes U = i<j u(|r i − r j |). In addition to the conservative force due to U a constant force f propels every particle in the direction In the minimal model studied here we assume that these particle orientations undergo free diffusion without explicit alignment. For spherical particles with diameter σ the rotational diffusion coefficient is We perform Brownian dynamics simulations for N = 1936 particles using periodic boundary conditions. Commensurable box dimensions L x /L y = 2/ √ 3 are chosen such that the suspension can crystallize into the hexagonal crystal without any defects. We fix the rotational diffusion coefficient to D r = 3.5 and the inverse screening length to λ = 3.5; leaving Γ and f as variable parameters. The time step for updating particle positions is ∆t = 10 −4 , while particle-particle interactions are cut off after an inter-particle distance of 7/λ = 2. We simulate cooling and melting runs for forces 0 f 15. For the cooling runs we use one long trajectory. We start from a random particle configuration with uniformly distributed orientations. After a sufficient large relaxation time (t = 25) we collect data for 500 time units. The coupling parameter Γ is then increased by 20 and the protocol of relaxation and recording data is repeated until we reach the maximal Γ. The melting runs for each pair {Γ, f } are independent starting out of the perfect hexagonal crystal albeit with random particle orientations. Again, we wait an adequate amount of time before collecting data for 50 time units. For both cooling and melting we record data from 5 independent runs for each {Γ, f }.
We monitor structural changes through the global bond-orientational order parameter [25] where N (i) is the set of the six nearest neighbors of the ith particle and θ ij is the angle between the bond vector pointing from particle i to j and an arbitrary fixed axis. This order parameter is practically zero in the disordered phase whereas in a perfect crystal ψ 6 = 1. In Fig. 1A we plot the ψ 6 values averaged over all runs for both the cooling and the melting protocol. There is a clear transition between a disordered liquid and an ordered crystalline phase even for self-propelled particles (f > 0). However, while the transition is rather abrupt for f = 0 the structural ordering is more gradual for higher propelling forces f . In Fig. 1A we do not resolve a possible hexatic intermediate phase [26]. However, we note that no hysteresis is observed in agreement with a second-order transition scenario. As a structural criterion for both the melting and freezing transition we determine Γ * S from the condition ψ 6 = 0.45. In particular, for f = 0 we find Γ * S 240, which agrees well with a previous estimate [27].
Cooling the suspension, a dynamical criterion for freezing is given by the precipitous drop of the long-time diffusion coefficient with ∆r i (t) ≡ r i (t) − r i (0). In Fig. 1B we plot the diffusion coefficient for different forces. The value Γ * D at which the suspension freezes is estimated from the condition D = 0.086 [23]. This gives an upper bound Γ < Γ * D to the liquid region, see the phase diagram Fig. 1C. Moreover, for not too large forces Γ * D Γ * S correlates well with the position of the structural ordering as observed in Fig. 1A. Hence, this dynamical criterion for freezing based on particle mobility extends only to weakly driven suspensions of self-propelled particles. Note that at large forces and small Γ the diffusion coefficient D exceeds 1, the diffusion coefficient of a free passive Brownian particle.
We next consider a dynamical criterion for melting starting in the solid state and decreasing Γ. In one of the first theories for melting Lindemann conjectured that melting is caused by atom vibrations that start to interpenetrate [28]. In our case it is natural to consider the vibrational displacements of particles with respect to their lattice position. The Lindemann criterion then states that melting commences once the vibrational displacements reach a certain fraction of the lattice spacing. However, in two dimensions fluctuations on long wavelengths eventually destroy long-ranged positional order in the crystal [29]. The mean-squared displacement, therefore, is not a good measure to distinguish the liquid from the crystal. Instead, one defines a Lindemann-like parameter [21,22] from the neighbor-neighbor displacements. Here, i and j denote two particles that are initially neighbors. The lat- tice spacing of the hexagonal crystal is ≡ 2 1/2 3 −1/4 1.075. In the liquid γ L (t) diverges for long times without a plateau, whereas in the solid one observes a well defined plateau with Lindemann parameter γ L . Hence, we determine the melting point from the smallest value Γ * L for which we still observe a plateau with value γ * L , see Fig. 2A. Above Γ > Γ * L the suspension is crystalline both with respect to orientational order and a vanishing diffusion coefficient.
Using a simplified picture to describe the process of melting we assume particles to move independently close to their lattice position. The linearized forces then read −∇ i U ≈ −k(r i − r 0 i ) with effective curvature k ∝ Γ. The initial positions r i (0) = r 0 i correspond to lattice positions in the hexagonal crystal. A straightforward calculation of Eq. (6) in the limit t → ∞ yields In Fig. 2B the plateau value γ * L as a function of force is plotted together with the prediction γ L (f, Γ * L ) from Eq. (7). Both values show excellent agreement. For the plot we have fixed the proportionality between k and Γ such that the values for f = 0 are equal. Moreover, γ * L (0) 0.026 agrees well with previous experiments [22]. While structural and dynamical criteria agree in equilibrium the phase diagram Fig. 1C shows that in nonequilibrium there is a transition region Γ * D < Γ < Γ * L between liquid and crystal, which widens for larger forces f . This region of parameter space is characterized by a high structural order but non-vanishing long-time diffusion. Moreover, the dynamical freezing and melting lines do not follow the orientational order but are shifted to higher Γ at higher forces. This implies that at high propelling speeds structural ordering occurs before dynamical freezing. While an effective temperature could be defined individually for each criterion, the resulting values as a function of force clearly do not agree. To gain further insight we define the order parameter q 6 (i) ≡ Re 1 6 j∈N (i) q 6 (i)q * 6 (j) (8) per particle in order to effectively describe the local environment of every particle [30]. The advantage of the neighbor-shell averaging compared to |q 6 | 2 is thatq 6 more sharply distinguishes liquid-like from ordered regions. In Fig. 3 probability distributions ofq 6 for f = 0 and f = 8 are plotted. For the crystal (ψ 6 0.8) no difference between the driven and the undriven suspension is discernible (see also the last row of Fig. 3). In the liquid (ψ 6 0) the driven suspension is somewhat less structured compared to equilibrium. This is caused by the larger effective diffusion due to the propulsion. A large difference can be seen in the distributions for suspensions with ψ 6 0.45, i.e., in the transition regime. Here the driven suspension is locally more ordered but with a long tail that extends down to unordered particles. The spatial distribution of order and disorder corresponding to this ψ 6 for a single snapshot is shown in the second row of Fig. 3. For f = 8 the suspension is overall more ordered but also more heterogeneous, i.e., small, well separated liquid "bubbles" remain. Interestingly, the ψ 6 = 0.67 iso-line crosses the melting line such that for f = 8 it is within the transition region. Two snapshots for this case are depicted in the third row of Fig. 3. Due to the crossing the two forces now also describe two different dynamic regimes: while diffusion has effectively ceased in equilibrium, some particles still move in the driven suspension.
In conclusion, we have shown by using Brownian dynamics computer simulations that self-motile colloidal particles crystallize at sufficiently high densities. As compared to the equilibrium freezing of passive particles there is a significant shift in the freezing density and additional large structural fluctuations appear caused by the self-propulsion. In principle, our predictions are verifiable in real-space experiments on colloidal model swimmers on a (quasi) two-dimensional substrate [12,13].
In future work, it would be interesting to generalize our model to one which embodies an explicit swarming behavior such as a self-propelled rod model [6]. Furthermore, since equilibrium freezing is different in two and three spatial dimensions, it would be very interesting to simulate a corresponding three-dimensional model. Last but not least the influence of self-motility of the glass transition has not yet been studied. Since glass formation competes with crystallization and is typically accompanied with dynamical heterogeneity [31][32][33], selfpropulsion may introduce an internal source of additional fluctuations which can help to form amorphous structures provided the density is large enough.