Active phase field crystal systems with inertial delay and underdamped dynamics

Active matter systems often are well approximated as overdamped, meaning that any inertial momentum is immediately dissipated by the environment. On the other hand, especially for macroscopic systems but also for many mesoscopic ones particle mass can become relevant for the dynamics. For such systems we recently proposed an underdamped continuum model which captures translationally inertial dynamics via two contributions. First, convection and second a damping time scale of inertial motion. In this paper, we ask how both of these features influence the collective behavior compared to overdamped dynamics by studying the example of the active phase field crystal model. We first focus on the case of suppressed convection to study the role of the damping time. We quantify that the relaxation process to the steady collective motion state is considerably prolonged with damping time due to the increasing occurrence of transient groups of circularly moving density peaks. Finally, we illustrate the fully underdamped case with convection. Instead of collective motion of density peaks we then find a coexistence of constant high and low density phases reminiscent of motility-induced phase separation.

We also propose a coarse-grained continuum approach to active systems which aims to incorporate translational inertia to the dynamics [23]. The derived underdamped equations of motion consider inertia via two qualitative contributions. First, inertial convection captures the persistence of mesoscopic particle flows even in the absence Supplementary material in the form of 2 .mp4 files available from the Journal web page at https://doi.org/10.1140/epje/i2020-11971-x a e-mail: michael.schmiedeberg@fau.de of driving forces. Second, the damping time introduces a delay time scale between acting forces and the resulting change in velocity. Also in [23], we apply this model to an example system with competing alignment interactions and find additional dynamical states in the system compared to the overdamped limit case. We trace the occurrence of these states back to the presence of convective flows in the underdamped system. Further, we find that the damping time has no impact on the observed state diagram in the long time limit. This raises the question whether the damping time as a parameter is relevant for the large scale collective dynamics of active systems. We approach this question here by studying another example system. To be specific, we use our model from [23] to extend the overdamped dynamics of the active phase field crystal model by Menzel et al. [24,25] to the underdamped regime. We briefly report on this system with fully underdamped translational dynamics. However, our main focus will be on the approximated case where inertial convection can be neglected or is absent because we want to study the influence of the damping time on the collective dynamics. The original phase field crystal (PFC) model for passive systems [26][27][28] is a continuum approach to crystalization processes. While its density mean field still resolves single particles spatially it averages over the diffusive time scale. Menzel et al. extend this approach to active systems by coupling the system to an additional vector field which captures the local orientation of self propulsion. The characteristic density peaks of the passive PFC model are still present and may still be associated with single particles, therefore allowing to directly analyze particle trajectories where diffusive noise is effectively averaged out. But we note that especially in the context of active systems other interpretations of density peak structures are reasonable as, e.g., periodic patterns of accumulated active particles are predicted to form [29]. Within their active PFC model Menzel et al. find that for low activity density peaks order hexagonally while remaining at rest (resting crystal state). However, above a critical active drive the peaks start to move due to self-propulsion. A local alignment interaction of orientations induces the formation of translationally migrating clusters of density peaks which coarsen over time. In the long time limit global collective motion is observed (traveling crystal state).
The paper is structured as follows: in sect. 2 we outline our previously introduced model for underdamped active matter [23], then apply it to active PFC systems. Furthermore, we motivate and introduce the case of suppressed convection and describe how we analyze our simulation results. The results are discussed in sect. 3 where we first focus on the simplified underdamped system with suppressed convection and then illustrate results of the fully underdamped case. Finally, we conclude in sect. 4.

Active matter with underdamped dynamics
We consider an ensemble of N interacting particles with mass m which self-propel with a constant force f 0 along their individual polar orientation given by a unit vector u i . In two-dimensional systems one angle ϕ i determines the orientationû i = (cos(ϕ i ), sin(ϕ i )) . The corresponding Langevin equations for the underdamped motion of particle i then read where γ = α/m is the damping rate with friction α, and D, D R are translational and rotational diffusion constants.
i (r N ) summarizes an external potential force F (1) (r i ) and interactions F (2) i (r N ) = − j ∇ ri V (2) (r i , r j ) with the pairwise interaction potential V (2) . Analogous to these translational forces particles change their orientation due to the torque (2) (r N , ϕ N ) with a one particle contribution G (1) describing a particle's ability to reorient independent of other particles and a two-body torque G (2) which sums up pairwise interactions.
Based on eqs. (1) we proposed a continuum model for underdamped active matter. For our derivation via Dynamical Density Functional Theory we refer to [23]. The found dynamical equations for the mean fields of number density ρ(r, t), velocity v(r, t) and polar orientational order P(r, t) read with the convective derivative D Dt = ∂ ∂t + (v · ∇). Pairwise interaction forces and torques are now captured via effective free energy functionals F and F P , respectively. Translational diffusion and external potentials also are incorporated in F. The one particle torque G (1) appears again, locally averaged over the distribution of particle momenta and orientations. The self-propulsion velocity v 0 = f 0 /α corresponds to the steady state velocity of a particle in an environment with friction constant α and accelerated by the active force f 0 .
We emphasize again that considering underdamped dynamics introduces two qualitative features in eqs. (2) not present in an analogous overdamped model. First, inertial convection in the velocity and orientation dynamics and second an additional time scale of inertial motion in the velocity dynamics given by the inverse damping rate γ −1 = m/α. We also note that our derived underdamped model eqs. (2) only considers translational inertia, since we assumed in the underlying microscopic eqs. (1) that the orientational dynamics is effectively overdamped. This reflects the situation of particles with their mass closely distributed around their center of mass, meaning that their moment of inertia is negligible. Further, note that in this approach the direction of the actual velocity and therefore of the inertial effects might be different from the direction of the intrinsic self-propulsion.
In [23] we study an example system without any preferred length scale in the density field, but competing alignment interactions at different length scales. For the underdamped dynamics in eqs. (2) we find additional states in this system compared to its overdamped limit case. In this article, we study another example system by choosing the unspecified interaction terms in eqs. (2) according to the original overdamped active PFC model by Menzel et al. [24,25]. The functionals and the mean onebody torque are set to F is given by the PFC functional [26,27]. We observe that the first term of the orientational interaction functional F P is conceptually equivalent to the Frank-Oseen free energy of liquid crystals with one bending constant C 1 [30]. For C 1 > 0 a homogeneous state of uniformly aligned particle orientations minimizes the free energy. Contrarily, the second term in F P , whose strength is given by the activity parameter v 0 , drives an instability from the isotropic state P = 0. We choose its form as in [24] as the simplest nontrivial coupling between the order parameter fields to induce structure formation. Physically, it describes an effective reorientation of particles away from higher densities.
Since the PFC functional only describes variations of a density around its constant bulk mean value, the parameter field ρ may be interpreted as such after rescaling which justifies a constant mobility approximation for the continuity equation in eqs. (2). After applying these simplifications we may now write down the dynamical equations for the active PFC model with underdamped dynamics in dimensionless form by inserting eqs. (3) in eqs. (2) and applying the rescaling rules in appendix A, leading to We will demonstrate in sect. 3.2 that this system, like the one in [23], forms different dynamical states compared to its overdamped limit due to convective flows. But for the main part of this work we instead focus on how the inertial damping time γ −1 alters the overdamped dynamics. Thus, in addition to the fully underdamped model described by eqs. (4) that is briefly studied in sect. 3.2, we simplify the dynamics in the next section to situations where an inertial delay is present but convective effects are absent or suppressed.

Active PFC systems with inertial delay
For the active system that we have studied in [23] we find qualitatively different states when comparing under-and overdamped dynamics. We trace this difference back to the presence of inertial convection in the underdamped model. Furthermore, we find that the damping parameter γ −1 is irrelevant for the observed state diagram in the long time limit.
Here, by studying the active PFC model as an example system, we aim to illustrate how this damping time of individual particles alters the collective dynamics in active systems. To do so, we focus in sect. 3.1 on the situation where convection in the orientation dynamics of eqs. (4) is effectively absent. In principle, we might want to additionally neglect convection of the velocity to assure that any difference to overdamped dynamics is caused by the inertial damping time. However, at least for the active PFC system we find that convection in the velocity alone does not alter the collective behavior. We refer to the case without orientational convection as inertial delay dynam-ics which for the active PFC sytem reads Physically, in the context of underdamped dynamics, this can correspond to the case of low velocities or a slowly varying orientation field, leading to a neglectable impact of convection. Another physical situation are active systems externally controlled via a feedback loop where the orientation is no intrinsic property of moving particles but determined at each fixed point in space [31]. The neglection of the convective term then corresponds to the change of the comoving reference frame to a laboratory frame. Further, when loosening from the background of underdamped motion, the dynamics eqs. (5) also is relevant for strongly damped systems where the velocity dynamics is modified by a memory time which introduces a delay between momentary acting force and thereby induced motion [32][33][34]. In this context, we restrict to delay for translational motion, measured here by the mass parameter, while our orientational dynamics is fully overdamped [35]. We emphasize again that the goal of our study of this model is to understand the influence of the inertial delay effect (depending on the inertial damping time) as a modification to the otherwise overdamped dynamics. Note that the overdamped limit of eqs. (5) is which coincides with the overdamped active PFC system studied by Menzel et al. [24,25] except for the additional prefactor |ρ| in front of the PFC interaction term. This prefactor occurs because we use the same rescaling rules as for eqs. (5) and can be dealt with via further rescaling (see appendix A).
In this dimensionless form of eqs. (4), (5) and (6) there are two independent PFC interaction parameters, the temperature-like ε as well as the mean density variation ρ. We always use ε = −0.98 and ρ = −0.4 such that for the usual passive PFC model we are in the hexagonal crystal phase [26]. Furthermore, we fix D R = 0.1 and C 1 = 0.2 to focus on the parameters associated with the system's dynamics, the active drive v 0 and mass m. All numerical calculations started from homogeneous initial conditions overlaid by small noise.

Methods
To analyze our numerical data we track the position of traveling density peaks like the ones shown in fig. 1 and assign an effective migration velocity to each. Especially for higher m and v 0 values the method in [24,25] used for this becomes unreliable here since then fluctuations in the density field increase and the peak structure is not as clear as in the overdamped limit. Instead, we first locate maxima in the density field and identify to each such found peak the area around it where the density field exceeds a threshold value. By averaging the velocity field v over this area we assign to each peak a discrete velocity of migration v net . As can be seen in fig. 1(a) the velocity field v at a density peak is maximal close to the peak maximum around which it is averaged to v net . Because of this v net overestimates the actual peak migration velocity in the sense of a displacement velocity. But for the following discussion of the results the absolute values of velocities are irrelevant since we only consider relative values. By restricting the area for velocity averaging to a given maximal distance from the corresponding peak position and by only considering peak positions with a minimal distance between them we suppress the influence of fluctuations in the density field.

Active PFC systems with inertial delay
Generally, we find that the active PFC systems with inertial delay and normal overdamped dynamics, eqs. (5) and (6), respectively, predict the same steady states in the long time limit. Hence, we only briefly recap these here and refer to the results for the overdamped model already studied in [24,25] for more details. For low activity v 0 , the system is in the resting crystal state. The kinetic energy input due to the particle's activity melts crystals close to the liquid-solid phase boundary and therefore shifts the latter to lower temperatures by Δ . The hexagonal peak structures remain at rest. Numerics and linear stability analysis of the resting crystal state for eqs. (5) predict Δ ∝ v 2 0 α/C 1 which is in accordance with the findings of Menzel et al. [24]. Above the critical activity v 0,c ≈ 0.3 the phase diagram is not shifted further to lower temperatures. Instead, the additional kinetic energy input from activity is used for directed self-propulsion of the density peaks as depicted in fig. 1. Both, inertial delay and overdamped dynamics, then predict the same traveling crystal state in the long time limit where all hexagonally ordered density peaks have aligned their migration direction and move with the same velocity.
The mass parameter m does not influence the traveling crystal state in the long time limit. Therefore, introducing the inertial damping time γ −1 ∝ m does not qualitatively change the non-equilibrium state. What indeed changes is the duration of the collective relaxation process to this steady state. As we describe in the following, the relaxation time considerably increases with m and cannot be simply identified with the damping or delay time γ −1 as it emerges as a collective effect from the interplay of interactions and dynamics. To explain this emergent behavior we now characterize the transition process from the first stage of randomly oriented density peaks to the final traveling crystal state with one global migration direction of aligned density peaks. As described in [24,25] for the overdamped active PFC model, randomly oriented density peaks start to align locally due to the alignment interaction given by F P . We refer to these regions of collectively moving density peaks as translationally moving clusters. Over time some clusters grow in size on the cost of others with different migration directions until in the long time limit one global orientation emerges. As shown in fig. 1 we observe this alignment of neighboring density peak orientations also within the active PFC model with inertial delay. However, for higher mass parameter m we increasingly observe regions where density peaks migrate circularly around a common center point. We refer to these regions as rotating clusters and provide a visual example thereof in the supplemented time lapse 1 of the density field which is also shown in fig. 2(a). These rotating clusters are transiently stable and delay with their existence the emergence of a global migration direction, i.e. the steady traveling crystal state. To quantify the amount of rotating clusters and their impact on the large scale dynamics we use the concept of circulation Γ . For a continuous velocity field v, it is defined as the closed line integral Γ c (r) = ∂A v · dl over the boundary of a surface A. Via Stokes theorem it can also be understood as the flux of the vorticity field ω = ∇ × v through A and thus measures the amount of circular motion in this area. Similarly, we compute a circulation Γ (r) from the discrete density peak velocities v net by averaging the velocity components in angular direction e l around any given point r over all density peaks within a radius R max around this point. Formally, the circulation then reads where the average runs over all velocities v net of density peaks with relative positions R from r where R < R max .
A typical example of a resulting circulation field is shown in fig. 2(b). We see that at the centers of rotating clusters Γ has a maximum or minimum depending on whether circular motion happens clockwise or counter-clockwise. Thus, we determine the number of rotating clusters at a given time N Γ by counting the number of extrema in the circulation field. To better distinguish rotating from translationally moving clusters, only extrema with values |Γ (r)| > 2 π v net (t) are counted. This threshold is the maximum circulation at the boundary between two translationally moving clusters where density peaks migrate anti-parallel in opposite directions with the momentary sample-averaged density peak velocity amplitude. Furthermore, we define the size of a rotating cluster as the area around an extremum in Γ where |Γ (r)| > 0. An upper bound for the distance between an area element and the extremum position is used to avoid including adjacent regions of dominant translational motion into the area estimation. In fig. 2(b) each rotating cluster is indicated with a circle of radius R Γ reflecting its size.
The presence of transiently stable rotating clusters prolongs the relaxation to the steady traveling crystal state. We quantify the relaxation time scale from the temporal evolution of the area fraction of rotating clusters. In the simulation box with extensions l x and l y the latter reads This number measures how strongly rotating clusters govern the momentary global dynamics. A value of η = π/4 corresponds to the extreme case of a frustration-free square lattice of equally sized clusters where each is surrounded by four contrary rotating ones. The time courses of η in fig. 3 generally show an initial rise, reflecting the formation and growth of rotating clusters over a period τ f , followed by a decay period τ d during which they break-up again. Afterward, the system approaches the steady traveling crystal state where N Γ = 0. Motivated by these observations we suggest a phenomenological expression for the area fraction. Treating τ f and τ d as exponential time scales leads to the expression We extract the mass and activity dependent formation and decay time, τ f and τ d respectively, by fitting this function to the values in fig. 3. η 0 = π/4 stays fixed and a variable time offset is used to account for the varying onsets of crystallization from the supercooled liquid.
The found values for τ f and τ d are shown in the inset of fig. 3. They represent the main result of this article since they summarize the impact of single particle inertial damping time γ −1 ∝ m on the collective dynamics. In order to discuss the dependence on m and activity v 0 we propose an idealized conception for the formation and break-up mechanism of rotating clusters in the following.
After the initialization from homogeneous and isotropic initial conditions small translationally moving clusters form. At the boundary between two contrarily moving ones, density peaks may switch from one of these common migration directions to the other by rotating their orientation of self-propulsion. If peaks from both clusters are doing so, they align with each other too while moving further past each other thereby having to rotate their orientations even further. This disturbance of the boundary between translationally moving clusters eventually drives the emergence of a rotating cluster which grows as more density peaks align to it. When a peak rotates its orientation, the change in motion happens instantly in the overdamped limit dynamics eqs. (6). On the other hand, for the inertial delay dynamics eqs. (5) the actual motion does not follow a change in orientation immediately. Instead of joining the contrarily moving cluster, a density peak might reorient back to its current one. Therefore, the probability of initializing a rotating cluster decreases with mass, meaning that cluster initializations are spread over a longer time, leading to the increasing τ f values in fig. 3 for v 0 = 0.35. Consistent with this, we find for this activity that the maximum number of rotating clusters observed at one time decreases with increasing m. We also find that rotating cluster radii saturate when equally many peaks align to and detach from it. Rotating clusters then further persist over the decay time scale τ d due to the sufficiently given local alignment. However, finally, the systems' preference for a hexagonal density peak structure prevails causing all rotating clusters to break apart. Then, translationally moving clusters with hexagonal symmetry further coarse-grain to the final traveling crystal state. For higher mass values m the translational PFC force becomes less relevant compared to the rotational alignment interaction meaning that the break-up of rotating clusters is delayed. This expresses in the increasing decay times τ d of rotating clusters observed in fig. 3 for v 0 = 0.35. In principle, the same mechanisms of formation and decay apply to the case v 0 = 0.4. However, the higher active drive increases the acceleration along a peak's orientation of self-propulsion. Judging from the higher N Γ values we obtain and their quicker rise during the formation time scale, this increased activity suffices to considerably increase the probability of initiating a rotating cluster. As a consequence, we solely observe low τ f values in fig. 3. Increasing the active drive relative to the strength of the PFC force leads to longer lifetimes of rotating clusters. Thus, τ d initially rises stronger with m for the higher activity. The following saturation of the decay time at high m values is explained with the high maximal area fraction close to η 0 . The corresponding high number of rotating clusters makes mutual destructive interactions inevitable thereby effectively limiting the decay time scale.

Underdamped active PFC systems
So far, we have studied the inertial delay dynamics for the active PFC system eqs. (5) to observe how the time scale of inertial motion γ −1 , as one qualitative feature of underdamped motion, modifies the collective dynamics relative to the overdamped limit. Inertial convection for the polarization has not been considered so far. Therefore, we now briefly report on the case of general underdamped dynamics which includes convection for the active PFC system eqs. (4).
For the inertial delay dynamics discussed so far we observe the same state diagram as for the normal overdamped case. Especially the same traveling crystal state is observed in both cases where in the long time limit density peaks order hexagonally and collectively migrate by globally aligning their directions of self-propulsion. In contrast, in the same parameter regime the underdamped active PFC system instead shows a more complex phase coexistence of constant high and low density regions which are dynamically reorganized, see fig. 4(a) and the corresponding time lapse (see footnote 1 ). The velocity and orientation field in these regions are isotropic. The constant high and low density values are a result of density accumulation or dispersion through active particle transport and the limiting cubic nonlinearity of the density in the PFC interaction. Active particle transport solely happens in regions of intermediate density where activity induces directed fluxes within structures of hexagonally ordered density peaks. By comparing the velocity and orientation fields in fig. 4(b), (c) with those in fig. 1 it can be seen that the inclusion of convection alters the migration mechanism of these density peaks. For inertial delay dynamics orientation and resulting velocity are induced close to the density peak maximum such that peaks migrate in this direction. The opposite happens for underdamped dynamics where convection modifies this steady state such that orientation and velocity are most expressed between density peaks and the associated particle flux is directed opposite to the motion of density peaks, which is now comparable to a counterpropagating wave within a traffic jam (cf. supplied time lapse (see footnote 1 ) where a constant high density region grows due to particle influx of surrounding density peak structures while peaks move away from the high density region). The observed phase separation into high and low density regions is a defining feature of motility-induced phase separation [36] whose occurrence depends on the particle mass in case of underdamped dynamics [22]. Whether such a dependence is present in the active PFC system or in a modification is an interesting question for future works.

Conclusion
In this article, we have studied an active PFC system with inertial delay dynamics for translational motion. The distinction of translationally moving and rotating density peak clusters is used to quantify the time scale of the collective relaxation process to the steady state. The suggested formation and break-up mechanism for rotating clusters qualitatively explain the observed parameter dependence. Generally, the relaxation process to the steady traveling crystal state is prolonged for increasing damping parameter γ −1 ∝ m. However, we emphasize that the duration of this process does not trivially scale with γ −1 which only captures the relaxation of individual particle velocities. Instead, we find that formation and decay times of rotating clusters capture the time scale of the transient collective dynamics. For the case of identifying the characteristic PFC density peaks with single particles, we note that underdamped particle based simulations with an effective alignment interaction also find transient rotating particle clusters before global collective motion emerges [37]. The studied inertial delay dynamics are relevant for, e.g., active systems with memory. In our case the damping time γ −1 can be seen as an effective delay time for the translational dynamics [35]. Other studies with delay for the angular dynamics also find an impact of the delay parameter on collective phenomena [33,34].
In conclusion, while convection may alter the long-time steady state of active systems with inertia, our results show in addition that the memory due to inertial delay induces transient phenomena that persist for time scales that are much longer than the inertial damping time.
Concerning an active PFC system with underdamped translational dynamics we find, as in another example system that we have studied previously in more detail [23], that inertial convection alters the collective dynamics relative to the overdamped limit, leading to different states in the long time limit. To be specific, instead of the traveling crystal state we find a coexistence of high and low density phases which are continually reorganized by directed particle fluxes which are induced by activity and modified by an interplay of structure forming interactions and inertial convection. This example demonstrates the relevance of our underdamped active matter model [23] for future works in the currently opening realm of underdamped active matter research [4]. Specifically, the emergence of motility induced phase separation in active systems is predicted to depend on the time scale of inertial motion, while high and low density phases have different kinetic temperatures [22]. Extending our underdamped model eqs. (2) to varying kinetic temperatures [38] might lead to a promising continuum approach for these systems.

Author contribution statement
MS conceptualized the project and guided the research. DA derived the model, carried out the simulations, and analyzed the results. Both authors wrote the manuscript.

Appendix A. Rescaling
A variable x * is rescaled to its dimensionless form x via the factor X as x * → X x. For simplicity of notation we use the same symbol for dimensional and rescaled variables. The shown dimensionless form of the model eqs. (4), (5), and (6) is obtained by applying the rescaling rules with T = αu 1/2 (λq 4 0 ) 3/2 q 2 0 . For convenience, we additionally rescale the orientation field |ρ|P → P in eqs. (6) with the dimensionless parameter ρ which is always fixed to the same value. Thus, we also might have used ρ as a numerical constant and write T /|ρ| instead of T in eqs. (A.1) to arrive at a form of the overdamped eqs. (6) equivalent to the form in [24,25]. However, this last rescaling is not done here.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.