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

TL;DR

This study investigates how topological defects influence the melting process and elastic properties of active Brownian particles, revealing deviations from classical theories and showing that activity shifts the solid regime without altering elastic constants.

Contribution

It demonstrates that active particle melting involves defect dynamics different from KTHNY theory and shows elastic constants are unaffected by activity, only shifting the solid phase.

Findings

Crystalline states melt into hexatic states without dislocation unbinding.

Elastic constants of active solids match those of passive solids.

Activity shifts the stable solid regime to higher densities.

Abstract

We find that crystalline states of repulsive active Brownian particles at high activity melt into a hexatic state but this transition is not driven by an unbinding of bound dislocation pairs as suggested by the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) 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.

Equations

$$\gamma {\dot{\mathbf{r}}}_i=-\sum _{j\ne i}{\nabla }_i\mathrm{U}(r_{ij})+\gamma v_0{\mathbf{e}}_i+\sqrt{2\gamma k_{B}T}{\mathbf{\Lambda }}_i^t,$$

$$\gamma {\dot{\mathbf{r}}}_i$$

$${\dot{\theta }}_i$$

$$⟨{\mathbf{\Lambda }}_i^t(t)⟩=0,$$

$$⟨{\mathbf{\Lambda }}_i^t(t){\mathbf{\Lambda }}_j^t(t^{\prime })⟩$$

$$⟨{\Lambda }_i^r(t){\Lambda }_j^r(t^{\prime })⟩$$

$$\mathrm{U}(\mathrm{r})$$

$$=0$$

$${\psi }_6({\mathbf{r}}_i)=\frac{1}{N_b}\sum _{j\in N_b}\exp (\iota 6{\theta }_{ij}),$$

$$g_T(\text{r})=⟨{\psi }_T^{\ast }({\mathbf{r}}^{\prime }+\mathbf{r}){\psi }_T({\mathbf{r}}^{\prime })⟩,$$

$${\psi }_T({\mathbf{r}}_i)=\exp (\iota {\mathbf{k}}_0\cdot {\mathbf{r}}_i).$$

$$P_{\alpha \beta }=P_{\alpha \beta }^{\text{vir}}+{\delta }_{\alpha \beta }(P_{\alpha \beta }^{\text{id}}+P_{\alpha \beta }^{\text{swim}}).$$

$$P_{\alpha \beta }^{\mathrm{vir}}=-\frac{1}{4V}⟨\sum _i^N\sum _{j\ne i}^N{\partial }_{{\mathbf{r}}_{i,\beta }}\mathrm{U}(r_{ij})\cdot ({\mathbf{r}}_{i,\beta }-{\mathbf{r}}_{j,\beta })⟩.$$

$$P^{\mathrm{swim}}=\frac{\gamma \rho v_0^2}{2D_r}-\frac{\gamma v_0}{4V{D}_r}⟨\sum _{i=1}^N\sum _{j\ne i}^N{\nabla }_i\mathrm{U}(r_{ij})\cdot {\mathbf{e}}_i⟩.$$

$$\mathbb{C}=\left[\begin{array}{cccc}{C}_{11}& C_{12}& 0& 0\\ & C_{22}& 0& 0\\ & & 0& 0\\ & & & C_{44}\end{array}\right]=\left[\begin{array}{cccc}\lambda +2\mu & \lambda & 0& 0\\ & \lambda +2\mu & 0& 0\\ & & 0& 0\\ & & & \mu \end{array}\right]$$

$$B_{\alpha \beta \gamma \delta }=C_{\alpha \beta \gamma \delta }-P({\delta }_{\alpha \gamma }{\delta }_{\beta \delta }+{\delta }_{\alpha \delta }{\delta }_{\beta \gamma }-{\delta }_{\alpha \beta }{\delta }_{\gamma \delta })$$

$$B_{11}=C_{11}-P,$$

$$B_{12}=C_{12}+P,$$

$$E=\lambda +\mu ,G=\mu ,K=\frac{4\mu (\lambda +\mu )}{\lambda +2\mu }=\frac{4EG}{E+G}.$$

$$\frac{1}{E}=\kappa =-\frac{1}{V}{\left(\frac{\partial V}{\partial P}\right)}_T=\frac{1}{\rho }{\left(\frac{\partial \rho }{\partial P}\right)}_T,$$

$$B_{11}=\frac{\partial P_{xx}}{\partial {ϵ}_{xx}},B_{22}=\frac{\partial P_{yy}}{\partial {ϵ}_{yy}},B_{12}=\frac{\partial P_{yy}}{\partial {ϵ}_{xx}},B_{44}=\frac{\partial P_{xy}}{\partial {ϵ}_{xy}},$$

$$p_d$$

$$=\exp \left(-\frac{2E_c}{k_{B}T}\right)\frac{2\pi \sqrt{3}}{\beta K/8\pi -1}{I}_0\left(\frac{\beta K}{8\pi }\right)\exp \left(\frac{\beta K}{8\pi }\right)$$

$$\frac{\partial }{\partial l}\left(\frac{8\pi }{\beta K(l)}\right)$$

$$\frac{\partial y(l)}{\partial l}$$

$$y=\exp \left(-\frac{E_c}{k_{B}T}\right).$$

Full text

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-ranged bond orientational order, and melts subsequently via a second continuous transition into a fluid state Kosterlitz and Thouless (1973); Nelson and Halperin (1979); Chui (1982). These transitions are triggered by the unbinding of topological defects, which are particles with a non-conforming 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 5-fold and 7-fold 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-ranged repulsive disks was only settled a few years ago, both in simulations Bernard and Krauth (2011); Engel et al. (2013); Qi et al. (2014); Kapfer and Krauth (2015); Anderson et al. (2017) and experiments Thorneywork et al. (2017), showing a two-stage melting scenario, that deviates from the KTHNY scenario, of a continuous solid-hexatic transition driven by an unbinding of dislocation pairs and a first-order hexatic-fluid transition due to a proliferation of grain boundaries.

Remarkably, non-equilibrium systems of self-propelled particles Fily and Marchetti (2012); Redner et al. (2013); Stenhammar et al. (2014); Siebert et al. (2018), which constantly convert energy from the environment into persistent motion, have recently also been shown to follow such a two-stage melting behaviour Digregorio et al. (2018); Klamser et al. (2018). It was found that the first-order nature of the liquid-hexatic transition persists upto a small degree of activity. The transition then becomes continuous until it reappears as a coexistence of dilute and dense states at high activity Digregorio et al. (2018); Klamser et al. (2018). 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:

[TABLE]

where $v_{0}\mathbf{e}_{i}$ is the self-propulsion speed, $\gamma$ is the damping coefficient, $k_{B}$ the Boltzmann constant, and $T$ the temperature of the solvent. The particle orientation $\mathbf{e}_{i}=(\cos\theta,\sin\theta)$ undergoes free rotational diffusion $\dot{\theta}_{i}=\sqrt{2D_{r}}\Lambda_{i}^{r}$, where $D_{r}$ is the rotational diffusion coefficient. The quantities $\boldsymbol{\Lambda}_{i}^{t}$ and $\Lambda_{i}^{r}$ are unit-variance Gaussian noise terms with zero mean. The particles interact with a pairwise repulsive WCA potential $\mathrm{U}(r)=4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]+\varepsilon$, and we set $k_{B}T/\varepsilon=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-$\rho$) representation, where the Péclet number is defined by $\mathrm{Pe}=v_{0}\tau/\sigma$. 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 $\eta_{T}\leq 1/3$ as described by the equilibrium KTHNY theory.

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 5-fold and 7-fold defects. In Fig. 2 we show typical configurations, highlighting these defects in the fluid, hexatic and crystalline state at $\mathrm{Pe}=0,4.8$ and $71.5$. At low activity $\mathrm{Pe}=4.8$, similar defect configurations are observed as for the passive hard-disk systems Qi et al. (2014). 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 of many defect clusters and 5-fold and 7-fold disclinations, both destroying the local bond orientational order. However, the hexatic states at $\mathrm{Pe}=71.5$, are significantly different in defect configurations. The low-density hexatic state at $\rho\sigma^{2}=1.420$ shows a high number fraction of dislocations with almost no dislocation pairs, whereas the high-density hexatic state ($\rho\sigma^{2}=1.600$) at the same Pe shows a complete absence of any defects.

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_{\mathrm{total}}/N$ increases upon decreasing the density for all Pe. At low $\mathrm{Pe}\leq 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_{\mathrm{quart}}/N$ increases slightly upon approaching the hexatic-solid transition, whereas the number fraction of dislocations $N_{\mathrm{pair}}/N$ increases more rapidly. At the hexatic-solid transition, $N_{\mathrm{quart}}/N$ is still higher than $N_{\mathrm{pair}}/N$ but at lower density this scenario is reversed. The high fraction of $N_{\mathrm{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_{\mathrm{free}}/N$ of 5- and 7-fold 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 2D passive systems of short-range repulsive particles, persists at low activity.

For higher activity $\mathrm{Pe}\geq 7.2$, the behavior of dislocations is different than for passive systems as shown in Fig. 2. For $\mathrm{Pe}\geq 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_{\mathrm{quart}}$ never exceeds $N_{\mathrm{pair}}$. This observation is more clearly visible in Fig. 3(b) where we plot the differences $(N_{\mathrm{pair}}-N_{\mathrm{quart}})/N$ and $(N_{\mathrm{pair}}-N_{\mathrm{free}})/N$. We find that for $\mathrm{Pe}\leq 2.4$ and for high densities $(N_{\mathrm{pair}}-N_{\mathrm{quart}})/N$ is negative whereas for $\mathrm{Pe}\geq 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 Fig. 3(a) and (b) we also mark the hexatic-solid transition densities as identified from the spatial decay of positional correlations SI . These densities agree closely with the minimum in $(N_{\mathrm{pair}}-N_{\mathrm{quart}})/N$ which corresponds to a reversal in the trend of defect concentrations and confirms our previous observation that for $\mathrm{Pe}\leq 2.4$ the hexatic-solid transition is driven by the unbinding of dislocation pairs. The difference $(N_{\mathrm{pair}}-N_{\mathrm{free}})/N$, however, is always positive for all densities and show 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 $\mathrm{Pe}$, marked by a cross on the $x$-axis and the corresponding density range by cross-hatching. As we increase the activity, this point crosses over from the dense crystal state to the hexatic state. For $\mathrm{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 kinds 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 non-equilibrium 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.

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 behaviour 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\times 10^{3}$ particles and measure the full stress tensor $P_{\alpha\beta}$, comprising of the ideal, virial and swim components Winkler et al. (2015), in the deformed box due to a fixed small linear strain $\epsilon_{xx}\in[-0.01,0.01]$. We calculate the Lamé elastic coefficients $\lambda$ and $\mu$ from the effective stiffness tensor $\mathbb{B}$ obtained from the slope of a linear fit to the stress versus strain curves Landau et al. (1989); Ray (1988); Frenkel and Smit (2001) (see Ref. SI for details).

In Fig. 4(a) we plot the Young’s modulus $K=4\mu(\lambda+\mu)/(\lambda+2\mu)$ obtained for $0.0\leq\mathrm{Pe}\leq 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 $\mathrm{Pe}$. This remarkable result can be explained from the fact that the swim contribution to the stress tensor is negligible in the solid, and hence the elastic moduli for active solids equals the ones of their passive counterparts. 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\times 10^{3}$ particles as a cross ($\times$) as well as the hexatic-solid transition densities obtained from the decay of positional correlations as a star ($\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 $\mathrm{Pe}\geq 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\pi$ for the renormalized Young’s modulus $\beta K_{R}\sigma^{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 $\mathrm{Pe}\leq 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\beta E_{c})Z(K)$ where $Z(K)$ is the ‘internal partition function’ of a dislocation Fisher et al. (1979); Sengupta et al. (2000). 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 reduces (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 reduces with increasing activity which eases the dissociation of dislocation pairs into dislocations. We then apply the recursion relations of the KTHNY theory Fisher et al. (1979); Sengupta et al. (2000) 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 $\mathrm{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 upto $\mathrm{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.

Conclusions: We found that at high activity the 2D melting of active Brownian particle solids into a hexatic state is not driven by an unbinding of dislocation pairs in contrast to passive systems. The hexatic state at high densities is completely devoid of any defects, and this defect-free region widens with activity. The solid-hexatic transition might be driven by a growing length scale of regions of cooperative motion, but this requires further investigation. Interestingly, we 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.

S1 Model

We consider a two-dimensional system of isotropic Brownian particles that exhibit a self-propulsion speed $v_{0}$ which is directed along the orientation vector $\mathbf{e}_{i}=(\cos\theta_{i},\sin\theta_{i})$ assigned to particle $i$. To describe the translational and rotational motion of the individual colloidal particle $i=1,\dots,N$ we employ the overdamped Langevin dynamics:

[TABLE]

where $\gamma$ is the damping coefficient due to the drag forces from the implicit solvent, $k_{B}$ is the Boltzmann constant, and $T$ is the bath temperature. $D_{r}$ is the rotational diffusion coefficient. The quantities $\boldsymbol{\Lambda}_{i}^{t}$ and $\Lambda_{i}^{r}$ are unit-variance Gaussian noise terms with zero mean:

[TABLE]

where $\mathbb{I}_{2}$ is the $2\times 2$ identity matrix. The angular brackets $\langle\cdots\rangle$ denote an average over different realizations of the noise. The particles interact with a short-range repulsive Weeks-Chandler-Andersen (WCA) potential given by:

[TABLE]

where $\mathrm{r}=|\mathbf{r}_{ij}|$ is the distance between the centers of particle $i$ and $j$, $\sigma$ is the particle diameter pertaining to the length scale in WCA potential, and $\varepsilon$ is the strength of the particle interactions.

We set the system temperature $k_{B}T/\varepsilon=1$, following Ref. Redner et al., 2013, and the damping coefficient $\gamma\sigma^{2}/\varepsilon=1$ fixing the translational diffusion coefficient to correspond to the free diffusion of particles given by the Stokes-Einstein relation $D_{t}=\gamma^{-1}k_{B}T$. This sets our time scale as $\tau=\gamma\sigma^{2}/k_{B}T$. The rotational diffusion coefficient is set to $D_{r}\tau=3$ and we use a time step size $dt=10^{-5}\tau$ for numerically integrating the equations of motion. We define a non-dimensional Péclet number $\mathrm{Pe}=v_{0}\tau/\sigma$ as the ratio of the persistence length of motion to the particle diameter and perform simulations using the HOOMD-Blue Anderson et al. (2008); Glaser et al. (2015) package in the range $0\leq\text{Pe}\leq 150$. We used $N=72\times 10^{3}$ particles in an approximately square 2D periodic simulation box with dimensions $L_{x},L_{y}\approx 250\sigma$ for identifying the fluid-hexatic-solid transitions and $N=2.8\times 10^{3}$ for calculating the elastic moduli. We used a regular hexagonally packed arrangement of particles at the overall system density as our initial configuration and we collect snapshots for $200\tau-500\tau$ at an interval of $\tau$ for analysis after allowing the system to achieve a stationary state for about $200\tau$.

As mentioned in the main text, the state diagram shown in Fig. 1 was obtained by measuring the equation of state (pressure-density curves), the density histograms and the decay of the orientational and the positional correlation functions to locate the boundaries as precisely as possible. Specifically, to locate the coexistence, we identified negative slope regions in the pressure-density curves as well as double-peaked structure of the density histograms similar to the analysis presented in Ref. Digregorio et al., 2018 and Ref. Cugliandolo et al., 2017. We describe below the orientational and positional order parameters used for identifying the hexatic and solid phases.

S1.1 Orientational and Positional order

We measure the local 6-fold orientational symmetry around particle $i$ using the hexatic order parameter $\psi_{6}({\mathbf{r}}_{i})$ given by:

[TABLE]

where $N_{b}$ denotes the number of nearest neighbors of particle $i$ and the bond angle $\theta_{ij}$ is measured as a deviation of the orientation of the vector $\mathbf{r}_{ij}$ from the reference global system orientation measured from $\Psi_{6}(L)$ averaged over all the particles. We identify the nearest neighbors $N_{b}$ of the particle by using a Voronoi construction.

To investigate the decay of positional order, we measure the positional correlation function

[TABLE]

where $\psi_{T}({\mathbf{r}}_{i})$ is the positional order parameter expressed as:

[TABLE]

Here $\mathbf{k}_{0}$ is the vector in reciprocal space denoting one of the first Bragg peaks in the 2D structure factor $S(\mathbf{k})$. The magnitude of this vector is equal to that of the reciprocal lattice vector i.e. $\mathbf{k}_{0}=(0,4\pi/a\sqrt{3})$ where $a=(2/\rho\sqrt{3})^{1/2}\sigma$ is the lattice spacing in a regular hexagonal packing at a number density $\rho$ in a 2D geometry. According to the KTHNY theory, the positional order of a two-dimensional solid decays algebraically as $g_{T}(r)\propto r^{-\eta_{T}}$ with an exponent $0\leq\eta_{T}\leq 1/3$. Upon melting, the decay of the positional correlations becomes exponential i.e. $g_{T}(r)\propto\exp(-r/\xi_{T})$ with a correlation length $\xi_{T}$, which decreases with decreasing density. We show the positional correlation functions $g_{T}(r)$ as a function of particle separation $r$ in Fig. S1 and extract the correlation lengths $\xi_{T}$ in the case of an exponential decay or the exponent $\eta_{T}$ in the case of an algebraic decay. We identify the hexatic-solid transition by locating the density at which the exponent $\eta_{T}$ becomes smaller than $1/3$. For $\mathrm{Pe}=0$, we find that the decay of $g_{T}(r)$ becomes algebraic with $\eta_{T}\approx 1/3$ at $\rho\sigma^{2}=0.926$, marking the hexatic-solid phase transition. We locate the transition densities for higher Pe in a similar manner.

S1.1.1 System-size dependence

To check the finite-size effects on the decay of the positional order we also simulate a few cases with $N=288\times 10^{3}$ particles. We show $g_{T}(r)$ for the two different system sizes at $\mathrm{Pe}=2.4,7.2$ and $71.5$ and with densities near the respective hexatic-solid transition in Fig. S2. The values for the exponents $\eta_{T,1}$ and $\eta_{T,2}$ of the power-law decay, for the small and large systems, respectively, are also listed in Fig. S2. For smaller system sizes, we observe clearly that $\eta_{T,1}<1/3$ for all three $\mathrm{Pe}$ values, which corresponds to the solid phase. For $N=288\times 10^{3}$, we find that the exponent $\eta_{T,2}$ is still close to $\eta_{T,1}$ for $\mathrm{Pe}=2.4$, but for higher Pe the exponent $\eta_{T,2}$ is significantly larger than $\eta_{T,1}$. Despite these differences in the decay of the positional correlations, the location of the density where the defects disappear does not change and the observed hexatic region devoid of defects is robust over these investigated system sizes.

S2 Elastic moduli

S2.1 Stress tensor

In order to measure the bulk pressure $P$ in our system consisting of $N$ active Brownian particles, we employ the expressions as introduced by Winkler et al.Winkler et al. (2015), but modified them to the 2D case. Specifically, the pressure for isotropic active particles in a periodic box with lateral dimensions $L_{x}$ and $L_{y}$ and 2D ‘volume’ $V=L_{x}L_{y}$ is calculated using $P=\mathrm{Tr}(\mathbb{P})$ where the full stress tensor $\mathbb{P}$ is given by:

[TABLE]

Here $P^{\mathrm{id}}$ is the ideal gas pressure given by $P^{\mathrm{id}}=\rho k_{B}T$ with $\rho=N/V$ the number density of the particles. The virial contribution $P^{\mathrm{vir}}$ is obtained using the standard virial expression

[TABLE]

The swim pressure contribution $P_{\mathrm{swim}}$ due to the self-propulsion is given by:

[TABLE]

S2.2 Stiffness tensor and Lamé elastic coefficients

In the linear elastic theory of isotropic solids, the elastic moduli relate the stress response of a system to an applied strain. In equilibrium, the elastic moduli are related to the free energy change due to such deformations Landau et al. (1989). Instead, for non-equilibrium systems we directly assume Hooke’s law which linearly relates the mechanical stress $P_{\alpha\beta}$ with the applied strain $\epsilon_{\gamma\delta}$ through a symmetric stiffness tensor $\mathbb{C}$ given by:

[TABLE]

where $\lambda$ and $\mu$ are the Lamé coefficients in equilibrium systems Landau et al. (1989). For conciseness, we follow the Voigt notation above for indexing $C_{\alpha\beta\gamma\delta}$ with $xx$=1, $yy$=2, and $xy$=4. If a system is under a uniform isotropic pressure, the stiffness tensor $\mathbb{C}$ can be rewritten in terms of an effective stiffness tensor $\mathbb{B}$ as Ray (1988); Frenkel and Smit (2001):

[TABLE]

[TABLE]

where $P=(P_{xx}+P_{yy})/2$ is the uniform pressure. The bulk modulus $E$, the shear modulus $G$ and the Young’s modulus $K$ are related to the Lamé coefficients in 2D as:

[TABLE]

Furthermore, from equilibrium statistical thermodynamics the isothermal compressibility $\kappa=1/E$, where $E$ is the bulk modulus, is expressed as:

[TABLE]

which can also be measured directly from the slope of $P-\rho$ curves. Once the uniform pressure of the system and the stiffness tensor (or the effective stiffness tensor $\mathbb{B}$) are known, we obtain the elastic moduli from $\lambda$ and $\mu$ using Eq. S11. In our simulations we apply the method of box deformations to numerically evaluate the stiffness tensor for a system of interacting particles in an $NVT$ ensemble Frenkel and Smit (2001); Clavier et al. (2017).

We extract the four non-zero elements of the stiffness tensor $\mathbb{C}$ by performing three kinds of deformations of the simulation box following Ref. Clavier et al., 2017. In the first kind of deformation, the box is elongated or compressed along the $x$-direction by a small factor $\epsilon_{xx}$ such that the particle coordinates in the $x$-direction become $x^{\prime}=x(1+\epsilon_{xx})$ and the box length also becomes $L_{x}^{\prime}=L_{x}(1+\epsilon_{xx})$. Similarly, the box can be elongated or compressed along the $y$-direction corresponding to imposing a small linear strain $\epsilon_{yy}$. Both these deformations correspond to a change in the overall density of the system but the magnitude is kept small in order to stay in the linear response regime. The third deformation is of a shearing type in which we change the shape of the box by keeping the volume constant. The angle between the $x$ and $y$ dimension box vectors, $\mathbf{a}$ and $\mathbf{b}$ respectively, is changed from $\pi/2$ to $\pi/2-\tan^{-1}(\epsilon_{xy})$. The particle positions are then transformed as $(x,y)\rightarrow(x+y\epsilon_{xy},y)$.

In our simulations, we start from a perfect hexagonal initial configuration with $N=2.8\times 10^{3}$ particles and deform the box corresponding to the applied strain. We then measure the full stress tensor $P_{\alpha\beta}$ after a sufficiently long equilibration time that allows the system to reach a steady state. We perform the measurements by applying fixed linear strain $\epsilon_{xx}\in[-0.01,0.01]$ in intervals of 0.004. For an isotropic solid only the first two elements $C_{11}$ and $C_{12}$ are sufficient to obtain the Lamé coefficients $\lambda$ and $\mu$, which can be measured just by applying a longitudinal strain $\epsilon_{xx}$. However, for some cases we also measure the values of $\mu$ obtained by imposing a shearing strain $\epsilon_{xy}$ and confirm that the two independent measurements agree. The effective stiffness tensor $\mathbb{B}$ is directly obtained from the slope of a linear fit to the stress vs. strain curves, as shown for $B_{11},B_{12}$ and $B_{44}$ in Fig. S4, using

[TABLE]

S2.3 Bulk and Shear elastic moduli

In Fig. S5(a) and S5(b) we plot the bulk modulus $E$ and the shear modulus $G$, respectively, as a function of density for various Pe obtained using the method described above. For $\mathrm{Pe}=0$ (magnified in the inset) we find that there is a distinct jump in both $E$ and $G$, as indicated in the figure by a blue arrow, at a density of $\rho\sigma^{2}=0.926$. This jump is indicative of the second order nature of the transition. Upon increasing Pe, we observe a similar jump appearing in both $E$ and $G$ at higher densities marked by arrows in the figure. The bulk modulus $E$ shows only a discontinuity for higher Pe but the shear modulus $G$ shows a sharp drop to very small values at this transition upon reducing the density. Such a small value of the shear modulus $G$ indicates that the system is not a solid anymore and undergoes plastic deformation upon shearing. Furthermore, in the same plots we also indicate the densities where we observe a finite number of defects in the simulations with $N=2.8\times 10^{3}$ particles by a plus marker ($+$) as in the main text Fig. 4(a). These points were determined by analyzing the sampled snapshots within our simulated time which show a complete absence of defects at densities higher than the marked points ($+$).

For $\mathrm{Pe}=0$ the defects disappear at a density of $\rho\sigma^{2}=0.950$ which is much higher than the point $\rho\sigma^{2}=0.926$ at which we observe the jump in the elastic moduli. For $\mathrm{Pe}\geq 7.2$ we find that the two transition points agree extremely well. This indicates that for active cases the system becomes plastic as soon as a finite number of defects, mainly dislocations, appear in the system. On the other hand, the elastic moduli of the active solid states as a function of density collapse onto a single master curve independent of Pe. This remarkable result can be explained by the fact that the swim contribution to the stress tensor is zero or negligible in the solid phase, and hence, the elastic constants of active solids become equal to those of passive solids at the same density. The sole effect of activity is that the stable solid regime shifts to higher densities with activity. A numerical fit of the form $E,G\propto\exp(a\rho^{3}+b\rho^{2}+c\rho+d)$ is shown as a black solid line in both Fig. S5(a) and Fig. S5(b), and agrees very well with the measurements.

S3 Renormalization procedure from KTHNY theory

In equilibrium systems, the KTHNY theory suggests that the melting transition is accompanied by a lowering of the Young’s modulus $\beta K$ below a critical value of $16\pi$. To correct for the interactions of defects present at a finite temperature a renormalization group analysis is applied to obtain the renormalized value $K_{R}$ of the Young’s modulus which can then be compared against the numerical value of $16\pi$ to identify the melting transition. The theory describes the dislocation defects in 2D systems associated with a ‘core energy’ $E_{c}$ Kosterlitz and Thouless (1973); Nelson and Halperin (1979). The probability of finding a bound pair of such dislocation defects is given by Fisher et al. (1979); Sengupta et al. (2000):

[TABLE]

where $Z(K)$ is the internal partition function of a dislocation, and $I_{0}$ is a modified Bessel function. The theory suggests a continuous transition from the solid to the hexatic state for large core energies $E_{c}\geq 2.8k_{B}T$ and a weakly to strongly first-order transition as $E_{c}$ approaches and becomes lower than a value of $2.8k_{B}T$ Strandburg (1988). Typically, for systems with hard-core interactions the value of $E_{c}$ near the solid-hexatic transition is $\sim 6k_{B}T$ as found in Ref. Qi et al. (2014) for monolayers of hard spheres.

The renormalization group recursion relations for the Young’s modulus $K$ are expressed as Kosterlitz and Thouless (1973); Nelson and Halperin (1979); Sengupta et al. (2000):

[TABLE]

where the fugacity $y$ of the dislocation-pair fluid is obtained from an estimate of the core energy $E_{c}$ as:

[TABLE]

The differential equations Eq. S14-S15 can be solved recursively for $l=0\dots\infty$ by using the unrenormalized (‘bare’) values $K(0)=K$ and $y(0)=\exp(-E_{c}(K(0)))$ as the initial guesses for $l=0$ and utilizing a trapezoidal (or higher order scheme) for performing the integration. The renormalized values are obtained from the renormalization-flow diagram of $y$-vs-$1/K$ (Fig. 1 in Ref. Sengupta et al. (2000)) for the separatrix and $K_{R}=K(\infty)$ when $y(\infty)=0$. Exactly at the transition, the renormalization-flow follows the separatrix and above ($T>T_{m}$) and below ($T<T_{m}$) the melting point goes to the end points $\infty$ and [math], respectively.