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

Most field theories for active matter neglect effects of memory and inertia. However, recent experiments have found inertial delay to be important for the motion of self-propelled particles. A major challenge in the theoretical description of these effects, which makes the application of standard methods very difficult, is the fact that orientable particles have both translational and orientational degrees of freedom which do not necessarily relax on the same time scale. In this work, we derive the general mathematical form of a field theory for soft matter systems with two different time scales. This allows to obtain a phase field crystal model for polar (i.e., nonspherical or active) particles with translational and orientational memory. Notably, this theory is of third order in temporal derivatives and can thus be seen as a spatiotemporal jerky dynamics. We obtain the phase diagram of this model, which shows that, unlike in the passive case, the linear stability of the liquid state depends on the damping coefficients. Moreover, we investigate sound waves in active matter. It is found that, in active fluids, there are two different mechanisms for sound propagation. For certain parameter values and sufficiently high frequencies, sound mediated by polarization waves experiences less damping than usual passive sound mediated by pressure waves of the same frequency. By combining the different modes, acoustic frequency filters based on active fluids could be realized.

Most field theories for active matter neglect effects of memory and inertia. However, recent experiments have found inertial delay to be important for the motion of self-propelled particles. A major challenge in the theoretical description of these effects, which makes the application of standard methods very difficult, is the fact that orientable particles have both translational and orientational degrees of freedom which do not necessarily relax on the same time scale. In this work, we derive the general mathematical form of a field theory for soft matter systems with two different time scales. This allows to obtain a phase field crystal model for polar (i.e., nonspherical or active) particles with translational and orientational memory. Notably, this theory is of third order in temporal derivatives and can thus be seen as a spatiotemporal jerky dynamics. We obtain the phase diagram of this model, which shows that, unlike in the passive case, the linear stability of the liquid state depends on the damping coefficients. Moreover, we investigate sound waves in active matter. It is found that, in active fluids, there are two different mechanisms for sound propagation. For certain parameter values and sufficiently high frequencies, sound mediated by polarization waves experiences less damping than usual passive sound mediated by pressure waves of the same frequency. By combining the different modes, acoustic frequency filters based on active fluids could be realized.

I. INTRODUCTION
Almost all field theories for active matter systems are derived in the overdamped limit and neglect the inertia of the active particles. However, there is a current increase of interest in the role of memory and (inertial) delay in active matter [1][2][3][4][5][6][7]. Experiments have shown inertial [1][2][3] and sensorial [4,5] delay to be important for the dynamics of self-propelled particles. Inertia can lead to interesting effects that are not present in overdamped active systems [1][2][3][8][9][10][11], such that inertial active matter models, as presented in Refs. [10,12], 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 [3,13]. Memory and inertia effects are intimately connected and sometimes equivalent, as can be seen from an analysis in the Mori-Zwanzig framework (see Section II).
The problem of memory effects has been discussed in other contexts, very notably in particle physics [14][15][16]. 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 * Corresponding author: raphael.wittkowski@uni-muenster.de the relevant dynamics itself is very fast, such that neglecting memory effects and causality constraints is an approximation that cannot be justified [15]. For this reason, more general theories have been derived that take these effects into account [14][15][16][17].
This problem can also arise in soft matter systems such as polymer solutions [17]. 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 [18,19], who suggested this as a topic of further investigation within classical dynamical density functional theory (DDFT) [20], 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.
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 [40][41][42][43]. Moreover, third-order dispersions have been studied in optics [44,45], where they naturally appear in delay systems [46].
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 [31] 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 behavior. 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 [47,48]. 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 as follows: In Section II, we explain how finite relaxation times can be incorporated into transport equations. We provide a brief introduction to jerky dynamics in Section III. A general model for systems with two time scales is developed in Section IV. In Section V, we use these results to obtain an active PFC model with memory. We perform a linear stability analysis of this model in Section VI. Sound propagation is studied in Section VII. We conclude in Section VIII.

II. CAUSAL FIELD THEORIES
For a conserved order-parameter field φ( r, t), where r and t denote position and time, respectively, the time evolution is given by the continuity equation with the current J( r, t). A typical example is traditional DDFT [20,[49][50][51][52], where the current is with the mobilityM (which for traditional DDFT is given byM = β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 [53,54], where the current is A third example is Active Model B [55], where the current is δφ( r, t) + λ( ∇φ( r, t)) 2 + Λ( r, t) (4) with the activity parameter λ and the noise Λ( r, t).
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 r does not depend on the value of φ at a different position r . This is not the most general case. A good starting point for a systematic analysis is the Mori-Zwanzig formalism [56][57][58][59][60][61]. It allows to derive macroscopic transport equations for an arbitrary set of relevant variables from the microscopic dynamics and thereby, as discussed in Ref. [61], gives very general insights into the general structure of these equations. As exploited in recent work by Meyer et al. [62][63][64][65], 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 [52,61].
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 φ( r, t) that is a function of the microscopic degrees of freedom. Then, we write the general time evolution as 2 where M is the memory kernel. We have ignored here an organized drift term, which often vanishes for reasons of symmetry if there is only one relevant variable, and we have dropped a noise term. 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 [52]), the exact transport equation (5) can be approximately written as with the diffusion tensor 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 [66]. In the Markovian limit, an H-theorem corresponding to an increase of entropy can be proven [67,68]. From a more practical point of view, the Markovian approximation corresponds to the assumption that the set of relevant variables we have chosen gives a complete description of the macroscopic state. For example, if φ( r, t) is the number density, making the Markovian approximation implies that the momentum density g( r, t) relaxes very rapidly, i.e., that we are working in the overdamped limit [52].
What is also relevant here is that Eq. (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 [69] (see Ref. [17] for a review). This issue is particularly relevant in the case of relativistic systems. If we wish to describe phase 2 In our notation, the operators ∇ and ∂t are always understood as acting on the whole term to the right of them. For example, in the expression ∂tf (t)g(t)+h(t) with time-dependent functions f , g, and h, the operator ∂t acts on the product f (t)g(t) and not just on f (t). Obviously, ∂t does not act on h(t). 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) [15]. Thus, when applying the Mori-Zwanzig formalism to relativistic systems, we generally have to be very careful in handling memory effects [14]. 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. [15] for describing phaseseparation 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 with the damping coefficient (inverse relaxation time) γ and the modified mobility M = γ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. In the Mori-Zwanzig formalism, the memory kernel is related to the correlation of the noise [60]. The assumption of white noise leads to a Markovian dynamics, whereas colored noise gives the form (7) [15]. The time derivative of Eq. (7) is given by Differentiating Eq. (1) with respect to time and inserting Eq. (8) then gives 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 [70][71][72][73]. There are two important limiting cases of Eq. (9) one can consider. The first one is the overdamped limit γ → ∞ with M/γ = const. (corresponding to fixedM ), in which case we recover the standard Cahn-Hilliard equation The second one is the underdamped limit γ → 0 at fixed M , which gives the inertial Cahn-Hilliard equation For the DDFT current (2), Archer [18,19] has noted that this procedure leads to a causal DDFT that has the same structure as the DDFT with particle mass m for particles with inertia. A similar result was obtained by Chavanis [74]. 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 Eq. (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 Eq. (12) then leads to traditional DDFT [19], whereas the underdamped limit leads to a generalized Euler equation (without convective term). From Eq. (12), we find that M = 1/m such that M = 1/(γm). Since we know from DDFT thatM = βD, we can infer γ = 1/(βmD). Therefore, varying γ at fixed M corresponds to changing 3 D, whereas varying γ at fixedM corresponds to changing m. 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 [61] (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 [20,59]. 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 ∂ t φ = − ∇ · g, such that using ∂ t φ as a relevant variable is equivalent to adding the momentum density g to the set of relevant variables.
with a time-dependent variable x(t) and a function J. Jerky dynamics is very important for chaos theory [41,75]. Mechanically, the "jerk" is the rate of change of acceleration [78]. Intuitively, one might expect jerks to be of no importance in classical mechanics, since Newton's equation of motionẍ with the force F is of second order. (Although a jerky dynamics can obviously be obtained by taking the time derivative of Eq. (14), it would not provide any physical insights.) However, if the force has memory, i.e., if it depends on values of x orẋ at previous times, taking the time derivative of Eq. (14) can lead to an interesting jerky dynamics that contains additive terms depending solely on x. Such models can allow for chaos [43]. 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 can be written in the form (13). An example is Sprott's model R [79] with the dynamical variables x, y, and z and the constants a and b, which is a simple model for chaos. The problem is addressed in detail in Ref. [40], where the conditions under which such a transformation is possible are discussed. A useful strategy [40] is to calculate the first time derivative of Eqs. (16) and (17) and the first and second time derivatives of Eq. (15). This gives seven coupled equations for x,ẋ,ẍ, ...

IV. SOFT MATTER WITH TWO TIME SCALES
In systems of active (or nonspherical) particles, one needs to take into account both translational and orien-tational degrees of freedom. As a toy model for describing the physics of such systems, we use the dynamical equation Later, we interpret T and 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 (19). We are assuming nothing about the form of T and 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 II, we generalize Eq. (19) towards a causal dynamics with time delay 4 : There is a very important difference to the models known from the literature which we have discussed in Section II: 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 Eq. (20), which gives For γ T = γ R = γ, we could write Eq. (21) as In general, however, the sum of the first two terms in Eq. (21) is not proportional to ∂ t φ( r, t), since the translational and rotational contributions appear with different prefactors due to the different time scales on which the contributions change. Therefore, Eq. (21) 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 φ( r, t) without time convolution.
This can be obtained using the procedure introduced in Section III.
It is easily verified that Eq. (42) also holds for γ T = γ R .

V. JERKY ACTIVE MATTER
Phase field crystal (PFC) models are a useful framework for the description of soft and active matter [21,31,[80][81][82][83]. They can be obtained as a limiting case of the more complex and more general case of DDFT [26,32,39]. Reviews are given by Refs. [21] (PFC models) and [20] (DDFT), both reviews discuss the derivation of PFC models from DDFT. The order parameters are the rescaled density ψ( r, t) and the polarization P ( r, t), which arise through an orientational expansion [84,85] of the one-body density ρ( r,û, t) depending on position r and orientationû.
The active PFC model reads [31,32] ∂ t ψ( r, t) =M ∇ 2 δF δψ( r, t) −ṽ 0 ∇ · P ( r, t), with the free-energy functional F , activity parameter (a rescaled self-propulsion velocity)ṽ 0 , and rotational diffusion constantD r . For generality and later convenience, we have introduced a constant mobilityM > 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 towards time-delay dynamics reads ∂ t P ( r, t) = where we have introduced the rescaled coefficients M = γ TM , v 0 = γ Tṽ0 , and D r = γ RDr . This ansatz assumes that the terms arising from convection and translational diffusion relax on a time scale γ −1 T , whereas the term arising from rotational diffusion relaxes on a time scale γ −1 R . Note that Eq. (45) still has the form of a continuity equation for ψ( r, t). Since all further steps are formally exact, the transformation to the jerky form therefore does not affect the fact that ψ( r, t) is conserved.
While Eq. (45) can be treated as a second-order dynamics in the standard way since it only depends on one time scale (see Section II), Eq. (46) gives rise to a thirdorder dynamics. To connect to the general structure presented in Section IV, we identify R( r, t) = −D r δF δ P ( r, t) .
We combine Eqs. (47)- (49) with Eq. (42) to the jerky active matter model As a consistency check, we confirm that the dynamics of P is of second order for γ T = γ R = γ. In this case, by the line of argument presented in Section II, it should read Equation (51) can, for γ T = γ R = γ, be written as If we impose Eq. (52) as an initial condition at t = 0, Eq. (53) ensures that Eq. (52) 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 fixedM ,ṽ 0 , andD r . Thereby, the standard active PFC model given by Eqs. (43) and (44) is recovered. The underdamped limit, on the other hand, corresponds to γ T , γ R → 0 at fixed M , v 0 , and D r . In this case, after integrating Eq. (51) over t, we find 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 0 , andD r fixed. We then obtain Note that the dynamics of P , given by Eq. (57), 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 Eq. (57) to the polarization dynamics given by Eq. (5) in Ref. [12], which reads (in our notation) with a constant C 1 . This equation, which is of first order in temporal derivatives, also arises in the context of an underdamped PFC model. However, Eq. (58) is derived by neglecting orientational convection and by introducing a gradient term proportional to ( ∇ P ) 2 in the free energy. Consequently, the gradient term in Eq. (58) has a different physical origin than the gradient terms in Eqs. (44) and (57), despite the fact that it looks very similar.

VI. LINEAR STABILITY
In the following, we drop the dependence on space and time and restrict ourselves to one spatial dimension. Following Ref. [80], we use the Swift-Hohenberg [86] freeenergy functional with the constant coefficients (shifted rescaled temperature [82]),ψ (mean density), and b. This leads to Next, we consider small deviations from a reference state (ψ, P ) = (0, 0) in the form with amplitudes ψ 1 and P 1 , growth rate λ, imaginary unit i, and wave number k. We insert Eqs. (60) and (61) together with the ansatz given by Eqs. (62) and (63) into Eqs. (50) and (51), linearize, and obtain When writing Eqs. (64) and (65) as an eigenvalue problem in the form with the matrix we get the characteristic polynomial determining the dispersion 5 : 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 Eq. (68) 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 Eq. (68) is a fifth-order polynomial, we have solved it numerically.
The results can be found in Fig. 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 = 0.5, D r = 0.5, = −1.5,ψ = −1, and b = 0.1, similar to the values used in Ref. [80]. There, it was argued that onlyψ < 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 |ψ|.
In our extended model, we find that for γ T = γ R = 0.2 (Fig. 1a), i.e., for low damping, the homogeneous state is unstable for the given parameters. However, it is stable for γ T = γ R = 0.4 (Fig. 1b), 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 Eq. (53)) factor out a contribution (γ + λ) in Eq. (68). 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 (Fig. 1c), 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 (Fig. 1d), the local maximum of Re(λ) also increases slightly above the stability boundary compared to the situation for γ T = γ R = 0.4 (Fig. 1b), but the effect is significantly smaller than in Fig. 1c. Consequently, changing γ R has a much smaller effect on the stability of the system than changing γ T .
The phase diagram, which is shown in Fig. 2, gives a broader picture. It shows the maximum of Re(λ(k)), 5 To be able to write Eqs. (64) and (65) as the eigenvalue problem (66), we have multiplied Eq. (64) by λ. In Eq. (68), we have then divided by λ again. Although this trick can only be applied for λ = 0, we can easily show that Eq. (68) is also correct for λ = 0: In this case, we have λ 2 = λ 3 such that Eqs. (64) and (65) form already an eigenvalue problem whose characteristic polynomial is easily found to be given by Eq. (68).
i.e., the maximal growth rate, as a function of γ T and γ R for (a) v 0 = 0.5, (b) v 0 = 0.75, and (c) v 0 = 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 0 shows that activity tends to destabilize the system. 6 It is a very remarkable observation that the values of the damping parameters affect the linear stability of the system. In the passive case, Eq. (64) leads to the dispersion relation with the function f (k) = + (1 − k 2 ) 2 + 3ψ 2 that is independent of γ T and γ R . The largest eigenvalue λ, obtained by choosing the "+" sign in Eq. (69), 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. 7

VII. SOUND WAVES
Next, we discuss sound propagation, which is a further possible application of underdamped PFC models [33,34]. Sound waves (propagating density perturbations) in active matter are known to differ from those in passive systems in interesting ways [87]. 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 [18,68]. However, in active systems, sound propagation is also possible in an overdamped system, as can be shown using the Toner-Tu model [88]. 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 travelling sound waves [47,48]. 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. 8 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 (68), obtained from the ansatz given by Eqs. (62) and (63). However, the physical interpretation and thus the procedure is slightly different: In Section VI, we have solved Eq. (68) 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 analyze 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 Eq. (68) for the complex wavenumber k = q + iκ (with q, κ ∈ 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 (68), 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 travelling to the right (to +∞), whereas solutions of types (3) and (4) describe waves travelling 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 = 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 = 0, Eqs. (50) and (51) describing ψ and P , respectively, decouple. Both equations have travelling 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 9 in a meaningful way for v 0 = 0. 9 For an active Brownian sphere, the orientation might only arise through the activity [84], such that there are no orientational degrees of freedom left for v 0 = 0 and P cannot be defined. In practice, however, an asymmetry can still be present even if it does not result in self-propulsion. As an example, consider a nonspherical particle that is propelled by ultrasound [89] -it still has a (geometrical) orientation if the ultrasound is switched off. If the polarization cannot be defined for v 0 = 0, modes linked to the polarization have a physical meaning only for v 0 = 0.
The results (ω(q) and κ(ω)) for parameters (a) =  Fig. 3, where the other parameters are fixed to M = 1, D r = 0.5, and ψ = −1. In the first row, ω(q) is shown (Figs. 3a(i)-3c(i)), the second row shows the corresponding curves for κ(ω) (Figs. 3a(ii)-3c(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).
When comparing the various plots, the physical origin of the different modes becomes transparent: The green mode is identical in Fig. 3a and Fig. 3b, whereas the other modes show a significant change: In Fig. 3b(i), they intersect at ω ≈ 0.6, which is not the case in Fig. 3a(i). At ω ≈ 0.6, one also observes a maximum for the orange and a minimum for the blue mode in Fig. 3b(ii), which is not the case in Fig. 3a(ii). On the other hand, the red, blue, and orange modes do not change between Fig. 3a and Fig. 3c, whereas the growth of q with ω in the green mode now starts at ω ≈ 0.5 in Fig. 3c(i) rather than ω ≈ 0.2 as in Fig. 3a(i). Since only , a parameter affecting ψ, has been changed between Figs. 3a and 3b, 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 (59). (We have verified that the orange and blue modes disappear if these terms are not present in Eq. (59).) 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 Fig. 3a and Fig. 3c, such that the green mode corresponds to polarization waves. The situation is quite different in the active case (v 0 = 0.5), which is shown in Fig. 4. First, we need to take into account here that, as shown in Section VI, activity tends to destabilize the system. Consequently, we have to increase the damping if we want sound waves in a fluid that does not spontaneously become a crystal. Here, we choose γ T = 0.9 and γ R = 1. Apart from the change in damping and activity, we use the same parameters as in Fig. 3.
When comparing Figs. 4a-4c, 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 Fig. 4c(i). Moreover, the red mode is now damped for larger ω, whereas it can propagate almost without damping at small ω ( Fig. 4a(ii)-Fig. 4c(ii)). When considering the green mode, strong variations can be found for both ω(q) and κ(ω). For = −1.5 (Figs. 4a and 4c), both functions grow monotonously. In contrast, for = −2.5, ω(q) is not even a function due to a local maximum in q(ω) (Fig. 4b(i)) and κ(ω) reaches a local maximum at ω ≈ 0.4 (Fig. 4b(ii)). A very significant change compared to Fig. 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 Figs. 4a(ii) and 4c(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 = 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 = 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 Fig. 4.

VIII. CONCLUSIONS
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 active matter systems in which inertia is important with greater accuracy than previous theories. 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 behavior that is not possible in secondorder 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 densitypolarization 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 Ref. [90]). Moreover, one could consider spherical domains (as done previously for simpler PFC models [91,92]) or more general memory kernels.