Active cluster crystals

We study the appearance and properties of cluster crystals (solids in which the unit cell is occupied by a cluster of particles) in a two-dimensional system of self-propelled active Brownian particles with repulsive interactions. Self-propulsion deforms the clusters by depleting particle density inside, and for large speeds it melts the crystal. Continuous field descriptions at several levels of approximation allow to identify the relevant physical mechanisms.


Introduction
The collective behavior of self-propelled particles is a fascinating topic both for its numerous applications and its intrinsic theoretical interest [1][2][3][4]. Many studies focused on the formation of clusters reporting two main different cases: for 'active crystals', self-propulsion leads to a modification of the properties of a pre-existing crystal that is usually induced by long-range attractive and short-range repulsive interactions [5][6][7]. On the other hand, for Mobility-Induced Phase Separation (MIPS), the system separates into two fluid phases of different densities. Clusters of the densest phase form by a purely non-equilibrium mechanism induced by the persistence of the motion of the particles, the nature of their interactions being not as crucial (this phenomenon was observed for various repulsive and attractive forces [8][9][10]). Methodologically, both mechanisms can be studied by considering an effective density-dependent velocity replacing the two-body interacting potential [11][12][13][14][15][16][17].
We analyze in this paper a different case of cluster formation with active objects (some additional models and experiments on clustering of self-propelled particles may be found, for example, in [10,18,19]). It is the nonequilibrium counterpart of the so-called cluster crystals [20][21][22][23][24][25], which appear in equilibrium systems interacting with soft-core repulsive potentials, and are solid-like structures where the unit cell is occupied by a closely packed cluster of particles. Here, clustering appears under a repulsive potential so that it is not a consequence of purely local effects but rather involves a global minimization of energy by intercluster interactions. Namely, for sufficiently soft repulsive potentials, repulsion from the neighboring clusters can exceed the intracluster repulsion and lead to an effective particle confinement in the cluster. For the activeparticle case discussed here, this mechanism is very different from the other cases of active crystals and it is unlikely that previous local arguments derived for active crystals or MIPS can describe it.
Similarly to previous works [5,6], we consider the simplest extension of an equilibrium system of repulsive Brownian particles in two dimensions by providing them with an internal degree-of-freedom, the orientation of a self-propulsion speed. We will address the following questions: Do active cluster crystals (ACC) with only repulsive interactions exist? Is the structure of the clusters modified by activity? Can we find continuum equations describing this active system? The answer to the first question has been partially given by observations of cluster crystals in [26,27], for example. But in these references the focus is in deformable-body and alignment interactions, so that the specific role of repulsive interactions needs further clarification.
To do this and answer the remaining questions we first present numerical simulations of self-propelled soft repulsive particles showing that ACC can be observed but that an increasing self-propulsion eventually leads to their destruction. For small diffusion, empty clusters are found as the particles tend to accumulate on their edges.
Then we provide a continuum field description of this interacting particle system and analyze some of its predictions.

Numerical results and pattern formation
We start with the Langevin equations of a two-dimensional overdamped system of active Brownian particles interacting in pairs via a potential V x ( ): ) and i q are respectively the position and orientation of the particle i, n U 0ˆi s the self-propelling velocity, of constant modulus U 0 and direction given by the unit vector n cos , sin ). The particles are subjected to Gaussian translational and rotational noises, , . D and D r are the translation and rotational diffusion coefficients. If they originate from the same thermal bath they are related to temperature by fluctuationdissipation relationships. But under general non-equilibrium conditions they can have different origins and then we will assume here that they are independent parameters.
Cluster crystals appear at equilibrium (U 0 0 = ) under repulsive soft-core potentials which have negative Fourier components [23,28]. A convenient class of such potentials is the generalized exponential model of , 0  > , and R is the interaction range. The Fourier transform V q ( ) of GEM-α potentials with 2 a > takes both positive and negative values, so that cluster crystals would appear in that case for sufficiently small D [23,28] and/or large average densities r . Here we consider 3 a = .
Besides the global parameters N and L, which fix the mean density N L r = , our model contains five parameters: D, D r , U 0 , and the two parameters in the potential (ò and R). The first parameter compares the strength of self-propulsion with that of translational diffusion, the second the strength of the interaction forces with selfpropulsion, and the third the persistence length with the interaction length. Unless otherwise stated, in the following we will fix the parameters in the potential (to 0.0333  = and R=0.1) and D r =0.1 (this last one being equivalent to fixing the units of time), and describe the behavior of the system by varying U 0 and D. Other ways to explore the transitions occurring in the system are possible, and they can be easily related to the one described here simply by looking at the dimensionless quantities stated above. For example, in the dimensionless parameters the dependence on D r is reciprocal to that on U 0 , meaning that the changes in behavior described below when increasing U 0 at D r fixed will also be observed when decreasing D r at U 0 fixed.
We integrate numerically (1) with the Euler-Maruyama method [29]. For small U 0 we observe a transition, when increasing the average density r or decreasing the diffusion coefficient D, from a homogeneous distribution of particles to a statistically steady hexagonal crystal of clusters (see figures 1(a) and (b)), similar to the behavior observed for passive particles (U 0 0 = ). If U 0 and D are small enough, there is no particle exchange between clusters. If U 0 is further increased, the cluster crystal remains but jumps between clusters eventually occur. For high values of U 0 clusters are destroyed and the steady-state becomes statistically homogeneous (see table 1 for precise values). Contrary to the passive case for which the distribution of particles inside a cluster is Gaussian to a good approximation [28], active particles tend to stay close to the edges. This is particularly obvious for small D and high U 0 . In that case, the clusters have a ring shape with their centers empty (see figure 1(b)). We confirmed this by computing the distribution of r D , the relative distance between the position of a particle and the center of its cluster. Figure 2(a) shows that when U 0 is large and D small, there is a clear depletion close to r 0 D = while the distributions look Gaussian for higher D (remember than the relevant dimensionless parameter is U D D r 0 , so that simple Gaussian clusters are also obtained for large enough D r ). A similar cluster-center depletion [30][31][32] has been observed in active systems in which aggregation occurs by confinement in an external potential. Our studies in section 3.2 below will show that indeed both phenomena are related. They arise from a purely non-equilibrium effect caused by the persistence of the particle velocity which disappears for small U 0 or large D r . Another consequence of persistence is that the orientation n of the particles inside a cluster is clearly radial, pointing to the exterior (inset of figure 1(b)). For larger D like in figure 1(a), this radial orientation of velocity vectors remains predominant but not as strong. From the definition of the model, and as we will show later on, the equilibrium Brownian non-active dynamics, for which clusters are no longer empty, is achieved for large D r .

Derivation and stability analysis
To have some analytical handle on the phenomena above we derive now the equivalent of a DK equation [28,33,34] for active particles. Note that a continuum description equivalent to the particle system should include a noise term arising from the random discrete particle dynamics. However, as in [28] we average out   (1). The function is an Itō process so that we can apply the Itō formula [35]. Using integration by parts we obtain the following deterministic equation for the average which is the expected value of the density of particles at location x and with orientation θ: )is the total particle density at point x¢ and n cos , sin ) is the unit vector in the θ direction. Equations similar to (2) have already been derived by other means to describe MIPS [16,36] but here we do not approximate the non-local interaction term in terms of a local density-dependent velocity, and instead we have introduced a mean-field approximation: ). For passive particles interacting with soft potentials, this is usually a good approximation because of the large number of particles within the interaction range [28].
An equivalent representation is obtained by introducing the angular Fourier transform: (2) we find a coupled set of equations for the angular modes n r : t n n n r n x n x n y n y n is the force induced at point x by the particle interactions. The homogeneous and isotropic state given by    = + ( ). The idea of self-propulsion being equivalent to an enhanced diffusion given by this D eff has been pointed out by several authors in the dilute limit [12,37,38], but we have not made any assumption regarding the density of the system so that our expressions remain valid at high densities. As mentioned, the Brownian dynamics is obtained for very large D r , or rather very small U 0 . The behavior for q 0  guarantees that for repulsive V x ( ) (for which ) there is no long-wavelength spinodal-decomposition-like instability as in the case of MIPS. Instead we can only have finite-wavelength instabilities. Note that our expressions for the instability of the homogeneous state can be used to elucidate if particular potentials beyond the GEM class explicitly discussed here lead to crystal formation, as long as the potential is sufficiently soft to justify the mean-field approximation which leads to equation (2).
We can now check the accuracy of truncating at different values of M. If enough modes M have been included, the predictions obtained for higher values of M should collapse on the same curve. It appears that the higher the self-propulsion, the higher M needs to be to reach convergence: indeed, figure 3 shows that M=3 modes are enough for U 0 =1.75, but we need at least M=5 modes for U 0 =2.4. Note that for U 3 0 = , M=13 modes are barely enough. Therefore, keeping only a small number of modes will accurately describe the structural transition only if this one takes place for a small value of U 0 , meaning that the system is already quite close to the critical point. However, we will see in the following subsection that they contain qualitatively the relevant mechanisms shaping the cluster crystals. Figure 4 shows the maximum real part of the growth rate as a function of q for several values of U 0 . The continuous fraction defining f q , l ( )has been truncated to M=13, which is sufficiently large for the results to remain unchanged when increasing M further. In the passive case [28] and for the GEM-3 potential the homogeneous state becomes first unstable with respect to a wavenumber q c depending only on the interaction range R and satisfying q R 5.0 c » . This corresponds to a cluster crystal of periodicity q R 2 1 . 2 6 c p » . For the parameters used here, for which R=0.1, we have q 50 c » . In figure 4 we see that a very similar wavenumber is the most unstable in the active case for small U 0 . When U 0 is sufficiently increased though, these positive growth rates eventually become negative. Activity thus stabilizes the homogeneous state and prevents clustering, in agreement with our particle simulations.
More quantitatively, one can compare the values of the transition threshold U 0 c obtained analytically for sufficiently large M, and from particle simulations. For the second case we start from an initial hexagonal crystal of clusters and increase U 0 step-by-step until we get a long-lived gas phase. Table 1 shows that the linear stability analysis gives us the right order of magnitude but it systematically underestimates the threshold. We attribute this to the nature of the transition, which in the passive-particle case is subcritical and subject to hysteresis [28]. We expect this subcritical nature to remain also for small U 0 . Also, the use of a deterministic version of the DK equation, as an approximation to the complete stochastic one, may be a source of error.

macroscopic equation for 2 modes
Equations (2) or (3) are quite involved and it would be desirable to work with a simpler set of equations such as a small-M truncation of (3). The linear analysis reveals that this is in general not an accurate approach, as values up to M=13 need to be considered to have convergence of the growth rates. This is so because, except for the largest values of D r or the smallest U 0 , the average velocity field at each point has a well-defined direction (see for example figure 1(b)) and many angular modes are needed to represent such localized distribution in θ space. Nevertheless a truncation to order M=1 leads to relatively simple equations which give insight into the physical mechanisms and favor qualitative understanding. Following Bertin et al [39], we consider a hypothetical situation in which )ˆ( ) is the momentum or polarization field. The product F P r appearing in (10) is a tensor product. Even if equations (10) are a rough approximation to (3) they still describe the system behavior qualitatively. Integrating numerically equation (10) with a spectral method using 512×512 grid points, we recover the crystal of clusters for small values of U 0 (see figure 1(c)) and a homogeneous state for higher U 0 . The polarization structure of the clusters is also in good agreement with our previous observations: figure 1(d) shows that the polarization field inside the clusters is radial. Empty clusters are also observed in one dimension for small D and high enough U 0 . Having to work at small values of D makes it more difficult to obtain empty clusters in two dimensions because of the need of a high numerical resolution. We have shown however than equations (10) support clusters with a depletion in their center (see figure 2(b)) if we replace the convolution product defining F r by an effective confinement potential justified by the approximation described in the next paragraph. Equations (10) serve also as a starting point to better understand the mechanisms shaping cluster structure. As in the passive case [28] we can focus on the case of small D so that a first approximation for the steady crystal density is a set of delta functions at the lattice points a i { }: . The number of particles per cluster, N p , can be expressed in terms of r and the intercluster distance a: N a 3 2 To consider the structure of a narrow cluster centered at x 0 = we keep in the lattice sum only the central and the six neighboring clusters, and expand around x 0 » . For a GEM-α potential with 2 a < the dominant term is the interparticle repulsion within the central cluster, so that F x r ( ) points outwards and, as it is observed, the aggregate disappears. But when 2 a > the repulsion from the neighboring clusters prevails and produces a confining effective force at x 0 » overcoming the local repulsion.
At first order in the distance x to the cluster center, we have F x x g »r ( ) . γ depends on N p , on the intercluster distance a and on the interaction potential parameters ò and α. Its precise expression is not particularly illuminating, but it can be found as To get the equation for ρ, a first integral has already been performed under the condition of zero net particle flux in or out of the cluster, as appropriate for the steady state. The system (11) should be solved in r 0, ], but the approximations used are valid only if it gives a cluster width much smaller than the interaction range R. The first equation in (11) shows that in the limit of small U 0 , this cluster width is of the order of w D g » . Equations (11) require three initial or boundary conditions, which could be taken as the values of ρ, p r and p r ¢ at which shows that there is a shape change, in agreement with the particle simulations, from a maximum of density at the origin ( 0 0 r < ( ) ) to a minimum ( 0 0 r > ( ) ) when the polarization slope at the origin p 0 r n = ¢ ( ) changes from U 0 0 n gr < ( ) to U 0 0 n gr > ( ), respectively. Finally, the condition P x x  r | ( )| ( )-resulting from the definition of P -imposes the value of this slope: ν must indeed cancel the prefactor of the slow decay at large r, p r r~m with D U D 2 2 r 0 2 m g g = --( ) ( ), arising from the large-r behavior of the hypergeometric function that solves the homogeneous part of the (11) for p r (r). This determination can only be done numerically. Figure 5 shows examples of solutions of the radial equation (11), displaying the two qualitatively different cluster shapes, with smaller/larger density at the center. Although the truncation to M=1 leading to equations (10) and (11) precludes quantitative agreement with particle simulations, figure 5 (see also figure 2(b)) shows that this 2-mode truncation contains the essential qualitative mechanisms needed for the cluster-crystal formation and the structure of the clusters: an effective potential which confines particles within clusters, originating from the repulsion by the neighboring clusters, and the presence of a polarization field P , increasing with U 0 , which pushes the particles towards the periphery of the clusters until destroying them.

Conclusion
We have shown that ACCs occur in systems of active Brownian particles with soft repulsive interactions. Selfpropulsion deforms the clusters by depleting particle density inside, and large self-propulsion stabilizes the homogenous-isotropic state. We have derived a continuous description and analyzed the crystal forming instability by linear analysis. Truncation to two angular modes, despite not being quantitatively accurate, retains the basic mechanisms. In particular it allows to understand crystal persistence as an effect of the confining forces arising from repulsion by neighboring clusters, and the cluster shape as a balance between the confining force and the tendency to radial escape driven by the polarization field.