Role of topological defects in the two-stage melting and elastic behavior of active Brownian particles

We find that crystalline states of repulsive active Brownian particles at high activity melt into a hexatic phase, but this transition is not driven by an unbinding of bound dislocation pairs as suggested by the Kosterlitz-Thouless-Halperin-Nelson-Young theory. Upon reducing the density, the crystalline state melts into a high-density hexatic state devoid of any defects. Decreasing the density further, the dislocations proliferate and introduce plasticity in the system, nevertheless maintaining the hexatic state, but eventually melting into a fluid state. Remarkably, the elastic constants of active solids are equal to those of their passive counterparts, as the swim contribution to the stress tensor is negligible in the solid state. The sole effect of activity is that the stable solid regime shifts to higher densities. Furthermore, discontinuities in the elastic constants as a function of density correspond to changes in the defect concentrations rather than to the solid-hexatic transition.


I. INTRODUCTION
According to the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory, a two-dimensional solid of passive particles melts via a continuous transition into an intermediate hexatic state of quasi-long-range bond orientational order and melts subsequently via a second continuous transition into a fluid state [1][2][3]. These transitions are triggered by the unbinding of topological defects, which are particles with a nonconforming number of neighbors with respect to that of the crystal lattice. The coordination number is 6 for an ideal triangular lattice. Hence, particles with a number of neighbors N b that deviates from 6 are classified as defects. One can distinguish defects that either exist freely as fivefold and sevenfold disclinations or cluster into dislocations (5-7 pairs), dislocation pairs (5-7-5-7 quartets), or higher-order clusters. The debate on the melting behavior of an equilibrium system of hard disks as well as short-range repulsive disks was only settled a few years ago, in both simulations [4][5][6][7][8] and experiments [9], showing a two-stage melting scenario that deviates from the KTHNY scenario of a continuous solidhexatic transition driven by an unbinding of dislocation pairs and a first-order hexatic-fluid transition due to a proliferation of grain boundaries.
Remarkably, nonequilibrium systems of self-propelled particles [10][11][12][13], which constantly convert energy from the environment into persistent motion, have recently also been shown to follow such a two-stage melting behaviour [14,15]. It was found that the first-order nature of the liquid-hexatic transition persists up to a small degree of activity. The transition then becomes continuous until it reappears as a coexistence of dilute and dense states at high activity [14,15]. The solid state was found to melt into a hexatic state via a continuous transition. In this work we further explore this hexatic-solid transition with respect to the role of topological defects discussed in the KTHNY theory and investigate whether the melting is driven by transitions in defect concentrations.
We numerically simulate a system of N active Brownian particles exhibiting overdamped Langevin dynamics where v 0 e i is the self-propulsion speed, γ is the damping coefficient, k B is the Boltzmann constant, and T is the temperature of the solvent. The particle orientation e i = (cos θ, sin θ ) undergoes free rotational diffusionθ i = √ 2D r r i , where D r is the rotational diffusion coefficient. The quantities t i and r i are unit-variance Gaussian noise terms with zero mean. The particles interact with a pairwise repulsive Weeks-Chandler- 6 ] + ε with a cutoff at r c = 2 1/6 σ , and we set k B T /ε = 1. We first determine the phase boundaries by measuring the equation of state (pressure-density curves), density histograms, and the decay of the orientational and positional correlation functions. We present the state diagram in Fig. 1 in the activity-density (Peρ) representation, where the Péclet number is defined by Pe = v 0 τ /σ as the ratio of the persistence length of self-propelled motion to the particle diameter and τ the timescale of thermal diffusion. The state diagram shows that the continuous hexatic-solid transition persists upon introducing activity, but shifts to higher densities due to the softness of the interaction potential. We hereby assume that the solid state transforms into a hexatic phase when the positional correlations decay with a power-law exponent η T 1/3 as described by the equilibrium KTHNY theory.

II. TOPOLOGICAL DEFECTS
We identify the topological defects by performing a Voronoi construction and calculating the number of neighbors N b for each particle. We classify them as N b -fold defects and account for fivefold and sevenfold defects. In Fig. 2 we show typical configurations, highlighting these defects in the fluid, hexatic, and crystalline states at Pe = 0, 4.8, and 71.5. At low activity Pe = 4.8, defect configurations are observed similar to those for the passive hard-disk systems [6]. A finite number fraction of dislocation pairs can be identified in the solid state which can exist due to thermal fluctuations without disturbing the positional order. In the hexatic state, the presence of unpaired dislocation defects causes the positional order to decay exponentially, but the orientational order decays algebraically. Finally, the fluid state comprises many defect clusters and fivefold and sevenfold disclinations, both destroying the local bond orientational order. However, the hexatic states at Pe = 71.5 are significantly different in defect configurations. The low-density hexatic state at ρσ 2 = 1.420 shows a high number fraction of dislocations with almost no dislocation pairs, whereas the high-density hexatic state (ρσ 2 = 1.600) at the same Pe shows a complete absence of any defects.
FIG. 2. Typical sections of active Brownian particle configurations (80σ × 80σ ) for Pe = 0.0, 4.8, and 71.5 at labeled states, showing fivefold (blue) and sevenfold (red) defects, other defects (black), and particles with N b = 6 (gray). Some vacancies are indicated by arrows in the hexatic and solid states. The configurations at Pe = 4.8 are similar to the corresponding passive states. For Pe = 71.5 we find that with increasing density the hexatic states show a decrease in the number fraction of defects, with a complete absence of defects at a density ρσ 2 1.560. The positional correlations become quasi-long-range around a density ρσ 2 1.900 for Pe = 71.5 as shown in Fig. 1. The insets in the top right and bottom right corners show the same configurations as the main panels but colored according to the hexatic and positional order parameters ψ 6 and ψ T , respectively, following the color mapping shown at the top.   We further quantify this observation by measuring the number fractions of the different types of defects, which are plotted in Fig. 3(a) as a function of density for various Pe. We clearly observe that the overall number fraction of defects N total /N increases upon decreasing the density for all Pe. At low Pe 2.4, we find that the solid phase contains mainly bound dislocation pairs, which move freely in the solid phase. The number fraction of bound dislocation pairs N quart /N increases slightly upon approaching the hexatic-solid transition, whereas the number fraction of dislocations N pair /N increases more rapidly. At the hexatic-solid transition N quart /N is still higher than N pair /N, but at lower density this scenario is reversed. The high fraction of N pair /N implies that the quasi-long-range positional order of the system is destroyed, suggesting that the solid-hexatic transition is induced by the unbinding of dislocation pairs. Upon reducing the density further, we find that the fraction N free /N of fivefold and sevenfold defects also starts to increase due to the unbinding of dislocations into disclinations and the bond orientational order decays exponentially in the stable liquid phase. To summarize, we find that the two-step melting scenario consisting of a solid-hexatic transition driven by the unbinding of dislocation pairs and a hexatic-fluid transition caused by defect clusters as observed for the two-dimensional (2D) passive systems of short-range repulsive particles persists at low activity.
For higher activity Pe 7.2, the behavior of dislocations is different from that for passive systems as shown in Fig. 2. For Pe 7.2 the system seems to support a higher fraction of dislocations as compared to dislocation pairs for all densities, as observed in Fig. 3(a), where N quart never exceeds N pair . This observation is more clearly visible in Fig. 3(b), where we plot the differences (N pair − N quart )/N and (N pair − N free )/N. We find that for Pe 2.4 and for high densities (N pair − N quart )/N is negative whereas for Pe 7.2 this difference is always positive, thereby demonstrating that a higher fraction of isolated dislocations over bound ones is favored at higher activity.
In Figs. 3(a) and 3(b) we also mark the hexatic-solid transition densities as identified from the spatial decay of positional correlations [16]. These densities agree closely with the minimum in (N pair − N quart )/N which corresponds to a reversal in the trend of defect concentrations and confirms our previous observation that for Pe 2.4 the hexatic-solid transition is driven by the unbinding of dislocation pairs. The difference (N pair − N free )/N, however, is always positive for all densities and shows similar trends for all values of Pe considered.
In Fig. 3(a) we also locate the densities above which we find a complete absence of defects in the system for all values of Pe, marked by a cross on the x axis and the corresponding density range by crosshatching. As we increase the activity, this point crosses over from the dense crystal state to the hexatic state. For Pe = 23.8, 47.7, and 71.5 we clearly see that a large region corresponding to the dense hexatic states is free from any kind of topological defects in contrast to the low-activity systems. This dense hexatic region, devoid of any defects, still has an exponentially decaying positional order as indicated by the background colors corresponding to the state diagram (Fig. 1). The activity-induced fluctuations decorrelate the particle positions at long range and are responsible for a faster decay of positional order in this density regime. The solid-hexatic transition in active systems is thus driven by a striking nonequilibrium feature. We now test the correspondence of changes in defect concentrations and the elastic response of the system with respect to the predictions of the KTHNY theory.

III. ELASTIC MODULI
The KTHNY theory for the melting of 2D equilibrium solids is based on the linear elastic properties of a continuum. Although the equilibrium theory relies on the elastic deformation energies to resolve the transition to a fluid state with a vanishing shear and renormalized Young's modulus, we can extend the notion of mechanical stress, which is well defined for an isotropic active fluid, to describe the elastic behavior in terms of the response to an externally imposed linear strain on the simulation box. In our simulations, we start from a perfect hexagonal initial configuration with N = 2.8 × 10 3 particles and measure the full stress tensor P αβ , comprising the ideal, virial, and swim components [17], in the deformed box due to a fixed small linear strain xx ∈ [−0.01, 0.01]. We calculate the Lamé elastic coefficients λ and μ from the effective stiffness tensor B obtained from the slope of a linear fit to the stress versus strain curves [18][19][20] (see Ref. [16] for details).
In Fig. 4(a) we plot the Young's modulus K = 4μ(λ + μ)/(λ + 2μ) obtained for 0.0 Pe 71.5 and identify the densities where K shows a discontinuous transition. Interestingly, we observe that K collapses onto a single master curve for all Pe. This remarkable result can be explained from the fact that the swim contribution P swim to the stress tensor [16] is negligible in the solid, and hence the elastic moduli for active solids equal the ones of their passive counterparts at the same densities. Activity only shifts the stability region of the solid to higher densities. Additionally, in the bottom part of the figure we mark the densities where the defects disappear for larger systems of N = 72 × 10 3 particles as a cross and the hexatic-solid transition densities obtained from the decay of positional correlations as a star, similar to the ones marked in Fig. 3. We find that in the passive case the discontinuity in K agrees well with the hexatic-solid transition. However, for Pe 7.2 the discontinuity in K agrees well with the (dis)appearance of defects.
In equilibrium systems, the KTHNY theory suggests a critical value of 16π for the renormalized Young's modulus βK R σ 2 , below which the solid is unstable to shear. The renormalization procedure corrects for the interactions of defects at finite temperature. We apply the renormalization procedure to the elastic constants up to Pe 2.4 for which there is a finite fraction of defects at the hexatic-solid transition.
To evaluate the renormalized Young's modulus we first explicitly measure the probability of dislocation pairs p d in simulations and calculate the core energy E c of defects using p d = exp(−2βE c )Z (K ), where Z (K ) is the internal partition function of a dislocation [21,22]. In equilibrium, due to thermal fluctuations there is a finite probability for the formation of dislocation pairs and the dislocation energy E c near the melting transition is small but finite. As we increase the activity the concentration of dislocation defects near the melting transition decreases (see Fig. 3). This observation hints that the energy needed to create a dislocation pair becomes higher as we increase activity. Conversely, we can interpret that the unbinding energy decreases with increasing activity, which eases the dissociation of dislocation pairs into dislocations. We then apply the recursion relations of the KTHNY theory [21,22] to obtain the renormalized Young's modulus K R , shown in Fig. 4(b) as a function of density for Pe = 0, 0.8, 1.6, and 2.4. In the same plot we also show the bare values K as in Fig. 4(a). For Pe = 0, we find that the renormalized K R differs significantly from the bare value. The density at which K R drops to zero agrees well with our estimate of the hexatic-solid transition. However, as we increase the activity up to Pe = 2.4 we find that this is no longer valid. Hence, even for a small Pe we can already see that the predictions of KTHNY theory based on the elastic constants deviate significantly from the transitions as obtained from the positional correlations of the particles. At higher activity there is a complete absence of defects at the hexatic-solid transition point and the transition is entirely driven by activity.

IV. CONCLUSION
We found that at high activity the solid-hexatic transition of 2D active Brownian particles cannot be mediated by topological defects as both the solid and hexatic states are devoid of any defects in contrast to passive systems. We speculate that the solid-hexatic transition may be driven by a growing length scale of regions of cooperative motion. To provide support of a growing length scale, we measure the particle displacement u i ( t ) = r i (t + t ) − r i (t ) for each particle and show in [16] that the regions of cooperative motion increase in size with increasing density. Interestingly, we also observed that the elastic constants of active solids are equal to those of the passive counterparts, as the swim contribution to the stress tensor is negligible in the solid state. The activity only shifts the stability regime of the solid state to higher densities. acknowledge funding from the Industrial Partnership Programme "Computational Sciences for Energy Research" (Grant No. 14CSER020) of the Foundation for Fundamental Research on Matter, which is part of the Netherlands Organization for Scientific Research. This research program is cofinanced by Shell Global Solutions International B.V.