Jerky active matter: a phase field crystal model with translational and orientational memory

Most field theories for active matter systems are derived in the overdamped limit, despite the fact that it was noted at an early stage of their development that incorporating inertia can lead to interesting modifications of the dynamics [1]. However, there is a current increase of interest in the role of memory and (inertial) delay in active matter [2–11]. Experiments have shown inertial [2–4] and sensorial [5, 6] delay to be important for the dynamics of self-propelled particles. Examples for such experiments include vibrobots driven by an electromagnetic shaker, where an inertial delay between orientation and velocity of the particles was observed [2], camphor surfers driven by surface tension gradients and hexbug crawlers driven by a vibrating motor, both exhibiting underdamped active motion [12], and phototactic robots responding to measured light intensities with a certain (sensorial) delay [5]. Inertia can lead to interesting effects that are not present in overdamped active systems [2–4, 13–16], such that inertial active matter models, which range from particle-based descriptions such as the underdamped active Ornstein–Uhlenbeck particle model [18] or the underdamped active Brownian particle model [2, 4, 10], to coarse-grained field theories [15, 17], form an important extension of usual models. A particularly interesting aspect that remains to be explored is the fact that active particles have translational and rotational degrees of freedom which both are associated with inertia, but might relax on different time scales [4, 19]. Memory and inertia effects are intimately connected and sometimes equivalent, as can be seen from an analysis in the Mori–Zwanzig framework (see section 2).

The problem of memory effects has been discussed in other contexts, very notably in particle physics [20–22]. Most dissipative transport equations (such as the diffusion equation) are acausal, since they assume that signals propagate with an infinite velocity. This is typically justified by arguing that microscopic relaxation processes occur infinitely rapidly compared to the relevant dynamics. However, this assumption is not consistent with special relativity where the maximum propagation speed of signals is the speed of light. Relativistic descriptions are required for phase transitions in heavy-ion collisions, where the relevant dynamics itself is very fast, such that neglecting memory effects and causality constraints is an approximation that cannot be justified [21]. For this reason, more general theories have been derived that take these effects into account [20–23].

This problem can also arise in soft matter systems such as polymer solutions [23]. Although relativistic effects are typically irrelevant here, it is still an approximation to assume that the macroscopic order parameters change much slower than the microscopic degrees of freedom. The time scale will then be set by the damping coefficients rather than by the speed of light. As noted by Archer [24, 25], who suggested this as a topic of further investigation within classical dynamical density functional theory (DDFT) [26], the mathematical structure of transport equations for relativistic heavy-ion collisions and underdamped soft matter is identical. There is, however, an additional difficulty that can occur in the latter case, namely that two different relaxational time scales are relevant. A typical example would be soft matter systems consisting of particles with orientational degrees of freedom—in particular active matter—where position and orientation relax on different time scales. However, different time scales can also be relevant for other systems. Here, we discuss this problem in a very general way and then specialize our results to the case of active field theories.

Phase field crystal (PFC) models [27], which provide a coarse-grained description of crystallization and pattern formation in materials, are a particularly suitable framework for such investigations. They were first proposed phenomenologically [28–30] and then derived from (dynamical) density functional theory [31, 32]. The early forms have been extended into a variety of directions including orientational degrees of freedom [33–36]. A particularly important variant is the active PFC model, which was developed by Menzel and Löwen [1]. The active PFC model has been investigated and extended by a variety of researchers, including Menzel et al [37], Alaimo et al [38, 39], Chervanyov et al [40], Praetorius et al [41], Ophaus et al [42–44], Arold and Schmiedeberg [15, 17], and Holl et al [45]. Moreover, PFC models have been extended to incorporate inertia. An important step in this direction was the 'modified PFC' (MPFC) model, developed by Stefanovic et al [46], which has been investigated by a variety of authors [46, 47, 49, 51]. It can also be obtained by introducing an exponential memory function into the standard PFC model [52] (see section 2). The MPFC model is a second-order ¹ model that allows for propagating sound waves. It aims to overcome a drawback of the original PFC model, namely that it only describes the dynamics on diffusive time scales and therefore does not incorporate fast elastic interactions, by adding an additional time derivative [46]. A microscopic theory involving also phonon dynamics was developed by Majaniemi and Grant [48], this model allows to recover the MPFC model as a limiting case [48, 53]. Archer [24] has derived an underdamped DDFT for atomic liquids that is similar in form to the MPFC model. However, the MPFC approach can lead to inaccuracies in the long-wavelength limit [48, 50, 54]. Another way of extending PFC models (or, more generally, DDFT) toward the underdamped case is coupling to a velocity field [25, 50, 55–58]. A difficulty that arises in such approaches is that the density fields of PFC models vary on microscopic length scales, whereas hydrodynamics considers macroscopic fields [50, 57]. Using the amplitude expansion method [59, 60], elastic equilibrium has been incorporated into PFC models [50, 61]. Possible inaccuracies of incorporating elastic equilibrium via constraints imposed on the phase evolution for defected crystals far from the melting point were pointed out in reference [62]. The time scale separation between diffusive and elastic distortions can be incorporated by adding an additional variable for static elastic distortions [62] (see also references [54, 63]). Consequently, PFC models are a suitable framework for incorporating delay into active field theories, as recently shown by Arold and Schmiedeberg [17]. Inertia has been incorporated into active PFC models by coupling them to velocity fields describing either a background fluid [1] or the momentum density of the active particles themselves [15, 17]. All these models assume the rotational dynamics to be overdamped. PFC models are reviewed in reference [27], their relation to DDFT is also discussed in references [26, 64].

While these models are of second order in temporal derivatives, memory effects have also been linked to third-order models. A typical example is jerky dynamics, which provides differential equations for the rate of change of the acceleration. These can be relevant in Newtonian mechanics if the force is memory-dependent. The resulting third-order differential equations are useful for obtaining simple models of chaos [65–68]. Moreover, third-order dispersions have been studied in optics [69, 70], where they naturally appear in delay systems [71, 72].

In this work, we discuss how to incorporate finite relaxation times in systems where two different time scales are present. This is done in a very general way, since this problem can occur in a variety of contexts. We then specialize to the case of active PFC models and obtain a generalization of the traditional active PFC model [1] that incorporates translational and rotational memory. Notably, the resulting theory is of third order in temporal derivatives and can thus be viewed as a spatiotemporal jerky dynamics. Field theories of this form are almost completely unexplored, and can—like their aforementioned counterparts from other areas of nonlinear dynamics—be expected to show very interesting behaviour. As shown in the present work, this is relevant for experiments on active matter. In contrast to the passive case, it is found that the damping can affect the linear stability of the fluid state in the PFC model.

An interesting application of the new PFC model is the investigation of sound waves in active matter. In passive systems, sound propagates due to translational inertia. However, it has been found that propagating density waves are also possible in overdamped active systems as a consequence of the coupling between density and polarization [73, 74]. Our model incorporates both translational inertia and activity. It therefore allows to describe two different mechanisms of sound propagation that are relevant for active matter. In particular, we show that at higher frequencies, the usual 'passive' sound can experience a stronger damping than 'active' sound waves, such that waves of the latter type can propagate further into the medium.

This article is structured into eight sections. In section 2, we explain how finite relaxation times can be incorporated into transport equations. We provide a brief introduction to jerky dynamics in section 3. A general model for systems with two time scales is developed in section 4. In section 5, we use these results to obtain an active PFC model with memory. We perform a linear stability analysis of this model in section 6. Sound propagation is studied in section 7. We conclude in section 8.

For a conserved order-parameter field $\phi \left( \overrightarrow {r},t\right)$, where $\overrightarrow {r}$ and t denote position and time, respectively, the time evolution is given by the continuity equation

$$\begin{equation}{\partial }_{t}\phi \left( \overrightarrow {r},t\right)=- \overrightarrow {\nabla }\cdot \overrightarrow {J}\left( \overrightarrow {r},t\right)\end{equation} \tag{ 1 }$$

with the current $\overrightarrow {J}\left( \overrightarrow {r},t\right)$. A typical example is standard DDFT [26, 75–78], where the current is

$$\begin{equation} \overrightarrow {J}\left( \overrightarrow {r},t\right)=-\tilde {M}\phi \left( \overrightarrow {r},t\right) \overrightarrow {\nabla }\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}\end{equation} \tag{ 2 }$$

with the mobility $\tilde {M}$ (which for standard DDFT is given by $\tilde {M}=\beta D$ with the thermodynamic beta β and the diffusion constant D), the density ϕ (which is the order-parameter field), and the free-energy functional F. Another example is the Cahn–Hilliard equation [79, 80], where the current is

$$\begin{equation} \overrightarrow {J}\left( \overrightarrow {r},t\right)=-\tilde {M} \overrightarrow {\nabla }\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}.\end{equation} \tag{ 3 }$$

A third example is active model B [81], where the current is

$$\begin{equation} \overrightarrow {J}\left( \overrightarrow {r},t\right)=-\tilde {M} \overrightarrow {\nabla }\left(\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}+\lambda {\left( \overrightarrow {\nabla }\phi \left( \overrightarrow {r},t\right)\right)}^{2}\right)+ \overrightarrow {{\Lambda}}\left( \overrightarrow {r},t\right)\end{equation} \tag{ 4 }$$

with the activity parameter λ and the noise $\overrightarrow {{\Lambda}}\left( \overrightarrow {r},t\right)$. All these currents share a common property: They are Markovian, i.e., they only depend on the value of ϕ at time t, but not on values of ϕ at times s < t. Moreover, they are local, i.e., the time derivative of ϕ at position $\overrightarrow {r}$ does not depend on the value of ϕ at a different position ${ \overrightarrow {r}}^{\;\prime }$.

This is not the most general case. A good starting point for a systematic analysis is the Mori–Zwanzig formalism [82–87]. It allows to derive macroscopic transport equations for an arbitrary set of relevant variables from the microscopic dynamics and thereby, as discussed in reference [87], gives very general insights into the general structure of these equations. As exploited in recent work by Meyer et al [88–91], the formalism is also useful for analyzing memory effects in a systematic way. In general, ignoring some degrees of freedom of a system leads to a transport equation in which the dynamics depends on the state of the system at previous times. Approximating this equation by a Markovian equation is possible if the set of relevant variables one has chosen captures the complete macroscopic dynamics [78, 87].

Let us assume that we know the microscopic equations governing a passive many-particle system, but are only interested in a conserved scalar order parameter $\phi \left( \overrightarrow {r},t\right)$ that is a function of the microscopic degrees of freedom. Then, we write the general time evolution as

$$\begin{equation}{\partial }_{t}\phi \left( \overrightarrow {r},t\right)= \overrightarrow {\nabla }\cdot \left({\int }_{0}^{t}\mathrm{d}s\enspace \int {\mathrm{d}}^{3}{r}^{\prime }\enspace \enspace \mathcal{M}\left(t,s, \overrightarrow {r},{ \overrightarrow {r}}^{\;\prime }\right){ \overrightarrow {\nabla }}^{\prime }\frac{\delta F\left[\phi \right]}{\delta \phi \left({ \overrightarrow {r}}^{\;\prime },s\right)}\right),\end{equation} \tag{ 5 }$$

where $\mathcal{M}$ is the memory kernel. We have ignored here the organized drift term, which is reversible and does not contribute to the change of the entropy [92]. This term typically vanishes for symmetry reasons if there is only one relevant variable—for example, in the derivation of DDFT from the Mori–Zwanzig formalism [78], the organized drift can be dropped since it corresponds to a phase-space integral over a function that is odd in the particle momenta. Moreover, we have dropped a noise term. The absence of noise terms in generalized transport equations for the mean values of the relevant variables can be linked to the assumption that the system's time evolution starts in a state of constrained equilibrium (see equation (95) in reference [26]). Apart from this, the form (5) is completely general.

The rather complicated general form (5) can be simplified significantly if the relevant variable ϕ is slow. In this case, it can be assumed that ϕ is constant on the time scales on which the microscopic degrees of freedom relax. For this to be possible, the relaxation of these microscopic degrees of freedom has to occur very rapidly. If we also assume that the nonlocality can be neglected (e.g., because the system is dilute [78]), the exact transport equation (5) can be approximately written as

$$\begin{equation}{\partial }_{t}\phi \left( \overrightarrow {r},t\right)= \overrightarrow {\nabla }\cdot \left(\mathcal{D}\left( \overrightarrow {r},t\right) \overrightarrow {\nabla }\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}\right)\end{equation} \tag{ 6 }$$

with the diffusion tensor $\mathcal{D}$. This is the so-called Markovian approximation.

Although it is made in almost all practical cases, the Markovian approximation is not innocent. From a foundational perspective, it introduces the thermodynamic irreversibility not present in the time-reversal-invariant microscopic laws of physics [93]. In the Markovian limit, an H-theorem corresponding to an increase of entropy can be proven [94, 95]. From a more practical point of view, the Markovian approximation corresponds to the assumption that the order parameter $\phi \left( \overrightarrow {r},t\right)$ gives a complete description of the macroscopic state. For example, if $\phi \left( \overrightarrow {r},t\right)$ is the number density, making the Markovian approximation implies that the momentum density $\overrightarrow {g}\left( \overrightarrow {r},t\right)$ relaxes very rapidly, i.e., that we are working in the overdamped limit [78].

What is also relevant here is that equation (6) is a diffusive equation which is of first order in temporal and second order in spatial derivatives. This leads, in general, to an infinite propagation speed of signals, i.e., to an acausal equation. Early treatments of this problem include inertial extensions of the heat equation [96] (see reference [23] for a review). This issue is particularly relevant in the case of relativistic systems. If we wish to describe phase transitions in the early Universe or in heavy-ion collisions, we generally have to take into account that signals can only propagate with a finite velocity (namely the speed of light) [21]. Thus, when applying the Mori–Zwanzig formalism to relativistic systems, we generally have to be very careful in handling memory effects [20]. However, as we shall see below, the same problem can arise in soft matter physics, even though the velocities and time scales are very different there.

We start discussing this issue by presenting the theory derived by Koide et al [21] for describing phase-separation processes in relativistic high energy physics. For a conserved order parameter, the Cahn–Hilliard equation is a very successful theory for phase separation. In order to take into account the finite propagation speed in relativistic systems, the Cahn–Hilliard current (3) is modified as

$$\begin{equation} \overrightarrow {J}\left( \overrightarrow {r},t\right)=-M{\int }_{0}^{t}\mathrm{d}s\enspace \enspace {\text{e}}^{-\gamma \left(t-s\right)}\enspace \overrightarrow {\nabla }\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},s\right)}\end{equation} \tag{ 7 }$$

with the damping coefficient (inverse relaxation time) γ and the modified mobility $M=\gamma \tilde {M}$. We have here inserted as the memory kernel a memory function γ exp(−γ(t − s)), which can be motivated by certain assumptions about the noise. For example, the motion of a noninteracting Brownian particle subject to white noise is governed by the usual acausal diffusion equation, whereas the assumption of coloured noise with an inverse correlation time γ gives a causal diffusion equation similar to the models considered here [21]. Such causal diffusion equations are used in the theory of heavy-ion collisions [97]. Moreover, exponential memory kernels have been linked to hyperbolic Master equations. There, additional time scales can require models also of third or higher order [98]. In the Mori–Zwanzig formalism, the memory kernel is related to the correlation of the noise [86].

The time derivative of equation (7) is given by

$$\begin{equation}\begin{aligned}{\partial }_{t} \overrightarrow {J}\left( \overrightarrow {r},t\right)& =-M \overrightarrow {\nabla }\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}+\gamma M{\int }_{0}^{t}\mathrm{d}s\enspace \enspace {\text{e}}^{-\gamma \left(t-s\right)}\enspace \overrightarrow {\nabla }\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},s\right)}\\ & =-M \overrightarrow {\nabla }\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}-\gamma \overrightarrow {J}\left( \overrightarrow {r},t\right).\end{aligned}\end{equation} \tag{ 8 }$$

Equation (8) shows that the assumption of an exponential relaxation of the memory kernel made in equation (7) corresponds to the assumption that the force that leads to the relaxation of the current (the second term on the right-hand side of equation (8)) is proportional to the current. If we assume that $\overrightarrow {J}$ is a momentum or velocity field, we can compare equation (8) to the underdamped DDFT derived by Archer [25] and see that such a term arises from a microscopic derivation under the assumption of a linear friction force (such as Stokes drag) acting on the individual particles. This indicates that a more complex relaxation process (as would arise, e.g., from a nonlinear Stokes force) would require a more general memory kernel.

Differentiating equation (1) with respect to time and inserting equation (8) then gives

$$\begin{equation}{\partial }_{t}^{2}\phi \left( \overrightarrow {r},t\right)+\gamma {\partial }_{t}\phi \left( \overrightarrow {r},t\right)=M{ \overrightarrow {\nabla }}^{2}\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}.\end{equation} \tag{ 9 }$$

This is a second-order causal Cahn–Hilliard equation that is similar in form to the telegrapher's equation. Theories of this form can be used to study hyperbolic spinodal decomposition [99–102]. There are two important limiting cases of equation (9) one can consider. The first one is the overdamped limit γ → ∞ with M/γ = const. (corresponding to fixed $\tilde {M}$), in which case we recover the standard Cahn–Hilliard equation

$$\begin{equation}{\partial }_{t}\phi \left( \overrightarrow {r},t\right)=\tilde {M}{ \overrightarrow {\nabla }}^{2}\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}.\end{equation} \tag{ 10 }$$

The second one is the underdamped limit γ → 0 at fixed M, which gives the inertial Cahn–Hilliard equation

$$\begin{equation}{\partial }_{t}^{2}\phi \left( \overrightarrow {r},t\right)=M{ \overrightarrow {\nabla }}^{2}\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}.\end{equation} \tag{ 11 }$$

For the DDFT current (2), Archer [24, 25] has noted that this procedure leads to a causal DDFT that has the same structure as the DDFT

$$\begin{equation}{\partial }_{t}^{2}\phi \left( \overrightarrow {r},t\right)+\gamma {\partial }_{t}\phi \left( \overrightarrow {r},t\right)=\frac{1}{m} \overrightarrow {\nabla }\cdot \left(\phi \left( \overrightarrow {r},t\right) \overrightarrow {\nabla }\frac{\delta F\left[\phi \right]}{\delta \phi \left( \overrightarrow {r},t\right)}\right)\end{equation} \tag{ 12 }$$

with particle mass m for particles with inertia. A similar result was obtained by Chavanis [103]. Within our above considerations about the Mori–Zwanzig formalism, we can give a physical explanation: for systems where inertia is relevant, such as atomic fluids, it is no longer possible to assume that the number density is the only degree of freedom that is relevant, such that we also require the momentum density. Hence, if we want to derive a transport equation for a system with inertia in which the number density is the only relevant variable (i.e., a DDFT for an atomic fluid), we can no longer make the approximation of infinitely fast relaxations. Instead, we need to take into account that in an underdamped system the velocities need a finite time to relax. From equation (12), we can also see why the over- and underdamped limits have to be taken in the form (10) and (11), respectively: the overdamped limit of equation (12) then leads to standard DDFT [25], whereas the underdamped limit leads to a generalized Euler equation (without convective term). From equation (12), we find that M = 1/m such that $\tilde {M}=1/\left(\gamma m\right)$. Since we know from DDFT that $\tilde {M}=\beta D$, we can infer γ = 1/(βmD). Therefore, varying γ at fixed M corresponds to changing ² D, whereas varying γ at fixed $\tilde {M}$ corresponds to changing m. Recent simulations by Caprini and Marini Bettolo Marconi [104] indicate that the difference between varying the particle mass at fixed drag coefficient and varying the drag coefficient at fixed particle mass is particularly important in active systems. (The drag coefficient is inversely proportional to the diffusion constant [104].)

Another way to look at this issue is to take into account that equations of motion containing memory effects only arise if we do not consider all degrees of freedom of a system, but only a reduced set [87] (Hamilton's equations or the Heisenberg equation of motion have no memory). The memory terms then incorporate (along with the noise) the dynamics of those parts of the system that we do not wish to model explicitly. This also implies that we can, instead of considering the memory, obtain a Markovian dynamics if we enlarge our set of relevant variables [26, 85]. In this case, we can, rather than using a transport equation for the density that contains memory, also write down coupled memoryless equations for mass and momentum density. Mathematically, this is reflected by the fact that the equations with memory are of second order. We could alternatively obtain two equations of motion with the usual first-order structure by using ∂_t ϕ as an additional relevant variable. The reason why ∂_t ϕ appears is that it is no longer negligible if memory terms are relevant. In the example presented here, we would have ${\partial }_{t}\phi =- \overrightarrow {\nabla }\cdot \overrightarrow {g}$, such that using ∂_t ϕ as a relevant variable is equivalent to adding the momentum density $\overrightarrow {g}$ to the set of relevant variables.

For analyzing the problem at hand, the theory of jerky dynamics [66–68, 105–107] will prove to be very useful. A jerky dynamics is an ordinary differential equation of the form [68]

$$\begin{equation} \dddot {x}=J\left(x,\dot {x},\ddot {x}\right)\end{equation} \tag{ 13 }$$

with a time-dependent variable x(t) and a function J. Jerky dynamics is very important for chaos theory [66, 105]. Mechanically, the 'jerk' is the rate of change of acceleration [108]. Intuitively, one might expect jerks to be of no importance in classical mechanics, since Newton's equation of motion

$$\begin{equation}\ddot {x}=\frac{1}{m}F\left(x,\dot {x}\right)\end{equation} \tag{ 14 }$$

with the force F is of second order. (Although a jerky dynamics can obviously be obtained by taking the time derivative of equation (14), it would not provide any physical insights.) However, if the force has memory, i.e., if it depends on values of x or $\dot {x}$ at previous times, taking the time derivative of equation (14) can lead to an interesting jerky dynamics that contains additive terms depending solely on x. Such models can allow for chaos [68].

This relation to memory is what makes jerky dynamics relevant for the present investigation. What will also be useful is that, in some cases, a three-dimensional dynamical system

$$\begin{equation}\dot {x}={f}_{x}\left(x,y,z\right),\end{equation} \tag{ 15 }$$

$$\begin{equation}\dot {y}={f}_{y}\left(x,y,z\right),\end{equation} \tag{ 16 }$$

$$\begin{equation}\dot {z}={f}_{z}\left(x,y,z\right)\end{equation} \tag{ 17 }$$

with dynamical variables x, y, and z and functions f_x , f_y , and f_z can be written in the form (13). The problem is addressed in detail in reference [65], where the conditions under which such a transformation is possible are discussed. A useful strategy [65] is to calculate the first time derivative of equations (16) and (17) and the first and second time derivatives of equation (15). This gives seven coupled equations for $x,\dot {x},\ddot {x}, \dddot {x},y,\dot {y},\ddot {y},z,\dot {z}$, and $\ddot {z}$. These can be used to eliminate $y,\dot {y},\ddot {y},z,\dot {z}$, and $\ddot {z}$, which gives a closed equation for $\dddot {x}$ that solely depends on x, $\dot {x}$, and $\ddot {x}$.

In systems of active (or nonspherical) particles, one needs to take into account both translational and orientational degrees of freedom. As a toy model for describing the physics of such systems, we use the dynamical equation

$$\begin{equation}{\partial }_{t}\phi \left( \overrightarrow {r},t\right)=\tilde {T}\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},t\right)\right\}\right)+\tilde {R}\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},t\right)\right\}\right)\end{equation} \tag{ 18 }$$

with $k\in {\mathbb{N}}_{0}$. Later, we interpret $\tilde {T}$ and $\tilde {R}$ as the contributions from translational and orientational degrees of freedom, respectively. However, the considerations in this section apply to any physical system whose dynamics can be written in the form (18). We are assuming nothing about the form of $\tilde {T}$ and $\tilde {R}$ here—it can be conserved or nonconserved, active or passive, and it can depend on ϕ as well as on spatial derivatives of ϕ. The field ϕ can be a scalar, a vector, or a tensor. As in section 2, we generalize equation (18) toward a causal dynamics with time delay:

$$\begin{align}\begin{aligned}{\partial }_{t}\phi \left( \overrightarrow {r},t\right)& ={\int }_{0}^{t}\enspace \mathrm{d}s\enspace \left({\text{e}}^{-{\gamma }_{T}\left(t-s\right)}\enspace T\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},s\right)\right\}\right)+{\text{e}}^{-{\gamma }_{R}\left(t-s\right)}\enspace R\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},s\right)\right\}\right)\right)\end{aligned}\end{align} \tag{ 19 }$$

with $T={\gamma }_{T}\tilde {T}$ and $R={\gamma }_{R}\tilde {R}$. There is a very important difference to the models known from the literature which we have discussed in section 2: in general, we cannot assume that translational and orientational degrees of freedom relax on the same time scale. Therefore, we have introduced two different relaxation parameters γ_T and γ_R for translation and orientation, respectively. This has important consequences for the resulting dynamics. To see this, we calculate the time derivative of equation (19), which gives

$$\begin{align}\begin{aligned}{\partial }_{t}^{2}\phi \left( \overrightarrow {r},t\right)& =-{\int }_{0}^{t}\enspace \mathrm{d}s\enspace \left({\gamma }_{T}\enspace {\text{e}}^{-{\gamma }_{T}\left(t-s\right)}\enspace T\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},s\right)\right\}\right)+{\gamma }_{R}\enspace {\text{e}}^{-{\gamma }_{R}\left(t-s\right)}\enspace R\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},s\right)\right\}\right)\right)\\ & \quad +T\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},t\right)\right\}\right)+R\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},t\right)\right\}\right).\end{aligned}\end{align} \tag{ 20 }$$

For γ_T = γ_R = γ, we could write equation (20) as

$$\begin{equation}\begin{aligned}{\partial }_{t}^{2}\phi \left( \overrightarrow {r},t\right)=& -\gamma {\partial }_{t}\phi \left( \overrightarrow {r},t\right)+T\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},t\right)\right\}\right)+R\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},t\right)\right\}\right).\end{aligned}\end{equation} \tag{ 21 }$$

However, for γ_T ≠ γ_R the sum of the first two terms in equation (20) is not proportional to ${\partial }_{t}\phi \left( \overrightarrow {r},t\right)$, since the translational and rotational contributions appear with different prefactors due to the different time scales on which the contributions change. Therefore, equation (20) cannot be written as a second-order partial differential equation without time convolution.

Physically, this is due to the fact that we require three rather than two variables for a complete description of the system—in addition to mass and momentum density the angular momentum density is needed. Therefore, we require a third-order partial differential equation to describe the dynamics of $\phi \left( \overrightarrow {r},t\right)$ without time convolution. This can be obtained using the procedure introduced in section 3.

First, we make the definitions

$$\begin{align}x\left( \overrightarrow {r},t\right)& =\phi \left( \overrightarrow {r},t\right),\end{align} \tag{ 22 }$$

$$\begin{align}y\left( \overrightarrow {r},t\right)& ={\int }_{0}^{t}\enspace \mathrm{d}s\enspace \enspace {\text{e}}^{-{\gamma }_{T}\left(t-s\right)}\enspace T\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},s\right)\right\}\right),\end{align} \tag{ 23 }$$

$$\begin{align}z\left( \overrightarrow {r},t\right)& ={\int }_{0}^{t}\enspace \mathrm{d}s\enspace \enspace {\text{e}}^{-{\gamma }_{R}\left(t-s\right)}\enspace R\left(\left\{{ \overrightarrow {\nabla }}^{k}\phi \left( \overrightarrow {r},s\right)\right\}\right).\end{align} \tag{ 24 }$$

For the remainder of this section, we drop the dependence on $\overrightarrow {r}$ since it is not important for the further calculations. With the definitions (22)–(24), we can obtain from equation (19) the dynamical system

$$\begin{align}\dot {x}& =y+z,\end{align} \tag{ 25 }$$

$$\begin{align}\dot {y}& =-{\gamma }_{T}y+T,\end{align} \tag{ 26 }$$

$$\begin{align}\dot {z}& =-{\gamma }_{R}z+R.\end{align} \tag{ 27 }$$

The problem has thus been reduced to the derivation of the jerky dynamics corresponding to the dynamical system given by equations (25)–(27). This can be done using the procedure discussed in section 3, namely, by computing $\ddot {x}$, $\ddot {y}$, $\ddot {z}$, and $\dddot {x}$, and solving the resulting systems of equations in order to express $\dddot {x}$ as a function of x, $\dot {x}$, and $\ddot {x}$. The derivation is shown in detail in appendix A and results in

$$\begin{equation} \dddot {x}=-\left({\gamma }_{T}+{\gamma }_{R}\right)\ddot {x}-{\gamma }_{T}{\gamma }_{R}\dot {x}+{\gamma }_{R}T+{\gamma }_{T}R+\dot {T}+\dot {R}.\end{equation} \tag{ 28 }$$

PFC models are a useful framework for the description of soft and active matter [1, 27, 42, 43, 45, 109]. They can be obtained as a limiting case of the more complex and more general case of DDFT [32, 37, 64]. Reviews are given by references [27] (PFC models) and [26] (DDFT), both reviews discuss the derivation of PFC models from DDFT. The order parameters are the rescaled density $\psi \left( \overrightarrow {r},t\right)$ and the polarization $\overrightarrow {P}\left( \overrightarrow {r},t\right)$, which arise through an orientational expansion [110, 111] of the one-body density $\rho \left( \overrightarrow {r},\hat{u},t\right)$ depending on position $\overrightarrow {r}$ and orientation $\hat{u}$. This expansion has the form

$$\begin{equation}\rho \left( \overrightarrow {r},\hat{u},t\right)=\bar{\rho }\left(1+\psi \left( \overrightarrow {r},t\right)+ \overrightarrow {P}\left( \overrightarrow {r},t\right)\cdot \hat{u}+\cdots \enspace \right),\end{equation} \tag{ 29 }$$

with the constant reference density $\bar{\rho }$, the rescaled density

$$\begin{equation}\psi \left( \overrightarrow {r},t\right)=\frac{1}{2\pi \bar{\rho }}\left({\int }_{0}^{2\pi }\enspace \mathrm{d}\varphi \enspace \enspace \rho \left( \overrightarrow {r},\hat{u},t\right)\right)-1,\end{equation} \tag{ 30 }$$

and the polarization

$$\begin{equation} \overrightarrow {P}\left( \overrightarrow {r},t\right)=\frac{1}{\pi \bar{\rho }}{\int }_{0}^{2\pi }\enspace \mathrm{d}\varphi \enspace \enspace \rho \left( \overrightarrow {r},\hat{u},t\right)\hat{u}\end{equation} \tag{ 31 }$$

with the polar angle φ determining the orientation $\hat{u}\left(\varphi \right)={\left(\mathrm{cos}\left(\varphi \right),\mathrm{sin}\left(\varphi \right)\right)}^{\mathrm{T}}$. The definitions (30) and (31) correspond to a two-dimensional system, the three-dimensional case (see reference [110]) is similar. An expansion of the form (29) is equivalent to a Fourier expansion in two and to an expansion in spherical harmonics in three dimensions [110]. One can derive the active PFC model by inserting the parametrization (29) (truncated after $\overrightarrow {P}\cdot \hat{u}$) into an active DDFT, separating the equations for ψ and $\overrightarrow {P}$ by integrating over φ, expanding the free energy functional, and assuming a constant mobility [37]. DDFT, in turn, can be derived from the microscopic dynamics of Brownian particles using the so-called 'adiabatic approximation', which is a closure relation that allows to relate the two-body density to the functional derivative of the equilibrium free energy functional [26, 75, 77].

The active PFC model reads [1, 37]

$$\begin{equation}{\partial }_{t}\psi \left( \overrightarrow {r},t\right)=\tilde {M}{ \overrightarrow {\nabla }}^{2}\frac{\delta F}{\delta \psi \left( \overrightarrow {r},t\right)}-{\tilde {v}}_{0} \overrightarrow {\nabla }\cdot \overrightarrow {P}\left( \overrightarrow {r},t\right),\end{equation} \tag{ 32 }$$

$$\begin{equation}{\partial }_{t} \overrightarrow {P}\left( \overrightarrow {r},t\right)=\left(\tilde {M}{ \overrightarrow {\nabla }}^{2}-{\tilde {D}}_{r}\right)\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},t\right)}-{\tilde {v}}_{0} \overrightarrow {\nabla }\psi \left( \overrightarrow {r},t\right)\end{equation} \tag{ 33 }$$

with the free-energy functional F, activity parameter (a rescaled self-propulsion velocity) ${\tilde {v}}_{0}$, and rotational diffusion constant ${\tilde {D}}_{r}$. For generality and later convenience, we have introduced a constant mobility $\tilde {M}{ >}0$, which is set to one in most treatments of the active PFC model, but which here will be useful for taking under- and overdamped limits.

A generalization toward time-delay dynamics reads

$$\begin{align}\begin{aligned}{\partial }_{t}\psi \left( \overrightarrow {r},t\right)& ={\int }_{0}^{t}\enspace \mathrm{d}s\enspace \enspace {\text{e}}^{-{\gamma }_{T}\left(t-s\right)}\left(M{ \overrightarrow {\nabla }}^{2}\frac{\delta F}{\delta \psi \left( \overrightarrow {r},s\right)}-{v}_{0} \overrightarrow {\nabla }\cdot \overrightarrow {P}\left( \overrightarrow {r},s\right)\right),\end{aligned}\end{align} \tag{ 34 }$$

$$\begin{align}\begin{aligned}{\partial }_{t} \overrightarrow {P}\left( \overrightarrow {r},t\right)& ={\int }_{0}^{t}\enspace \mathrm{d}s\enspace \enspace {\text{e}}^{-{\gamma }_{T}\left(t-s\right)}\left(M{ \overrightarrow {\nabla }}^{2}\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},s\right)}-{v}_{0} \overrightarrow {\nabla }\psi \left( \overrightarrow {r},s\right)\right)\qquad \\ & \quad \enspace -{\int }_{0}^{t}\enspace \mathrm{d}s\enspace \enspace {\text{e}}^{-{\gamma }_{R}\left(t-s\right)}{D}_{r}\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},s\right)},\end{aligned}\end{align} \tag{ 35 }$$

where we have introduced the rescaled coefficients $M={\gamma }_{T}\tilde {M}$, ${v}_{0}={\gamma }_{T}{\tilde {v}}_{0}$, and ${D}_{r}={\gamma }_{R}{\tilde {D}}_{r}$. This ansatz assumes that the terms arising from convection and translational diffusion relax on a time scale ${\gamma }_{T}^{-1}$, whereas the term arising from rotational diffusion relaxes on a time scale ${\gamma }_{R}^{-1}$. A microscopic derivation of the active PFC model from DDFT [37] shows that the dynamics of ψ and the conserved part of the dynamics of $\overrightarrow {P}$ arise from translational motion, whereas the nonconserved part of the dynamics of $\overrightarrow {P}$ is due to rotations. As discussed in section 2, the exponential relaxation of the memory kernel is related to a linear restoring force. Hence, this ansatz corresponds to the assumption that the particles are subject to a linear damping force and a linear damping torque with damping coefficients γ_T and γ_R , respectively. This assumption was also made in the microscopic model for underdamped active matter from reference [2], where exponential relaxations with different time scales for translational and rotational contributions were also observed. Note that equation (34) still has the form of a continuity equation for $\psi \left( \overrightarrow {r},t\right)$. Since all further steps are formally exact, the transformation to the jerky form therefore does not affect the fact that $\psi \left( \overrightarrow {r},t\right)$ is conserved.

While equation (34) can be treated as a second-order dynamics in the standard way since it only depends on one time scale (see section 2), equation (35) gives rise to a third-order dynamics. To connect to the general structure presented in section 4, we identify

$$\begin{equation} \overrightarrow {x}\left( \overrightarrow {r},t\right)= \overrightarrow {P}\left( \overrightarrow {r},t\right),\end{equation} \tag{ 36 }$$

$$\begin{equation} \overrightarrow {T}\left( \overrightarrow {r},t\right)=M{ \overrightarrow {\nabla }}^{2}\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},t\right)}-{v}_{0} \overrightarrow {\nabla }\psi \left( \overrightarrow {r},t\right),\end{equation} \tag{ 37 }$$

$$\begin{equation} \overrightarrow {R}\left( \overrightarrow {r},t\right)=-{D}_{r}\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},t\right)}.\end{equation} \tag{ 38 }$$

We combine equations (36)–(38) with equation (28) to the jerky active matter model

$$\begin{align}\begin{aligned}{\partial }_{t}^{2}\psi \left( \overrightarrow {r},t\right)& =-{\gamma }_{T}{\partial }_{t}\psi \left( \overrightarrow {r},t\right)+M{ \overrightarrow {\nabla }}^{2}\frac{\delta F}{\delta \psi \left( \overrightarrow {r},t\right)}-{v}_{0} \overrightarrow {\nabla }\cdot \overrightarrow {P}\left( \overrightarrow {r},t\right),\end{aligned}\end{align} \tag{ 39 }$$

$$\begin{align}\begin{aligned}{\partial }_{t}^{3} \overrightarrow {P}\left( \overrightarrow {r},t\right)& =-\left({\gamma }_{T}+{\gamma }_{R}\right){\partial }_{t}^{2} \overrightarrow {P}\left( \overrightarrow {r},t\right)-{\gamma }_{T}{\gamma }_{R}{\partial }_{t} \overrightarrow {P}\left( \overrightarrow {r},t\right)-\left({\gamma }_{T}+{\partial }_{t}\right){D}_{r}\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},t\right)}\\ & \quad \enspace +\left({\gamma }_{R}+{\partial }_{t}\right)\left(M{ \overrightarrow {\nabla }}^{2}\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},t\right)}-{v}_{0} \overrightarrow {\nabla }\psi \left( \overrightarrow {r},t\right)\right).\end{aligned}\end{align} \tag{ 40 }$$

As a consistency check, we confirm that the dynamics of $\overrightarrow {P}$ is of second order for γ_T = γ_R = γ. In this case, by the line of argument presented in section 2, it should read

$$\begin{equation}\begin{aligned}{\partial }_{t}^{2} \overrightarrow {P}\left( \overrightarrow {r},t\right)& =-\gamma {\partial }_{t} \overrightarrow {P}\left( \overrightarrow {r},t\right)+\left(M{ \overrightarrow {\nabla }}^{2}-{D}_{r}\right)\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},t\right)}-{v}_{0} \overrightarrow {\nabla }\psi \left( \overrightarrow {r},t\right).\end{aligned}\end{equation} \tag{ 41 }$$

Equation (40) can, for γ_T = γ_R = γ, be written as

$$\begin{equation}\begin{aligned} \overrightarrow {0}& =\left({\partial }_{t}+\gamma \right)\left(-{\partial }_{t}^{2} \overrightarrow {P}\left( \overrightarrow {r},t\right)-\gamma {\partial }_{t} \overrightarrow {P}\left( \overrightarrow {r},t\right)\right.\\ & \quad \left.\qquad \qquad \enspace +\left(M{ \overrightarrow {\nabla }}^{2}-{D}_{r}\right)\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},t\right)}-{v}_{0} \overrightarrow {\nabla }\psi \left( \overrightarrow {r},t\right)\right).\end{aligned}\end{equation} \tag{ 42 }$$

If we impose equation (41) as an initial condition at t = 0, equation (42) ensures that equation (41) is satisfied at all times, such that we have a second-order dynamics. (The necessity of imposing an additional initial condition arises because a differential equation of third order in temporal derivatives requires an initial condition for the second temporal derivative.)

We can also consider the over- and underdamped limits. The overdamped limit corresponds to γ_T , γ_R → ∞ with fixed $\tilde {M}$, ${\tilde {v}}_{0}$, and ${\tilde {D}}_{r}$. Thereby, the standard active PFC model given by equations (32) and (33) is recovered. The underdamped limit, on the other hand, corresponds to γ_T , γ_R → 0 at fixed M, v₀, and D_r . In this case, after integrating equation (40) over t, we find

$$\begin{equation}{\partial }_{t}^{2}\psi \left( \overrightarrow {r},t\right)=M{ \overrightarrow {\nabla }}^{2}\frac{\delta F}{\delta \psi \left( \overrightarrow {r},t\right)}-{v}_{0} \overrightarrow {\nabla }\cdot \overrightarrow {P}\left( \overrightarrow {r},t\right),\end{equation} \tag{ 43 }$$

$$\begin{equation}{\partial }_{t}^{2} \overrightarrow {P}\left( \overrightarrow {r},t\right)=\left(M{ \overrightarrow {\nabla }}^{2}-{D}_{r}\right)\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},t\right)}-{v}_{0} \overrightarrow {\nabla }\psi \left( \overrightarrow {r},t\right).\end{equation} \tag{ 44 }$$

Other limits are also possible. For example, we can consider the case in which translational degrees of freedom are underdamped, whereas rotational degrees of freedom are overdamped. This corresponds to γ_T → 0 and γ_R → ∞ with keeping M, v₀, and ${\tilde {D}}_{r}$ fixed. We then obtain

$$\begin{equation}{\partial }_{t}^{2}\psi \left( \overrightarrow {r},t\right)=M{ \overrightarrow {\nabla }}^{2}\frac{\delta F}{\delta \psi \left( \overrightarrow {r},t\right)}-{v}_{0} \overrightarrow {\nabla }\cdot \overrightarrow {P}\left( \overrightarrow {r},t\right),\end{equation} \tag{ 45 }$$

$$\begin{equation}{\partial }_{t}^{2} \overrightarrow {P}\left( \overrightarrow {r},t\right)=\left(M{ \overrightarrow {\nabla }}^{2}-{\tilde {D}}_{r}{\partial }_{t}\right)\frac{\delta F}{\delta \overrightarrow {P}\left( \overrightarrow {r},t\right)}-{v}_{0} \overrightarrow {\nabla }\psi \left( \overrightarrow {r},t\right).\end{equation} \tag{ 46 }$$

Note that the dynamics of $\overrightarrow {P}$, given by equation (46), is still of second order in the case of completely overdamped orientational dynamics as long as the translational dynamics is underdamped. This, however, is plausible since the local polarization changes not only due to rotations of the individual particles, but also due to translational motion.

It is interesting to compare equation (46) to the polarization dynamics given by equation (5) in reference [17], which reads (in our notation)

$$\begin{equation}{\partial }_{t} \overrightarrow {P}\left( \overrightarrow {r},t\right)=\left({C}_{1}{ \overrightarrow {\nabla }}^{2}-{\tilde {D}}_{r}\right) \overrightarrow {P}-{\tilde {v}}_{0} \overrightarrow {\nabla }\psi \left( \overrightarrow {r},t\right)\end{equation} \tag{ 47 }$$

with a constant C₁. This equation, which is of first order in temporal derivatives, also arises in the context of an underdamped PFC model. However, equation (47) is derived by neglecting orientational convection and by introducing a gradient term proportional to ${\left( \overrightarrow {\nabla } \overrightarrow {P}\right)}^{2}$ in the free energy. Consequently, the gradient term in equation (47) has a different physical origin than the gradient terms in equation (33), despite the fact that it looks very similar.

A further approach that our result can be compared to is the hydrodynamic active PFC model presented in the supplemental material of reference [1], the first article on active PFC models. In the standard active PFC model, it is assumed that the particles are immersed in an overdamped background fluid that is at rest, such that it exerts a viscous drag force on the particles but is not affected by them. Menzel and Löwen [1] incorporated hydrodynamic interactions by adding convective terms and a hydrodynamic equation for the velocity field $\overrightarrow {v}\left( \overrightarrow {r},t\right)$ describing the incompressible fluid flow. This approach is very different from ours in that it assumes that the active particles are overdamped while the background fluid is not, whereas our model considers translationally and rotationally underdamped active particles in a background fluid that is not modeled explicitly. An even more general approach would be to couple the equations of motion for underdamped active particles to those of an underdamped background fluid. This, of course, is only required if hydrodynamic interactions are relevant for the system one wants to describe.

In the following, we drop the dependence on space and time and restrict ourselves to one spatial dimension. Following reference [42], we use the Swift–Hohenberg [112] free-energy functional

$$\begin{equation}F=\int \enspace \mathrm{d}x\enspace \left(\frac{1}{2}\left(\psi \left({\epsilon}+{\left(1+{\partial }_{x}^{2}\right)}^{2}\right)\psi \right)+\frac{1}{4}{\left(\psi +\bar{\psi }\right)}^{4}+\frac{b}{2}{P}^{2}\right)\end{equation} \tag{ 48 }$$

with the constant coefficients , $\bar{\psi }$, and b. The parameter is a rescaled shifted temperature that some authors refer to as 'undercooling' [109]. Together with the parameter $\bar{\psi }$, which is the mean density, it determines whether the liquid phase is stable. Therefore, the Swift–Hohenberg free energy allows to describe the liquid-solid transition [42]. Finally, the parameter b, which we assume to be positive, leads to the relaxation of the polarization by making the state $\overrightarrow {P}= \overrightarrow {0}$ energetically favorable. In the simplest case, it is just proportional to the temperature, but it can be modified by interactions (see reference [35] for a microscopic expression). It is also possible that b < 0, which leads to spontaneous polarization. In this case, one would add a stabilizing term proportional to ${ \overrightarrow {P}}^{4}$ to the free energy [42]. The free energy (48) leads to

$$\begin{equation}\frac{\delta F}{\delta \psi }=\left({\epsilon}+{\left(1+{\partial }_{x}^{2}\right)}^{2}\right)\psi +{\left(\psi +\bar{\psi }\right)}^{3},\end{equation} \tag{ 49 }$$

$$\begin{equation}\frac{\delta F}{\delta P}=bP.\end{equation} \tag{ 50 }$$

Next, we consider small deviations from a reference state (ψ, P) = (0, 0) in the form

$$\begin{equation}\psi ={\psi }_{1}\enspace \mathrm{exp}\left(\lambda t-\mathrm{i}kx\right),\end{equation} \tag{ 51 }$$

$$\begin{equation}P={P}_{1}\enspace \mathrm{exp}\left(\lambda t-\mathrm{i}kx\right)\end{equation} \tag{ 52 }$$

with amplitudes ψ₁ and P₁, growth rate λ, imaginary unit i, and wave number k. We insert equations (49) and (50) together with the ansatz given by equations (51) and (52) into equations (39) and (40), linearize, and obtain

$$\begin{equation}{\lambda }^{2}{\psi }_{1}=-{\gamma }_{T}\lambda {\psi }_{1}-M{k}^{2}\left({\epsilon}+3{\bar{\psi }}^{2}+{\left(1-{k}^{2}\right)}^{2}\right){\psi }_{1}+\mathrm{i}{v}_{0}k{P}_{1},\end{equation} \tag{ 53 }$$

$$\begin{equation}{\lambda }^{3}{P}_{1}=-{\gamma }_{T}{\gamma }_{R}\lambda {P}_{1}-\left({\gamma }_{T}+{\gamma }_{R}\right){\lambda }^{2}{P}_{1}-\left({\gamma }_{T}+\lambda \right){D}_{r}b{P}_{1}-\left({\gamma }_{R}+\lambda \right)M{k}^{2}b{P}_{1}+\mathrm{i}{v}_{0}\left({\gamma }_{R}+\lambda \right)k{\psi }_{1}.\end{equation} \tag{ 54 }$$

When writing equations (53) and (54) as an eigenvalue problem in the form

$$\begin{equation}{\lambda }^{3}\left(\begin{matrix}\hfill {\psi }_{1}\hfill \\ \hfill {P}_{1}\hfill \end{matrix}\right)=\underline{M}\left(\begin{matrix}\hfill {\psi }_{1}\hfill \\ \hfill {P}_{1}\hfill \end{matrix}\right)\end{equation} \tag{ 55 }$$

with the matrix

$$\begin{equation}\underline{M}=\left(\begin{matrix}\hfill -{\gamma }_{T}{\lambda }^{2}-M{k}^{2}\left({\epsilon}+3{\bar{\psi }}^{2}+{\left(1-{k}^{2}\right)}^{2}\right)\lambda \hfill & \hfill \mathrm{i}{v}_{0}\lambda k\hfill \\ \hfill \mathrm{i}{v}_{0}\left({\gamma }_{R}+\lambda \right)k\hfill & \hfill -{\gamma }_{T}{\gamma }_{R}\lambda -\left({\gamma }_{T}+{\gamma }_{R}\right){\lambda }^{2}-\left({\gamma }_{T}+\lambda \right){D}_{r}b-\left({\gamma }_{R}+\lambda \right)M{k}^{2}b\hfill \end{matrix}\right),\end{equation} \tag{ 56 }$$

we get the characteristic polynomial determining the dispersion ³ :

$$\begin{equation}\begin{aligned}0& =\left(-{\gamma }_{T}\lambda -M{k}^{2}\left({\epsilon}+3{\bar{\psi }}^{2}+{\left(1-{k}^{2}\right)}^{2}\right)-{\lambda }^{2}\right)\\ & \quad \enspace {\times}\left(-{\gamma }_{T}{\gamma }_{R}\lambda -\left({\gamma }_{T}+{\gamma }_{R}\right){\lambda }^{2}-\left({\gamma }_{T}+\lambda \right){D}_{r}b-\left({\gamma }_{R}+\lambda \right)M{k}^{2}b-{\lambda }^{3}\right)\\ & \quad \enspace +{v}_{0}^{2}\left({\gamma }_{R}+\lambda \right){k}^{2}.\end{aligned}\end{equation} \tag{ 57 }$$

To investigate the properties of the jerky active matter model, we started by performing a linear stability analysis of the homogeneous state (ψ, P) = (0, 0). This can be done by solving equation (57) for the complex growth rate λ for different values of the real wavenumber k. If Re(λ) > 0 for any value of k, the system is unstable. Since equation (57) is a fifth-order polynomial, we have solved it numerically.

The results can be found in figure 1, which shows the dispersion relation for various values of the damping parameters γ_T and γ_R . The other parameters have been chosen as M = 1, v₀ = 0.5, D_r = 0.5, = −1.5, $\bar{\psi }=-1$, and b = 0.1, similar to the values used in reference [42]. There, it was argued that only $\bar{\psi }{< }0$ is physical for colloidal systems and that freezing of the liquid state is possible if the temperature parameter is negative. In this case, the liquid state is stable for larger $\vert \bar{\psi }\vert$.

Zoom In Zoom Out Reset image size

In our extended model, we find that for γ_T = γ_R = 0.2 (figure 1(a)), i.e., for low damping, the homogeneous state is unstable for the given parameters. However, it is stable for γ_T = γ_R = 0.4 (figure 1(b)), i.e., for higher damping. This shows that stronger damping can lead to stability of an otherwise unstable state. Consequently, in an active system where inertial delay plays a role, instabilities are more likely to occur at lower damping. One also finds a solution λ = −γ, since for γ_T = γ_R = γ, one can (due to equation (42)) factor out a contribution (γ + λ) in equation (57). Since γ > 0, this does not affect the stability.

A particularly interesting aspect of our model (in which the jerky form is relevant) is the case γ_T ≠ γ_R . For γ_T = 0.2 and γ_R = 0.4 (figure 1(c)), the liquid state is also unstable. Since γ_T = γ_R = 0.4 leads to a stable liquid state, reducing only one damping parameter (in this case γ_T ) while fixing the other one can make the system unstable. Finally, for γ_T = 0.4 and γ_R = 0.2 (figure 1(d)), the local maximum of Re(λ) also increases slightly above the stability boundary compared to the situation for γ_T = γ_R = 0.4 (figure 1(b)), but the effect is significantly smaller than in figure 1(c). Consequently, changing γ_R has a much smaller effect on the stability of the system than changing γ_T .

The phase diagram, which is shown in figure 2, gives a broader picture. It shows the maximum of Re(λ(k)), i.e., the maximal growth rate, as a function of γ_T and γ_R for (a) v₀ = 0.5, (b) v₀ = 0.75, and (c) v₀ = 1. As can be seen, stronger translational damping (larger γ_T ) stabilizes the homogeneous state, whereas rotational damping (measured by γ_R ) has only weak effects on the linear stability. Moreover, comparing the phase diagrams for different values of v₀ shows that activity tends to destabilize the system ⁴ .

Zoom In Zoom Out Reset image size

It is a very remarkable observation that the values of the damping parameters affect the linear stability of the system. In the passive case, equation (53) leads to the dispersion relation

$$\begin{equation}\lambda \left(k\right)=-\frac{{\gamma }_{T}}{2}{\pm}\sqrt{\frac{{\gamma }_{T}^{2}}{4}-M{k}^{2}f\left(k\right)}\end{equation} \tag{ 58 }$$

with the function $f\left(k\right)={\epsilon}+{\left(1-{k}^{2}\right)}^{2}+3{\bar{\psi }}^{2}$ that is independent of γ_T and γ_R . The largest eigenvalue λ, obtained by choosing the '+' sign in equation (58), has a positive real part if and only if f(k) < 0 (and k ≠ 0). This reflects the fact that the passive system with damping always approaches the minimum of the free-energy functional F, implying that whether the homogeneous state is stable solely depends on whether it is a minimum of F. The damping coefficient in the passive case can only determine how and how rapidly equilibrium is approached, but not the equilibrium state itself. The active system, however, does not obey such a minimization principle ⁵ .

These observations, which indicate an influence of inertial delay on the properties of an active matter system in its steady state, are in agreement with recent results by other authors: Caprini and Marini Bettolo Marconi [104] have observed that the spatial velocity correlations in an active system are affected by the translational damping. They also report that it makes a difference in active systems whether one changes the damping by varying the mass at fixed drag coefficient or by varying the drag coefficient at fixed mass. In contrast, the relaxation time in a passive system is governed by their relation, which does not affect the equilibrium properties. Moreover, Mandal et al [14] have found that motility-induced phase separation in underdamped active systems can be suppressed by increasing the inertial time scale determined by m.

Next, we discuss sound propagation, which is a further possible application of underdamped PFC models [46, 47]. Sound waves (propagating density perturbations) in active matter are known to differ from those in passive systems in interesting ways [113]. In the passive case, sound waves propagate due to the presence of inertia, and studying them requires hydrodynamic equations in which the momentum density is a relevant variable [24, 95]. However, in active systems, sound propagation is also possible in an overdamped system, as can be shown using the Toner–Tu model [114]. Here, the density is coupled to a dynamic equation for the polarization, which allows to derive a second-order equation for linear density variations and thus for traveling sound waves [73, 74]. This effect can also be described in our model, since it also couples the density to the polarization. However, since our model can also incorporate effects of inertia, it still contains sound modes in the passive limit, which then correspond to usual sound waves ⁶ . Hence, underdamped active systems, as described by the jerky active matter model, should allow for different types of sound waves.

Sound can also be studied by considering a linear perturbation, i.e., by using the characteristic polynomial (57), obtained from the ansatz given by equations (51) and (52). However, the physical interpretation and thus the procedure is slightly different. In section 6, we have solved equation (57) for the complex growth rate λ, given a fixed real value of k. Physically, this means that the system is subject to a perturbation ∝exp(ikx) at t = 0, and we then analyse whether this perturbation grows or decays in time. By considering all values of k (dispersion relation), we then obtain the linear stability, since all linear perturbations can be decomposed in Fourier modes exp(ikx). In the study of sound waves, on the other hand, we are interested in the way in which an oscillation exp(iωt) with angular frequency ω propagates through the system. This implies that we solve equation (57) for the complex wavenumber ⁷ k = q + iκ ($\text{with}\enspace q,\kappa \in \mathbb{R}$), given a fixed λ = iω that is purely imaginary. The imaginary part κ of k then gives the inverse of the length scale on which the sound wave decays (attenuation length).

Equation (57), which contains only even powers of k, then produces four types of solutions, namely (1) ω, q > 0, (2) ω, q < 0, (3) ω > 0 and q < 0, and (4) ω < 0 and q > 0. Here, solutions of types (1) and (2) describe waves traveling to the right (to +∞), whereas solutions of types (3) and (4) describe waves traveling to the left (to −∞). For reasons of symmetry, we can restrict ourselves to solutions of type (1). In this case, solutions with κ < 0 describe waves that start from x = 0 and travel to the right with an amplitude that decreases in the direction of propagation. On the other hand, a solution with κ > 0 would correspond to a wave that travels to the right with an amplitude that increases in the direction of propagation. The latter case is not a reasonable description of a propagating sound wave, such that we focus our discussion on solutions with κ < 0 (although other solutions are also shown for completeness).

First, we discuss a passive system (v₀ = 0) with very weak damping (γ_T = 0.01, γ_R = 0.02). This corresponds to the conditions under which sound propagation is usually investigated: for v₀ = 0, equations (39) and (40) describing ψ and P, respectively, decouple. Both equations have traveling waves as a solution. The ψ-waves correspond to usual passive sound, i.e., to density perturbations that propagate due to the coupling of density and momentum. On the other hand, the P-waves are propagating deviations from the equilibrium polarization. We here assume that the polarization can still be defined ⁸ in a meaningful way for v₀ = 0.

The results (ω(q) and κ(ω)) for parameters (a) = −1.5 and b = 0.1, (b) = −2.5 and b = 0.1, and (c) = −1.5 and b = 0.5 are shown in figure 3, where the other parameters are fixed to M = 1, D_r = 0.5, and $\bar{\psi }=-1$. In the first row, ω(q) is shown (figures 3(a.i)–(c.i)), the second row shows the corresponding curves for κ(ω) (figures 3(a.ii)–(c.ii)). Four modes can be found: first, there is a red mode that propagates almost without damping (κ ≈ 0). Second, there is a green mode for which q(ω) is small at small ω, followed by a region with stronger growth at larger ω. This mode is damped at small ω and propagates further into the medium at larger ω. Third, there are a blue mode and an orange mode where ω(q) has (for these parameters) a strong negative slope and reaches ω = 0 for q ≠ 0. The damping of the orange mode has no strong dependence on the frequency, whereas the blue mode has positive κ and is therefore not a physically reasonable sound wave (we plot it as a dashed line to indicate this).

Zoom In Zoom Out Reset image size

When comparing the various plots, the physical origin of the different modes becomes transparent: the green mode is identical in figures 3(a) and (b), whereas the other modes show a significant change: in figure 3(b.i), they intersect at ω ≈ 0.6, which is not the case in figure 3(a.i). At ω ≈ 0.6, one also observes a maximum for the orange and a minimum for the blue mode in figure 3(b.ii), which is not the case in figure 3(a.ii). On the other hand, the red, blue, and orange modes do not change between figures 3(a) and (c), whereas the growth of q with ω in the green mode now starts at ω ≈ 0.5 in figure 3(c.i) rather than ω ≈ 0.2 as in figure 3(a.i). Since only , a parameter affecting ψ, has been changed between figures 3(a) and (b), it follows that the red, blue, and orange modes correspond to density waves. The red mode is usual passive sound, whereas the orange and blue modes arise from the higher-order spatial gradients in the free energy (48). (We have verified that the orange and blue modes disappear if these terms are not present in equation (48).) As far as ω(q) is concerned, the orange and blue modes coincide, which is not the case for κ(ω). On the other hand, only the parameter b, which affects P, has been changed between figures 3(a) and (c), such that the green mode corresponds to polarization waves.

The situation is quite different in the active case (v₀ = 0.5), which is shown in figure 4. First, we need to take into account here that, as shown in section 6, activity tends to destabilize the system. Consequently, we have to increase the damping if we want sound waves in a stable fluid. Here, we choose γ_T = 0.9 and γ_R = 1. Apart from the change in damping and activity, we use the same parameters as in figure 3.

Zoom In Zoom Out Reset image size

When comparing figures 4(a)–(c), it is found that a change of both and b affects all four modes. The reason is that, in the active case, oscillations of ψ and P are coupled. Moreover, the shapes of the modes change. For the red mode, the most significant change in ω(q) is a sharp bend observed at ω ≈ 0.8 in figure 4(c.i). Moreover, the red mode is now damped for larger ω, whereas it can propagate almost without damping at small ω (figures 4(a.ii)–(c.ii)). When considering the green mode, strong variations can be found for both ω(q) and κ(ω). For = −1.5 (figures 4(a) and (c)), both functions grow monotonously. In contrast, for = −2.5, ω(q) is not even a function due to a local maximum in q(ω) (figure 4(b.i)) and κ(ω) reaches a local maximum at ω ≈ 0.4 (figure 4(b.ii)). A very significant change compared to figure 3 is also observed for the orange mode: now, ω(q) has a 'hockey-stick shape', where for larger frequencies the mode grows to the right (and not to the left as in the passive case). Moreover, ω(q) does not coincide for the blue and orange mode in the active case.

Consequently, if we define 'sound' as 'propagating density oscillation', we find two different mechanisms for sound propagation in active matter: first, there are oscillations due to translational inertia as in the passive case. These are still present in active matter, although their properties are modified by the coupling to polarization degrees of freedom. Second, there are oscillations that arise from the coupling of the density to the polarization. This mechanism is not present in fluids consisting of passive spherical particles. As can be seen in figures 4(a.ii) and (c.ii), the damping of the green mode decreases for larger frequencies if we set = −1.5, whereas the 'passive' sound modes are more strongly damped in this case. Consequently, sound can propagate further into the medium via the second, 'active' mechanism at higher frequencies at these parameter values.

We summarize our observations for sound propagation in active matter:

• For v₀ = 0, there is one mechanism for sound propagation (coupling of density and momentum), as is well known from passive fluids.
• For active systems, stronger damping is required to ensure that the system is still stable against perturbations.
• A second type of sound waves exists in the active case, arising from the coupling of density and polarization. (For v₀ = 0, these waves only affect the polarization, not the density.)
• For certain parameter values, waves of the first type propagate further into the medium at small frequencies, whereas waves of the second type propagate further at higher frequencies.

The last effect is of particular interest for possible applications, since it allows for the development of acoustic frequency filters: if one excites oscillations in the red mode, smaller frequencies will propagate further into the medium, i.e., the system constitutes a low-pass filter. On the other hand, if the green mode is excited, higher frequencies will experience less damping, such that we have a high-pass filter. The properties of the filters can be tuned by changing the system parameters as shown in figure 4.

We have shown how the presence of two relaxational time scales changes a general first-order soft matter model into a third-order spatiotemporal jerky dynamics. This general structure has then been used to derive a general active PFC model in which, in contrast to the standard case, relaxational and orientational time scales are both finite and different.

The resulting model as well as the method by which it was obtained allow to describe systems of underdamped active particles with greater accuracy than the standard active PFC model. Since our model is of third order in temporal derivatives (unlike usual theories which are of first or second order), it can be expected to show a number of interesting effects that are not present in other models. This is also plausible by comparison to the case of ordinary dynamical systems, where jerky dynamics allows for chaotic behaviour that is not possible in second-order models.

From the dispersion relation, we were able to investigate the linear stability and the resulting phase diagram, showing that in active systems, where the phase boundaries are not determined by the free energy, changes of the damping coefficients can affect the linear stability of the fluid phase. The fluid state is destabilized by increased activity. Moreover, it is found that active fluids can exhibit two different mechanisms for sound propagation: waves arising from density–momentum coupling as in the passive case, which can propagate with small damping at low frequencies, and waves arising from density–polarization coupling, which can propagate further at high frequencies. This effect has potential applications in the development of acoustic frequency filters based on active fluids.

Possible extensions of this work include the consideration of translational–rotational coupling (as discussed for particles with general shape in reference [118]). Moreover, one could consider spherical domains (as done previously for simpler PFC models [41, 119]) or more general memory kernels.

We thank Tobias Frohoff-Hülsmann, Svetlana V Gurevich, Stefan J Linz, Lukas Ophaus, Alice Rolf, and Uwe Thiele for helpful discussions. MtV thanks the Studienstiftung des deutschen Volkes for financial support. RW is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—WI 4170/3-1.

The data that support the findings of this study are available upon request from the authors.

In this appendix, we discuss the derivation of equation (28) in more detail.

We wish to re-write the dynamical system

$$\begin{equation}\dot {x}=y+z,\end{equation} \tag{ A1 }$$

$$\begin{equation}\dot {y}=-{\gamma }_{T}y+T,\end{equation} \tag{ A2 }$$

$$\begin{equation}\dot {z}=-{\gamma }_{R}z+R\end{equation} \tag{ A3 }$$

as a single third-order dynamic equation (jerky dynamics) for the variable x. For this purpose, we compute the derivatives

$$\begin{equation}\ddot {x}=\dot {y}+\dot {z},\end{equation} \tag{ A4 }$$

$$\begin{equation}\ddot {y}=-{\gamma }_{T}\dot {y}+\dot {T},\end{equation} \tag{ A5 }$$

$$\begin{equation}\ddot {z}=-{\gamma }_{R}\dot {z}+\dot {R},\end{equation} \tag{ A6 }$$

$$\begin{equation} \dddot {x}=\ddot {y}+\ddot {z}.\end{equation} \tag{ A7 }$$

One can now use equations (A1)–(A6) to express the unknown variables $\ddot {y}$ and $\ddot {z}$ in equation (A7) in terms of x, $\dot {x}$ and $\ddot {x}$.

From equations (A2) and (A5), we get

$$\begin{equation}\ddot {y}={\gamma }_{T}^{2}y-{\gamma }_{T}T+\dot {T}.\end{equation} \tag{ A8 }$$

Similarly, equations (A3) and (A6) give

$$\begin{equation}\ddot {z}={\gamma }_{R}^{2}z-{\gamma }_{R}R+\dot {R}\end{equation} \tag{ A9 }$$

and equations (A2)–(A4) give

$$\begin{equation}\ddot {x}=-{\gamma }_{T}y+T-{\gamma }_{R}z+R.\end{equation} \tag{ A10 }$$

Equations (A7)–(A9) lead to

$$\begin{equation} \dddot {x}={\gamma }_{T}^{2}y+{\gamma }_{R}^{2}z-{\gamma }_{T}T-{\gamma }_{R}R+\dot {T}+\dot {R}.\end{equation} \tag{ A11 }$$

Moreover, from equation (A1) we obtain

$$\begin{equation}y=\dot {x}-z.\end{equation} \tag{ A12 }$$

Equation (A10) gives

$$\begin{equation}y=-\frac{1}{{\gamma }_{T}}\ddot {x}+\frac{1}{{\gamma }_{T}}T-\frac{{\gamma }_{R}}{{\gamma }_{T}}z+\frac{1}{{\gamma }_{T}}R,\end{equation} \tag{ A13 }$$

and equating equations (A12) and (A13) allows to find

$$\begin{equation}z=-\frac{{\gamma }_{T}}{{\gamma }_{R}-{\gamma }_{T}}\dot {x}-\frac{1}{{\gamma }_{R}-{\gamma }_{T}}\left(\ddot {x}-T-R\right).\end{equation} \tag{ A14 }$$

At this point, we have made the assumption γ_R ≠ γ_T . From equations (A12) and (A14), we get

$$\begin{equation}\begin{aligned}y& =\dot {x}+\frac{{\gamma }_{T}}{{\gamma }_{R}-{\gamma }_{T}}\dot {x}+\frac{1}{{\gamma }_{R}-{\gamma }_{T}}\left(\ddot {x}-T-R\right)\\ & =\frac{{\gamma }_{R}}{{\gamma }_{R}-{\gamma }_{T}}\dot {x}+\frac{1}{{\gamma }_{R}-{\gamma }_{T}}\left(\ddot {x}-T-R\right).\end{aligned}\end{equation} \tag{ A15 }$$

Finally, combining equations (A11), (A14), and (A15) leads to

$$\begin{equation}\begin{aligned} \dddot {x}& =\frac{{\gamma }_{T}^{2}{\gamma }_{R}}{{\gamma }_{R}-{\gamma }_{T}}\dot {x}+\frac{{\gamma }_{T}^{2}}{{\gamma }_{R}-{\gamma }_{T}}\left(\ddot {x}-T-R\right)-{\gamma }_{T}T+\dot {T}\\ & \quad \enspace -\frac{{\gamma }_{R}^{2}{\gamma }_{T}}{{\gamma }_{R}-{\gamma }_{T}}\dot {x}-\frac{{\gamma }_{R}^{2}}{{\gamma }_{R}-{\gamma }_{T}}\left(\ddot {x}-T-R\right)-{\gamma }_{R}R+\dot {R},\end{aligned}\end{equation} \tag{ A16 }$$

which can be simplified to the final result

$$\begin{equation} \dddot {x}=-\left({\gamma }_{T}+{\gamma }_{R}\right)\ddot {x}-{\gamma }_{T}{\gamma }_{R}\dot {x}+{\gamma }_{R}T+{\gamma }_{T}R+\dot {T}+\dot {R}\end{equation} \tag{ A17 }$$

that equals equation (28). One can easily verify that equation (28) also holds for γ_T = γ_R .