Crystallisation and polymorph selection in active Brownian particles

We explore crystallisation and polymorph selection in active Brownian particles with numerical simulation. In agreement with previous work (Wysocki et al. in Europhys Lett 105:48004, 2014), we find that crystallisation is suppressed by activity and occurs at higher densities with increasing Péclet number (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${ Pe }$$\end{document}Pe). While the nucleation rate decreases with increasing activity, the crystal growth rate increases due to the accelerated dynamics in the melt. As a result of this competition, we observe the transition from a nucleation and growth regime at high \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${ Pe }$$\end{document}Pe to “spinodal nucleation” at low \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${ Pe }$$\end{document}Pe. Unlike the case of passive hard spheres, where preference for FCC over HCP polymorphs is weak, activity causes the annealing of HCP stacking faults, thus strongly favouring the FCC symmetry at high \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${ Pe }$$\end{document}Pe. When freezing occurs more slowly, in the nucleation and growth regime, this tendency is much reduced and we see a trend towards the passive case of little preference for either polymorph.


Introduction
The field of Active Matter may be said to consider systems of organisms or artificial bodies that consume energy for self-propulsion [1]. On mesoscopic length scales (nm to µm), it is concerned with describing the dynamics of biological microswimmers [2] such as bacteria and motile cells [3]. The dynamics of active matter in unbounded, homogeneous, and low-Reynolds number environments are well described by active Brownian motion and run and tumble dynamics [4]. Observation of matter behaving according to this description has led to the discovery of unique dynamical phenomena such as motility-induced phase separation (MIPS) [5], where bodies packed at densities greater than a critical volume fraction and with sufficient propulsion will separate into a dense phase and a dilute phase in the absence of attraction [6].
Key to the development of a better theoretical understanding is to use simple models of active particles which capture some of the complex behaviour observed experimentally, for example collective motion and demixing [5,[7][8][9][10][11][12][13][14][15][16]. In this context, simple model systems, such as active colloids, play an important role a e-mail: fergus.moore@bristol.ac.uk (corresponding author) b e-mail: paddy.royall@espci.fr c e-mail: john.russo@uniroma1.it and these may be modelled with the use of active Brownian particles (ABPs) [6]. Besides being the source of novel phenomena, activity can also fundamentally alter the nature of behaviour already observed in passive systems, such as crystallisation. Although certain aspects of crystallisation in passive colloidal systems, such as the nucleation rate at low supersaturation, are still poorly understood [17][18][19][20][21][22][23][24][25][26][27], at higher supersaturation, where crystallisation occurs on the timescales accessible to brute force computer simulations, very good agreement is found between experiment and simulation [28]. At higher colloidal volume fraction still, the barrier to nucleation falls so much that rather than conventional nucleation-and-growth, the system undergoes "spinodal nucleation", where, relative to the intrinsic structural relaxation time τ α , the timescale for crystallisation falls dramatically such that it is well below the relaxation time [28][29][30][31][32].
Another property of crystallising systems is polymorphism, i.e. the ability of a material to nucleate different crystalline phases, and whose understanding is fundamental to predict the structure of the growing nuclei. So far, our understanding of polymorphism is based on equilibrium thermodynamic principles, such as the Ostwald step rule of phases [33], stating that the first solid formed is not the thermodynamically most stable, but the state nearest in free energy to the original state. For hard spheres, which may be said to constitute the passive equivalent of the ABP system that we con-sider, two different crystalline polytypes are observed during nucleation: either face-centered-cubic (FCC) or hexagonal-close-packed (HCP), and the difference in all thermodynamic relevant quantities (such as free-energy, nucleation barrier, and stacking-free energy) between the competing polymorph are negligibly small (within 10 −3 k B T per particle for all cases) [34,35]. Thermodynamics thus dictate that the early stages of nucleation should produce an almost equal amount of FCC and HCP for hard-spheres. We will show that activity can significantly alter this result. Studying the effect of activity is thus an important step towards understanding polymorphism in out-of-equilibrium situations.
The effect of activity upon crystallisation has been studied in the context of the effect on the state diagram [36]. In both two [37,38] and three dimensions [12,[39][40][41], the freezing line is found to move to higher area or volume fraction as a function of activity. In this sense, activity may be said to suppress crystallisation. The effect of activity on the process of nucleation has been studied via classical nucleation theory, in which a renormalised surface tension was found to provide reasonable agreement with simulation [42]. At higher activity in dimension d = 3, the active fluid that coexists with a low-density active fluid through MIPS has a very high volume fraction [40,41] and crystal nucleation requires rare fluctuations that exhibit the nearly closepacked volume fraction of the solid [40]. One intriguing and unexpected effect of activity upon crystallisation was the observation of annealing of grain boundaries in the case of the addition of a small quantity of active particles to an otherwise passive system [43].
To date, there have been relatively few experiments with active colloids at high density where crystallisation due to excluded volume interactions is seen [44]. This is due in no small part to the difficulties in stabilising active colloids at high density against aggregation. However, recently this has begun to change and excluded volume interactions have driven ordering in a few experiments in two dimensions [45][46][47][48]. The study of 3d active colloids is in its infancy; however, one system that has emerged of active multi-polar colloids [49] does exhibit crystallisation to a variety of polymorphs also exhibited by related passive dipolar colloids [50,51].
In this work, we consider crystallisation regimes in a system of active Brownian particles in three dimensions. In particular, we investigate analogous behaviour to the nucleation-and-growth and spinodal regimes observed in passive colloidal systems. Furthermore, we find an unexpected polymorph selection phenomenon that is uniquely distinct from those observed in passive systems.
This article is organised as follows. In Sect. 2, we describe the methodology used for the simulation runs and the analysis of topological clusters in the fluid. Results are presented in Sect. 3, with subsections dedicated to the state diagram (Sect. 3.1), the dynamical and structural properties of the active fluid (Sect. 3.2), and the nucleation and crystal growth behaviour (Sect. 3.3). We summarise our findings in Sect. 4.

Computer simulations
We model active colloids as active Brownian particles, which propel with a constant velocity V 0 , along their individual direction vectors e, which in turn are subject to rotational diffusion. We implement this model through molecular dynamics simulations using a customised version of the open source LAMMPS package [52], which integrates the following equations of motion:ṙ Here,ṙ is the particle velocity, V 0 is the magnitude of the constant applied active velocity, and F is the inter-particle force. The thermal fluctuations promoting translational diffusion are included in the Gaussian white-noise term η, where η = 0, and D t is the translational diffusion coefficient. Thermal noise driving rotational diffusion of the direction vector e is represented by ξ, where ξ = 0, and D r is rotational diffusion coefficient. The two diffusion coefficients are related via D t = D r σ 2 /3. For all simulations in this work β = 5, m = 1, σ = 1. Our measure of time is the characteristic rotational diffusion time τ r = 1/(2D r ) [39].
The active particles are modelled as being similar to hard spheres and to achieve this we include a Weeks-Chandler-Andersen (WCA) inter-particle potential in the force term in equation (1), which takes the form: where ε is the interaction energy, r is the inter-particle distance.
Since we use the WCA interaction, we cannot assume the hard particle diameter σ to define a volume fraction. Furthermore, methods that determine an effective particle diameter such as Barker-Henderson effective hard sphere diameter [53] may not hold outside of equilibrium systems. Therefore, as in ref. [54], we use the total density ρ = N/V , where N is the number of particles and V is the volume of the system.
We use the Péclet number to refer to the relative strength of the activity in the system, which we define as: Pe = V 0 /σD r . Throughout this work, we keep D r constant at D r = 1x and vary Pe by changing the propulsion velocity V 0 . Previous studies [12,54] have shown that when the propulsion force from the activity increases relative to the repulsive force from the WCA interaction, the particles become softer. This change was observed to manifest in a shifting of the MIPS phase boundary to higher Pe and ρ, an effect that can be reduced by using an inter-particle potential that more closely resembles that of the hard sphere interaction [54]. However, as we focus only on crystallisation in this work, we perform simulations at low Pe and well below the MIPS phase boundary, and thus, this effect in our WCA system is negligible.
To prepare the high density initial configurations for the simulations, we use the Lubachevsky-Stillinger algorithm [55]. This takes a given box size and number of particles and slowly grows and displaces the particles from σ = 0.1 until they reach σ = 1 with minimal overlaps. To guard against the presence of any small but sufficient remaining particle overlaps, we perform a pre-run simulation with a soft potential: where the constant A is ramped from 0 to 100 over and r c = 2 1 6 . This is run for 1.2τ R without activity. Following this, we perform our data collection runs with particles following the equations of motion outlined in (1) and (2), and the Weeks-Chandler-Anderson (WCA) inter-particle potential (3).
We perform all simulations in a periodic cubic box of dimension length L = 27.5σ and vary the number of particles from 18,000 to 24,000 to explore a range of densities. This system size is such that the largest critical nucleus observed in this work was comprised of less than 4% of the particles in the system, to avoid the finite size effects that have been studied for seeded nucleation in the NVT ensemble [56]. For the determination of structure in longtime steady states, we run for 7200τ R and average over 10 independent configurations. For analysis of nucleation dynamics, we run for 600τ r and average over 20 independent configurations.

Dynamical analysis
The structural relaxation time τ α provides a useful metric through which we can understand the effects of active systems on the verge of crystallisation. We compute τ α for various φ and Pe, through calculation of the self-part of intermediate scattering function: where k is the wavevector k = |k|, taken as 2π/σ. We define τ α as F s (k, τ α ) = e −1 .

Topological cluster classification analysis
Local structures identified in this work are identified by the Topological Cluster Classification (TCC) [57].
The TCC algorithm analyses structure through clusters. To identify a cluster, the TCC uses a modified Voronoi construction to identify a bond network with a cutoff r c = 1.8σ and a four-membered ring parameter f c = 0.82. We identify clusters through calculation of the shortest path 3, 4, and 5 membered rings in the bond network. For non-crystalline clusters, we consider only the minimum energy clusters of the Lennard-Jones interaction, specifically: 5A, 6A, 7A, 8B, 9B, 10B, 11C, 12B, and 13A [58]. Here, the numbers denote the number of particles in each cluster and the lettering signifies the cluster geometry [59]. Furthermore, we use the TCC to identify crystal structure, where 13 particle FCC or HCP clusters are determined through a central particle and its 12 nearest neighbours. We quantify the degree to which a particular structure appears in a configuration as the cluster population N c /N , where N c is the number of particles in a given cluster, and N is the number of particles in the system. It is important to note that a particle can belong to more than once cluster. For example, under certain conditions a particle can belong to both an FCC cluster and an HCP cluster, and when comparing such cluster populations, N c /N will not sum to 1.

Results
In this section, we study the dynamical features of ABP at high density, where an ordered crystalline phase is found to spontaneously form in simulations. We will trace the boundaries of the crystal region, distinguishing between state points that nucleate through a nucleation and growth mechanism, and those that display spinodal nucleation. Through the Topological Cluster Classification (TCC) method, we will distinguish between the FCC and HCP structures and consider the effects of activity on polymorph selection.

State diagram
In Fig. 1, we show the state points we consider and the results of our simulations. While other work has addressed the phase diagram of active Brownian particles in two and three dimensions with respect to MIPS [60,61], or reported the full phase diagram of active disks, in which the freezing line is affected by activity [38]; here, we distinguish the crystallisation regimes of nucleation and growth and spinodal by inspection of the crystallinity as a function of time data (see Sect. 3.3).
In particular, we identify behaviour compatible with the passive WCA system for Pe = 0. Recall that we carry out brute force simulations of N ≥ 18,000 and for a run time of 7200τ R . Therefore, we do not obtain the equilibrium phase diagram [23,24]. Rather, for the passive case, we find nucleation and growth at density ρ = 0.89 and spinodal crystallisation at ρ = 0.91. Moreover, since in the passive case we observe an equilibrium system, we can convert these densities to effective volume fractions φ, via the Barker-Henderson method [53]. This comes out as φ = 0.56 for nucleation and growth and φ = 0.57 for spinodal crystallisation, which is consistent with previous work [28][29][30][31][32], noting that for numerical work such as this system size and runtime have significant consequences. We find similar trends to previous work which considers the effect of activity in two [37,38] and three dimensions [12,[39][40][41] in which the freezing line moves to higher volume (or area) fraction as a function of activity. We also note that the boundary between nucleation-and-growth and spinodal crystallization depends on system size, as large system sizes have a lower nucleation time shifting the transition between the two regimes to lower volume fractions (or generally to lower supercoolings [62]).

Dynamical and structural response of the active WCA fluid to activity
The intrinsic dynamics play an important role in setting the timescale of crystallisation. In this context, the dynamical response of supercooled liquids to activity has been found to be highly complex and to exhibit qualitatively different responses to activity, from accelerating to slowing down and even nonmonotonic behaviour [63][64][65][66]. At the densities, we consider (ρ = 0.72 to ρ = 1.15), in Fig. 2, we see that upon increasing activity, the system accelerates and the structural relaxation time drops. Note that since we consider a monodisperse system, strong supercooling is not possible as crystallisation intervenes. In Fig. 3a, we show the two-body structure of the active fluids via the radial distribution function g(r). The relationship between two-body structure and dynamics has been analysed in some detail in active systems [63][64][65][66] and we see a familiar trend here, of a weakening of the strength of correlations as the relaxation time falls, which here is driven by an increase in activity (Fig. 2). This is particularly evident in the first minima and subsequent maxima and minima in the inset of Fig. 3a.
We also consider the response of higher-order structure to activity in Fig. 3b. Previously this has been found to develop as an increase in the population of locally favoured structures in the Wahnström binary Lennard-Jones model, with activity induced via an Ornstein-Uhlenbeck process [66]. However, in Fig. 3b, we find that the population of all local structures that we consider (those pertinent to the Lennard-Jones model [58,67]) decreases with increasing activity for ρ = 0.87 in which we do not find any crystallisation and simply focus on the liquid local structure. This decrease in higher-order structure is in marked contrast to the previous work with the Wahsntröm model; however, the latter, a model glassformer was much more deeply supercooled and the dynamics, like the higher-order structure exhibited the opposite response to activity noted here, suggesting these two systems are in different regimes according to the categorisation introduced in ref. [65]. Our work shows qualitatively similar behaviour to passive systems when the temperature is increased [58,68,69].

Crystal growth and activity
We know from previous work [39,40] that monodisperse suspensions of active particles can crystallise at high density for activities below a certain Pe. The state diagram for this system is displayed in Fig. 1. We now consider in more detail the mechanism of crystal nucleation and polymorph selection at the examined state points.
In Fig. 4, we look at the size of some selected cluster populations as a function of Peclet number for two  Here, we emphasise the 10-membered defective icosahedron among the amorphous local structures detected by the TCC because it is a locally favoured structure in the hard sphere system [70,71]. The dashed vertical lines signal the transition from spinodal growth to nucleation and growth, followed by the transition to the fluid regime as Pe increases densities in the region of stability of the solid phase, ρ = 1.06 (a) and ρ = 1.15 (b). In particular, the fivemembered triangular bipyramid consists of two tetrahedra (the simplex for spheres in 3d). We also consider the defective icosahedron which is a locally favoured structure of the hard sphere system [70,71]. All curves are obtained by averaging the final state of the simulation runs over 10 independent trajectories. We observe the following common trends with increasing activity (Pe): For the passive case (Pe = 0), all trajectories crystallise into a mixture of FCC and HCP crystals, with a small preference for the FCC phase. This behaviour was observed in event-driven simulations of hard-spheres and is explained by finite-size structural fluctuations that favour FCC-rich nuclei compared to HCP-rich nuclei, due to the higher stacking entropy of cubic phases compared to hexagonal phases.
Crystallisation at these high densities occurs spinodally, i.e. it is characterised by the appearance of multiple nucleation events, and where crystal growth is controlled by the annealing of stacking faults. Spinodal crystallisation persists when activity is introduced in the system. Looking at the FCC and HCP populations, we observe that the effect of activity on polymorphism is to increase the fraction of FCC crystals with increasing Pe, at the expense of the HCP population, which decreases with increasing Pe. To explain the preference towards FCC, we recall that the formation of hardsphere crystals is subject to a mechanical instability under the effect of an external force which promotes the rearrangement of HCP layers into FCC layers [72]. This is confirmed in our simulations, where we observe the annealing of HCP stacking faults in favour of FCC environments promoted by the persistent motion of the active particles. Interestingly, the polymorph composition of the nuclei changes behaviour at a finite value of Pe: for example, for ρ = 1.06 (ρ = 1.15) the FCC population reaches a maximum in relative composition at Pe ∼ 3 (Pe ∼ 8).
This change of polymorphic behaviour coincides with a change in the crystallisation channel from spinodal to a nucleation-and-growth regime. In Fig. 4, the onset of the nucleation-and-growth regime is indicated with the dashed vertical line. Here, nucleation is a rare event Further increasing the activity causes the nucleation rate to drop, until crystallisation is no longer observed. In Fig. 4, this transition is represented with the vertical dash-dotted line where, not only do crystalline environments rapidly decay, but defective icosahedra environments increase to signify the transition to a fluid regime. At high Pe, this higher-order structure weakens, similar to the fluid case (ρ = 0.87) (Fig. 3b).
In Fig. 5, we focus on the state points displaying nucleation-and-growth and plot the mean first passage time t fp (n) , defined as the average elapsed time until the appearance of a nucleus of size n, at ρ = 1.06 and at three different Pe numbers. Over a wide range of n, t fp (n) can be fitted with the expression where k is the nucleation rate, n c is the critical nucleus size, erf is the error function, and c is a constant which in the equilibrium case (Pe = 0) is proportional to the curvature of the nucleation barrier ΔF at the critical size, c = ΔF (n c )/k B T . The curves show that the mean first passage time increases with increasing activity, with a consequent drop in the nucleation rates K extracted from the functional fits of Eq. 6 and plotted in the inset. Activity hinders nucleation, and even more so if the nucleation rates are scaled by the relaxation time in the fluid τ α which, as plotted in Fig. 2, drops faster than exponentially with increasing Pe. Interestingly, the critical nucleus size n c , as indicated on the right axis of Fig. 5 inset, also decreases with increasing activity. These critical sizes are considerably larger than what is typically observed in the passive case (at the same density), owing to the acceleration of the underlying fluid dynamics with activity, that allows the observation of longer nucleation induction times. In Fig. 6, we focus on the spinodal nucleation regime, i.e. when the non-equilibrium nucleation barrier is low enough for multiple nucleation events to occur simultaneously in the simulation box, and crystallisation is an activated process controlled by the rate of addition of new crystals on the nuclei. Panel (a) shows the fraction of crystalline particles as a function of time for Pe = 0, 2, 4 (continuous, dashed and dotted curves, respectively) and distinguishing between FCC (green curves) and HCP (black curves). For both FCC and HCP, we observe an increase in the crystal growth rate as a function of activity. The fraction of HCP at Pe = 4 shows a marked decrease at long times, which is due to the annealing of HCP stacking faults that we observed also in Fig. 4. To analyse the growth regime in panel (b), we fit the curves for the FCC phase with the Avrami equation.
where Y is the crystal fraction Y = (N c − N 0 )/(N − N 0 ), with N c the number of crystalline particles, N 0 the starting number of crystalline particles, and N the total number of particles. K = πkĠ 3 /3 is the Avrami constant proportional to the nucleation rate k and the growth rateĠ, and n is the Avrami exponent. From the fits in panel (b), we obtain n 1 and an increase in the growth constant K with increasing Pe. In panel (c) we show how the crystalline growth is rescaled by a characteristic time t Av = K −1/n (plotted in the inset). This timescale decreases with activity, signalling the increase in the growth rate with Pe. What appears to be an increase in the growth rate of nuclei in units of the active particles rotational time τ R is still a significant slowing down if measured instead in units of the relation time τ α .

Conclusion
We have considered the crystallisation behaviour of a suspension of active Brownian Particles that interact with a hard-sphere like interaction. We showed that the freezing line is strongly affected by activity and moves to higher densities as we increase the Pe number of the active particles consistent with previous work [12,39]. This is accompanied by a reduction in the nucleation rate, nucleation barriers and critical nucleus size with increasing activity. Despite the suppression of nucleation, the growth of nuclei is enhanced by the accelerated dynamics of the melt. This allows us to observe spinodal nucleation, where the growth is controlled by the rate of particle attachment, and thus speeded-up with activity. We observe a decrease in pair-and higher-order structure in the fluid with increasing activity. The former is compatible with certain dynamic regimes observed previously [63][64][65]. This is intriguing as one may enquire as to the nature of the higher-order structure approaching the MIPS phase boundary [40,41]. Very recently, comparisons have been made between MIPS and criticality in passive systems [41], and in the case of passive systems, approaching criticality, the population of higherorder structure detected by the TCC increases [60], in marked contrast to our findings here. In the future, it would be interesting to investigate whether the trend we have observed changes closer to the MIPS boundary or whether the response of the higher-order structure is profoundly different to passive systems.
Remarkably, activity also has a strong effect on polymorph selection. While the passive system crystallises in an equimolar mixture of FCC and HCP, active particles progressively favour the FCC phase at higher Pe. We observe this as annealing of HCP stacking faults, especially close to the crystal boundaries.