A microscopic field theoretical approach for active systems

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–3]. We will here introduce a continuum 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 also for circular or spherical particles collective phenomena can be observed if inelastic behavior in the interaction rules is considered, see e.g. [7]. In [8] it has been shown that even deformation of self-propelled particles can lead to alignment, without any explicit alignment rule. Besides these agent-based approaches also microscopic phase field models have been proposed in the context of collective cell migration [9, 10]. Each cell is thereby approximated by a deformable phase field variable and the physical processes behind cell motility are considered using an active polar gel theory [11–14]. Due to the complexity of these models the number of cells, which can be considered numerically is limited. However, collective alignment could be identified also with these models as a result of inelastic interactions.

Field theoretical macroscopic modeling of active systems allows us to describe emerging phenomena on larger scales. Most approaches are based on the Toner-Tu model [15] and consider only orientational ordering. Collective motion has been addressed within such models, see e.g. [16]. 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 [17, 18]. Such systems arise for high densities, where particle interactions dominate the propulsion.

In this article we propose an extension of the model considered in [17] to also allow for low densities. The proposed microscopic field theoretical approach can be considered as a minimal continuum model to describe generic properties of active systems and a coarse grained model of the detailed descriptions in [9, 10].

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 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 [17], 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 [19]. Another remarkable property of active systems is giant number fluctuations. In contrast to equilibrium systems, where the standard deviation ${\rm{\Delta }}N$ in the mean number of particles N scales as $\sqrt{N}$ for $N\to \infty$, in active systems ${\rm{\Delta }}N$ can become very large and scales as ${N}^{\alpha }$, with α an exponent as large as 1 in two dimensions. This theoretical prediction is often associated with elongated particles and a broken orientational symmetry [20–23], but it has also been verified in simulations of agent-based models for disks with no-alignment rule, see [24], and was demonstrated by experiments and simulations in [25]. We use large scale simulations to show giant number fluctuations in the proposed microscopic field theoretical approach.

Our starting point is the active phase field crystal model derived in [17, 18], which reads in scaled units

$$\begin{eqnarray}{\partial }_{t}\psi & = & {M}_{0}{\rm{\Delta }}\displaystyle \frac{\delta {{ \mathcal F }}_{{\rm{pfc}}}}{\delta \psi }-{v}_{0}\triangledown \cdot {\boldsymbol{P}}\\ {\partial }_{t}{\boldsymbol{P}} & = & {\rm{\Delta }}({\alpha }_{2}{\boldsymbol{P}}+{C}_{2}{{\boldsymbol{P}}}^{3})-({\alpha }_{4}{\boldsymbol{P}}+{C}_{4}{{\boldsymbol{P}}}^{3})-{v}_{0}\triangledown \psi ,\end{eqnarray} \tag{ 1 }$$

see the supplement in [17] for a detailed non-dimensionalization of the equations, which can be derived from microscopic dynamic density functional theory. The model combines the phase field crystal (PFC) model, introduced in [26, 27] to model elasticity in crystalline materials, with a polar order parameter and a self-propulsion term. The energy functional ${{ \mathcal F }}_{{\rm{pfc}}}$, the PFC model is based on, is a Swift-Hohenberg energy [28]

$$\begin{eqnarray}&&{{ \mathcal F }}_{{\rm{pfc}}}=\int \left[\displaystyle \frac{{\psi }^{4}}{4}+\displaystyle \frac{1}{2}\psi (r+{(1+{\rm{\Delta }})}^{2})\psi \right]d{\boldsymbol{r}},\end{eqnarray} \tag{ 2 }$$

for a one-particle density field $\psi ({\boldsymbol{r}},t)$, which is defined with respect to a reference density $\bar{\psi }$. The parameter r is related to an undercooling and together with $\bar{\psi }$ 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 [29, 30]. The second quantity in equation (1) is the polar order parameter ${\boldsymbol{P}}({\boldsymbol{r}},t)$, which is related to a coarse-grained velocity field with a typical magnitude v₀ of the self-propulsion velocity. The remaining parameters are: M₀ mobility, v₀ self-propulsion determining the strength of the activity, ${\alpha }_{2}$ and ${\alpha }_{4}$ two parameters related to relaxation and orientation of the polarization field, and C₂ and C₄ are parameters which govern the local orientational ordering. The model is used in [17, 18] 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 and analysed in [31–33] and also considered in [34]. Instead of the Swift-Hohenberg functional equation (2) we consider a density field with positive density deviation ψ only, which leads to a modification of the particle-interaction and allows us to phenomenologically describe single particles, see [33] for a detailed analysis. The new energy functional ${{ \mathcal F }}_{{\rm{vpfc}}}$ is:

$$\begin{eqnarray}&&{{ \mathcal F }}_{{\rm{vpfc}}}={{ \mathcal F }}_{{\rm{pfc}}}+\int H(| \psi {| }^{3}-{\psi }^{3})d{\boldsymbol{r}},\end{eqnarray} \tag{ 3 }$$

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}{\boldsymbol{P}}$ for the local density field ψ. This modification is more general and turns out to be more stable in comparison to the term used in equation (1), if considered for individual particles. The second modification ensures the polar order parameter ${\boldsymbol{P}}$ to be a local quantity that is different from zero only inside the particles. This allows us to interpret ${\boldsymbol{P}}$ in the context of cells as a coarse grained orientation of actin filaments, similar to [9, 10]. The new set of dynamical equations we obtain is:

$$\begin{eqnarray}{\partial }_{t}\psi & = & {M}_{0}{\rm{\Delta }}\displaystyle \frac{\delta {{ \mathcal F }}_{{\rm{vpfc}}}}{\delta \psi }-{v}_{0}\triangledown \cdot (\psi {\boldsymbol{P}})\\ {\partial }_{t}{\boldsymbol{P}} & = & {\rm{\Delta }}({\alpha }_{2}{\boldsymbol{P}}+{C}_{2}{{\boldsymbol{P}}}^{3})-({\alpha }_{4}{\boldsymbol{P}}+{C}_{4}{{\boldsymbol{P}}}^{3})-{v}_{0}\triangledown \psi -\beta {\boldsymbol{P}}{{\mathbb{1}}}_{\psi \leqslant 0},\end{eqnarray} \tag{ 4 }$$

with β a parameter, which is typically larger than the other terms entering the ${\boldsymbol{P}}$ equation.

The aim of this section is to show that the microscopic field theoretical model equation (4) can be used to simulate active particles and allows the recovery of 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]. 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₂ and C₄ are set to zero, providing ${\alpha }_{2}\gt 0$ and ${\alpha }_{4}\gt 0$. Therefore, we simplify our model by restricting ourselves to the case ${C}_{2}={C}_{4}=0$, which only allows gradients in the density field ψ to induce local polar order. M₀ and v₀ will be specified for each simulations. The other parameters are $({\alpha }_{2},{\alpha }_{4},\beta ,H,r)=(0.2,0.1,2,1500,-0.9)$ unless otherwise specified in the figure captions. The second test concerns the emergence of collective phenomena in confined geometries with systems with a number of particles $\simeq 100$.

All the simulations are in two dimensions. We consider a finite element approach to solve the evolution equations. An operator splitting approach is used, where the equation for ψ is discretized according to [35] and the equation for ${\bf{P}}$ using a semi-implicit Euler-scheme with all nonlinear terms linearized around the value of the last time step. The maxima in the local one-particle density field ψ are always tracked for post-processing and evaluation, and every maximum is interpreted as a particle. The morphology of the particle is obtained from a contour line of the density field, see below for details, and its velocity is computed as the discrete time derivative of two successive maxima. The computational domain varies for the different examples and is specified below. The initial condition is chosen for ψ using a one-mode approximation for each particle, with the centers placed randomly, see [34] for details. Only for high densities a perturbation from a hexagonal ordered state is chosen as initial condition. The ${\bf{P}}$ field is set to zero initially.

3.1. 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₀ required for the onset of movement, as it has been observed in [17] for active crystals.

In figure 1(a) the particle velocity is plotted as a function of the activity v₀ and we can see that for small activities (${v}_{0}\lt 0.5$) the particle does not move at all. After a certain threshold value ${v}_{0}\simeq 0.5$ the particle starts to move with a constant velocity, which approximately linearly increases for increasing v₀. This is exactly the same behavior as observed in [17] for the sample-averaged magnitude of the crystal peak velocities.

Zoom In Zoom Out Reset image size

An important new feature of our model has to do with the mobility term M₀ entering equation (4). ${M}_{0}=1$ is used in [17]. This value would lead to strong numerical instabilities for the modified model equation (4). Larger values for M₀ can suppress this numerical instability. While the mobility does not particularly change the particle velocity, we observe that M₀ directly influences the shape of the particle: for small (but still greater than 1) M₀ the particle shows an elliptic form, whereas further increasing M₀ restores a circular shape for the particle. The dependency of the particle shape on the value of M₀ is shown in figure 1(b) for different values of v₀ above the threshold value. Two examples of the morphology are shown in figures 1(c–f) together with the ${\boldsymbol{P}}$ field.

3.2. 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. [9, 36]. 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₀ and activity v₀ are shown in figures 2(a) and (b). The elastic deformation of a single particle during a collision can be seen in figure 2(c), where the eccentricity is plotted as a function of time, whereas a time series of a single collision is shown in figure 2(d). Collisions of deformable particles have also been considered in agent-based models [8], with a qualitatively similar behavior. However, the deformations in our approach strongly depend on M₀ 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₀.

Zoom In Zoom Out Reset image size

All results indicate the particle alignment to be not instantaneous. There is an initial oscillatory phase, whose length and magnitude depends on M₀ and v₀. 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.

3.3. Collective motion in confined geometries

As already analyzed using agent-based, e.g. [7, 8], and phase field models [9, 10] the interaction of a moderate number of particles can lead to collective motion. We here consider simulations with $\simeq 100$ particles to recover these results. To analyze the phenomena we define the translational order parameter ${\phi }_{T}$ and the rotational order parameter ${\phi }_{R}$ as

$$\begin{eqnarray}&&{\phi }_{T}(t)=\displaystyle \frac{1}{n}\left|\displaystyle \sum _{i=1}^{n}{\hat{{\boldsymbol{v}}}}_{i}(t)\right|,\qquad {\phi }_{R}(t)=\displaystyle \frac{1}{n}\displaystyle \sum _{i=1}^{n}{\hat{{\boldsymbol{e}}}}_{{\theta }_{i}(t)}\cdot {\hat{{\boldsymbol{v}}}}_{i}(t)\end{eqnarray} \tag{ 5 }$$

where ${\hat{{\boldsymbol{v}}}}_{i}(t)$ is the unit velocity vector of particle i at time t, ${\hat{{\boldsymbol{e}}}}_{{\theta }_{i}(t)}=(-\sin ({\theta }_{i}(t)),\cos ({\theta }_{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 ${\phi }_{T}=1$ and rotational migration leads to ${\phi }_{R}=\pm 1$.

3.3.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. Figure 3 shows the resulting behavior: at the beginning (figure 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 (figure 3(b)). If these blocks collide they change their direction until all particles in the system are traveling in the same direction (figure 3(c)). This behavior is further confirmed by the translational order parameter ${\phi }_{T}\simeq 1$ after a certain time, see figure 3(d). We also observe that this migration state is not particularly affected by the value of the mobility M₀.

Zoom In Zoom Out Reset image size

3.3.2. Vortex formation and oscillatory motion

We now consider a confined geometry and specify $\psi =0$ and ${\boldsymbol{P}}=0$ at the boundary, which ensures that particles cannot leave the domain. The first geometry is a disk. Here again at the beginning the particles move in a chaotic way, see figure 4(a). Then a vortex is formed and most of the particles follow an anti-clockwise trajectory. In the center some particles move in the opposite direction, see figure 4(b). Eventually also these particles are forced to align with the rest of the system, see figure 4(c). This behavior is confirmed by the rotational order parameter ${\phi }_{R}$. If clockwise or anti-clockwise rotation is observed depends on the initial condition, we therefore plot $| {\phi }_{R}|$ instead of ${\phi }_{R}$ in figure 4(d). Similarly to the collective migration case studied above, the mobility M₀ does not play a major role in the formation of the vortex.

Zoom In Zoom Out Reset image size

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\leqslant 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 figure 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$. Figure 5(c) shows the scaled mean position and mean velocity of the particles together with the computed energy ${{ \mathcal F }}_{{\rm{pfc}}}$. Maxima in the energy thereby correspond to turning point of direction, whereas minima are associated with situations, where the mean particle velocity is constant. If the particles hit the high curvature region the particles get jammed, which results in an increase in ${{ \mathcal F }}_{{\rm{pfc}}}$, resampling the elastic effects. After a certain energy value is reached an energy barrier can be overcome and the system starts to relax by moving in the opposite direction.

Zoom In Zoom Out Reset image size

3.3.3. 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 [37] that allows us to use a direct solver locally. The approach is implemented in AMDiS [38, 39] and shows good scaling properties, which allows the consideration of systems with ≃15,000 particles on the available hardware.

We again consider a system with periodic boundary conditions for which collective migration and cluster formation is expected also for larger numbers of particles. 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 $\phi =n\sigma /| {\rm{\Omega }}|$, with number of particles n, domain size $| {\rm{\Omega }}|$ and σ the area occupied by a single particle, which is equal to $\sigma =\pi {(d/2)}^{2}$, with $d=4\pi /\sqrt{3}$ the lattice constant determined by the free energy equation (2). For $\phi \gt 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 figure 6(b).

Zoom In Zoom Out Reset image size

If we decrease the volume fraction the behavior changes. For very low density, $\phi =0.03$, figure 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, $\phi =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 figure 7(b). Further increasing the density, $\phi =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 figure 7(c). This behavior is similar to the results in [40] for (quasi-) two-dimensional colloidal suspensions of self-propelled particles. For these systems we compute the standard deviation ${\rm{\Delta }}N$ as a function of the mean number of particles N. For active systems it is theoretically predicted that ${\rm{\Delta }}N$ scale as ${N}^{\alpha }$, with $0.5\lt \alpha \leqslant 1$ and giant number fluctuations occur if $\alpha \approx 1$, see [15].

Zoom In Zoom Out Reset image size

We compute α by considering different subregions of our computational domain. The results for ≃600 particles are shown in figure 8, demonstrating an increase of α with increasing ϕ, with the largest value reached being $\alpha =0.79$. This continuous increase in α, as well as the obtained values are consistent with the behavior found in [24] for moderate numbers of particles but larger volume fractions. However, there is a significant difference, the formed clusters in [24] are stationary, our clusters are mobile, similar to the light activated living crystals in [25]. The experimental results as well as the simulations in [25] lead to similar values of α as ours already for smaller ϕ.

Zoom In Zoom Out Reset image size

In summary, we extended the phase field crystal model for active crystals [17, 18], which combines the classical phase field crystal model of Elder et al [26, 27] with a polar order parameter and a self-propulsion velocity, by considering individual active particles. This can be realized by penalizing a locally vanishing one-particle density, as considered in [31]. The resulting microscopic 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 ${\rm{\Delta }}N$ scales as ${N}^{\alpha }$ 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 [25].

The proposed microscopic 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 [41] and [34, 42–44], respectively. Other variants of the phase field crystal model have also been used to simulate the dynamics of epithelial cell colonies [45]. Together with efficient numerical algorithms, see e.g. [37], 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 [46, 47].

This work is funded by the European Union (ERDF) and the Free State of Saxony via the ESF project 100231947 (Young Investigatorts Group Computer Simulations for Materials Design—CoSiMa). We used computing resources provided by JSC within project HDR06.