A mesoscopic field theoretical approach for active systems

We introduce a mesocopic modeling approach for active systems. The continuum model allows to consider microscopic details as well as emerging macroscopic behavior and can be considered as a minimal continuum model to describe generic properties of active systems with isotropic agents. The model combines aspects from phase field crystal (PFC) models and Toner-Tu models. The results are validated by reproducing results obtained with corresponding agent-based microscopic models. We consider binary collisions, collective motion and vortex formation. For larger numbers of particles we analyze the coarsening process in active crystals and identify giant number fluctuation in the cluster formation process.


Introduction
Active systems can exhibit a wide range of collective phenomena. Depending on the density and interaction details, energy taken up on the microscopic scale can be converted into macroscopic, collective motion. Theoretically, such behavior can be addressed either from the microscopic scale, taking the interactions into account or from the macroscopic scale, focusing on the emerging phenomena. For reviews on both theoretical descriptions see e.g. [1,2,3]. We will here introduce a mesoscopic modeling approach which combines aspects from both scales.
Vicsek-like models [4] are examples for the microscopic viewpoint. These models consider particles, which travel at a constant speed to represent self-propulsion, whose direction changes according to interaction rules which comprise explicit alignment and noise. However, explicit alignment rules are not necessary for the emergence of collective phenomena. For elongated particles already the shape gives rise to alignment mechanisms due to steric interactions, see e.g. [5,6], and even for circular or spherical particles collective phenomena can be observed if inelastic behavior in the interaction rules is considered, see e.g. [7]. Besides these agent-based microscopic approaches also field theoretical macroscopic modeling of active systems has been widely applied. Most approaches are based on the Toner-Tu model [8] and consider only orientational ordering. Collective motion has been addressed within such models, see e.g. [9]. More detailed continuum modeling approaches which address besides orientational ordering also positional ordering and thus include also aspects from microscopic models are rare and have so far only been considered for active crystals [10,11]. Such systems arise for high densities, where particle interactions dominate the propulsion. In this article we propose an extension of the model considered in [10] to allow also for low densities. The proposed mesoscopic field theoretical approach can be considered as a minimal continuum model to describe generic properties of active systems with isotropic agents. After introducing the model, we validate it by reproducing results obtained with corresponding agent-based microscopic models [7]. We consider binary collisions, collective motion and vortex formation. Considering larger systems, which can be accessed using the introduced mesoscopic field theoretical approach, the formation of collective motion can be analyzed. For high densities we observe a coarsening process of regions of different directions of collective motion. This was already mentioned in [10], but not analyzed. In a broader context the observations can also be related to defects in active crystals. For orientational ordering this was e.g. analyzed in [12]. Another remarkable property of active systems is giant number fluctuations. In contrast to equilibrium systems, where the standard deviation ∆N in the mean number of particles N scales as √ N for N → ∞, in active systems ∆N can become very large and scales as N α , with α an exponent as large as 1 in two dimensions [8]. This theoretical prediction is often associated with elongated particles and a broken orientational symmetry [13,14,15,16], but it has also been verified in simulations of agent-based models for disks with no-alignment rule, see [17], and was demonstrated by experiments and simulations in [18]. We use large scale simulations to show giant number fluctuations in the proposed mesoscopic field theoretical approach.

The model
Our starting point is the active phase field crystal model derived in [10,11] for a one-particle density field ψ(r, t), which is defined with respect to a reference densitȳ ψ. The parameter r is related to an undercooling and together withψ determines the phase diagram. The functional arises by splitting the Helmholtz free energy in an ideal gas contribution and an excess free energy, rescaling and shifting ψ, expanding the ideal gas contribution in real-space and the excess free energy in Fourier-space, and simplifications by removing constant and linear terms that would vanish in the dynamical equations. A detailed derivation of the energy and its relation to classical density functional theory can be found in [22,23]. The second quantity in eq. (1) is the polar order parameter P (r, t) and the remaining parameters are: M 0 mobility, v 0 selfpropulsion determining the strength of the activity, α 2 and α 4 two parameters related to relaxation and orientation of the polarization field, and C 2 and C 4 are parameters which govern the local orientational ordering. The model is used in [10,11] to study crystallization in active systems. To allow for a description of individual particles, we consider a variant of the PFC model, the vacancy PFC (VPFC) model, introduced in [24,25] and also considered in [26]. Instead of the Swift-Hohenberg functional eq. (2) we consider a density field with positive density deviation ψ only, which leads to a modification of the particle-interaction and allows to handle single particles, as well as many individual particles. The new energy functional F vpfc is: with a penalization parameter H. As in other, more coarse-grained models for active systems [3] we use a classical transport term with advection velocity v 0 P for the local density field ψ. This modification turns out to be more stable in comparison to the term used in eq. (1), if considered for individual particles. The second modification ensures the polar order parameter P to be a local quantity that is different from zero only inside the particles. The new set of dynamical equations we obtain is: with β a parameter, which is typically larger than the other terms entering the P equation.

Model verification
The aim of this Section is to show that the mesoscopic field theoretical model eq. (4) can be used to simulate active particles and allows to recover known phenomena. The versatility of the model thereby allows us to apply it to different physical situations that have been previously studied using agent-based models, see e.g. [7,27]. We first consider the situation of one particle, followed by studying the interaction of two particles. These simple computations allow for detailed parameter studies. We found no qualitative difference in the results of our simulations when the parameters C 2 and C 4 are set to zero. Therefore, we simplify our model by restricting ourselves to the case C 2 = C 4 = 0, which allows only gradients in the density field ψ to induce local polar order. The remaining parameters to specify are M 0 , v 0 , α 2 , α 4 , β, H and r. The second test concerns the emergence of collective phenomena in confined geometries with systems with a number of particles 100. All the simulations are in two dimensions. The maxima in the local one-particle density field ψ(r, t) are always tracked for post-processing and evaluation, and every maximum is interpreted as a particle (see [28] for a discussion about the validity of this interpretation). The particle's velocity is computed as the discrete time derivative of two successive maxima.

Onset of movement and particle shape
We at first want to understand what happens in a minimal system, where a single active particle is free to move. In particular, we are interested to know if there is a critical value for the activity v 0 required for the onset of movement, as it has been observed in [10] for active crystals.
In Fig. 1(a) the particle velocity is plotted as a function of the activity v 0 and we can see that for small activities (v 0 < 0.5) the particle does not move at all. After a certain threshold value v 0 0.5 the particle starts to move with a constant velocity, which approximately linearly increases for increasing v 0 . This is exactly the same behavior as observed in [10] for the sample-averaged magnitude of the crystal peak velocities.
An important new feature of our model has to do with the mobility term M 0 entering eq. (4). M 0 = 1 is used in [10]. This value would lead to strong numerical instabilities for the modified model eq. (4). Larger values for M 0 can suppress this numerical instability. While the mobility does not particularly change the particle velocity, we observe that M 0 directly influences the shape of the particle: for small (but still greater than 1) M 0 the particle shows an elliptic form, whereas further increasing M 0 restores a circular shape for the particle. The dependency of the particle shape on the value of M 0 is shown in Fig. 1(b) for different values of v 0 above the threshold value.

Binary collisions and elastic deformation
The study of binary collisions between particles is often used as a benchmark problem to predict how larger systems evolve, see e.g. [27,29]. In particular it has been observed [7] that completely inelastic collisions lead to a force that aligns the particles direction. We here consider only perfectly symmetric collisions meaning that the incidence angle and the initial velocity are the same for both particles. Different particle trajectories obtained using different mobility M 0 and activity v 0 are shown in Fig. 2(a) and (b). The elastic deformation of a single particle during a collision can be seen in Fig. 2(c), where the eccentricity is plotted as a function of time, whereas a time series of a single collision is shown in Fig. 2(d). Collisions of deformable particles have also been considered in agent-based models [30], with a qualitatively similar behavior. However, the deformations in our approach strongly depend on M 0 and are negligible for large values. We therefore do not analyze this effect further and interpret the particles as being spherical in the coming simulations, which are all done for large M 0 .
All results indicate the particle alignment to be not instantaneous. There is an initial oscillatory phase, whose length and magnitude depends on M 0 and v 0 . Small activity and large mobility lead to an almost instantaneous alignment, whereas large activity and small mobility lead to oscillations for a certain period of time, before the particles finally align and travel together.

Collective motion in confined geometries
As already analyzed using agent-based, e.g. [7,30], and phase field models [29,31] the interaction of a moderate number of particles can lead to collective motion. We here consider simulations with 100 particles to recover these results. To analyze the phenomena we define the translational order parameter φ T and the rotational order parameter φ R as wherev i (t) is the unit velocity vector of particle i at time t,ê θ i (t) = (− sin(θ i (t)), cos(θ i (t))) is the unit angular direction vector of particle i at time t, and n is the number of particles. In case of translational migration we obtain φ T = 1 and rotational migration leads to φ R = ±1.

Collective migration
We consider a square domain with periodic boundary conditions, i.e. simulating an infinite plane where particles are free to move without obstacles. Fig. 3 shows the resulting behavior: At the beginning (Fig. 3(a)) there is no specific order and particles move towards different directions. After the first collisions take place some particles start to align with each other and small blocks of particles, in which particles orient in the same direction are formed (Fig. 3(b)). If these blocks collide they change their direction until all particles in the system are traveling in the same direction (Fig. 3(c)). This behavior is further confirmed by the translational order parameter φ T 1 after a certain time, see Fig. 3(d). We also observe that this migration state is not particularly affected by the value of the mobility M 0 .

Vortex formation and oscillatory motion
We now consider a confined geometry and specify ψ = 0 and P = 0 at the boundary, which serve as an approximation for reflecting boundary conditions. The first geometry is a disk. Here again at the beginning the particles move in a chaotic way, see Fig. 4(a). Then a vortex is formed and most of the particles follow an anticlockwise trajectory. In the center some particles move in the opposite direction, see Fig. 4(b). Eventually also these particles are forced to align with the rest of the system, see Fig. 4(c). This behavior is confirmed by the rotational order parameter φ R , which is shown as a function of time in Fig. 4(d). Similarly to the collective migration case studied above, the mobility M 0 does not play a major role in the formation of the vortex.
Ellipses provide a more interesting geometry and confining active particles inside them can give rise to different kind of collective phenomena, where the ellipse aspect ratio A/B, where A, B are the length of the semi-major and semi-minor axis, respectively, plays an important role. An ellipse with a small A/B ≤ 3 show a similar behavior as the disc shape. Such geometries produce once again a vortex, where the particles move along the boundaries. More elongated shapes with A/B = 10 dramatically change the behavior. As already shown in [7] particles move collectively along the major axis with oscillating direction. The same behavior could be observed with our model, see Fig. 5. All particles move in one direction until they hit the high curvature region. This produces an impulse that propagates fast along the whole system and reverts the direction of the particles. The whole process repeats every time a boundary is  reached and an oscillatory motion is the result. To obtain this result and ensure a constant particle number it is necessary to use a particularly high value for the mobility, M 0 = 500.

Validation and numerical issues
These examples demonstrate the validity of our continuous modeling approach. All known qualitative properties which have been shown using agent-based simulations could be reproduced. Until now a sequential finite element approach has been used to solve the evolution equation. For larger systems we must work in a parallel environment with multiple processors. We adopt a block-Jacobi preconditioner [32,33] that allows us to use a direct solver locally. The approach is implemented in AMDiS [34,35] and shows good scaling properties, which allows to consider systems with 15, 000 particles on the available hardware.

Results
We again consider a system with periodic boundary conditions for which collective migration and cluster formation is expected also for larger numbers of particles. Average and standard deviation of the data for 6 different simulations started with different initial conditions and a tolerance parameter equal to π/10 are shown in blue. Each of the shaded curves has been obtained as the average of the six simulations, but using different values for the tolerance parameter, ranging from π/8 (purple curve) to π/12 (red curve). Simulation parameters are (v 0 , M 0 ) = (2,60) However, how the collective migration state is reached is not well understood and will be analyzed in detail. We start with a system with high volume fraction. The volume fraction is defined as φ = nσ/|Ω|, with number of particles n, domain size |Ω and σ the area occupied by a single particle, which is equal to σ = π(d/2) 2 , with d = 4π/ √ 3 the lattice constant determined by the free energy eq. (2). For φ > 0.6 the system shows behavior of active crystals. Figure 6(a) shows various snapshots of the evolution with regions of particles moving in the same direction color coded. The regions can be identified as active grains, which undergo a coarsening process. The black particles determine orientational defects. They are identified as particles where the change in orientation from one particle to its neighbors is above a certain threshold. The number of black particles certainly depends on the threshold value, however its decrease is independent on the value. The data is not sufficient to identify a scaling law. However, the robustness of the coarsening process is shown in Fig. 6(b).
If we decrease the volume fraction the behavior changes. For very low density, Here the color coding of each particle corresponds to the cluster it belongs to initially, which highlights to dynamics during cluster formation. φ = 0.03, Fig. 7(a) shows the tendency of the particles to group together, but the formed clusters are very small and there remain many isolated particles in the system. Increasing the density, φ = 0.12, increases the size of the formed clusters, but the average size of a cluster remains very small if compared to the total number of particles, see Fig. 7(b). Further increasing the density, φ = 0.25, leads to the formation of large mobile clusters and a drastic reduction of the number of particles which do not belong to any cluster, see Fig. 7(c). This behavior is similar to the results in [36] for (quasi-) two-dimensional colloidal suspensions of self-propelled particles. For these systems we compute the standard deviation ∆N as a function of the mean number of particles N . For active systems it is theoretically predicted that ∆N scale as N α , with 0.5 < α ≤ 1 and giant number fluctuations occur if α ≈ 1, see [8]. We compute α by considering different subregions of our computational domain. The results for 600 particles are shown in Fig. 8, demonstrating an increase of α with increasing φ, with the largest value reached being α = 0.79. This continuous increase in α, as well as the obtained values are consistent with the behavior found in [17] for moderate numbers of particles but larger volume fractions. However, there is a significant difference, the formed clusters in [17] are stationary, our clusters are mobile, similar to the light activated living crystals in [18]. The experimental results as well as the simulations in [18] lead to similar values of α as ours already for smaller φ.

Summary
In summary, we extended the phase field crystal model for active crystals [10,11], which combines the classical phase field crystal model of Elder et al. [19,20] with the approach of Toner and Tu [8], by considering individual active particles. This can be realized by penalizing a locally vanishing one-particle density, as considered in [24]. The resulting mesoscopic field theoretical model has been validated against known results obtained with minimal agent-based models. We found a threshold value for the activity, necessary to induce motion for a particle. Collective motion and vortex formations have been identified, as well as oscillatory motion, depending on the considered confinement. All these results are in agreement with the results in [7]. For larger systems, we analyze the formation of a traveling crystal if prepared from an initially disordered state. The traveling crystals emerge through a coarsening process from a multidomain texture of domains traveling collectively in different directions. For lower volume fractions we could identify giant number fluctuations. As theoretically predicted the standard deviation ∆N scales as N α for active systems. The computed exponent α as a function of volume fraction is in agreement with experimental and simulation results obtained for light activated colloidal particles [18].
The proposed mesoscopic field theoretical model can be extended from two to three spatial dimensions. Other possible extensions consider binary mixtures and hydrodynamic interactions, which are already considered within the phase field crystal model for passive systems, e.g. in [37] and [38,39,26,40], respectively. Together with efficient numerical algorithms, see e.g. [32,41], this provides the possibility to study emerging macroscopic phenomena in active systems with microscopic details, e.g. to validate coarse grained approaches, as considered in [42,43].