Mean field approach of dynamical pattern formation in underdamped active matter with short-ranged alignment and distant anti-alignment interactions

Many active matter systems, especially on the microscopic scale, are well approximated as overdamped, meaning that any inertial momentum is immediately dissipated by the environment. On the other hand, especially for macroscopic active systems but also for many mesoscopic systems the time scale of translational inertial motion can become large enough to be relevant for the dynamics. This raises the question how collective dynamics and the resulting states in active matter are influenced by inertia. Therefore, we propose a coarse-grained continuum model for underdamped active matter based on a mean field description for passive systems. Furthermore, as an example, we apply the model to a system with interactions that support an alignment on short distances and an anti-alignment on longer length scales as known in the context of pattern formation due to orientational interactions. Our numerical calculations of the under- and overdamped dynamics both predict a structured laning state. However, activity induced convective flows that are only present in the underdamped model destabilize this state when the anti-alignment is weakened, leading to a collective motion state which does not occur in the overdamped limit. A turbulent transition regime between the two states can be characterized by strong density fluctuations and the absence of global ordering.


Introduction
Active matter denotes a whole plethora of biological and arti cial many-body systems composed of self-propelled, interacting particles [1][2][3]. Often, microscopic systems are considered in which the constituents move in a viscous 1 Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. environment. Typical examples include protein laments in motility assays [4], cell colonies [5][6][7], biological or articial microswimmers [8][9][10], all surrounded by a viscous background uid at low Reynolds numbers or binding to an adhesive substrate. These systems are overdamped, meaning the motion of particles is solely governed by the momentary forces acting on them. However, Newton's rst law formulates that massive objects resist any change in momentum due to their inertia. While not relevant on the microscopic length scale in solution, inertial effects become relevant in typical macroscopic realizations of active matter like ocks of birds [11,12] or arti cially made massive robots moving on a twodimensional plate [13][14][15]. For the latter, the in uence of inertia on the long time statistics was observed recently for a one particle system [16]. Another interesting study proposes that the pressure of active particles exerted on a moving wall varies with inertia [17].
There is a wide range of particle based active matter models which explicitly incorporate the microscopic interactions [18][19][20][21][22]. Another approach to such systems is the formulation of continuum models which typically describe systems on a length scale above the particle resolution. Such coarse-grained models are already able to predict experimental observations of microscopic systems like bacterial colonies [23], ensembles of microswimmers [9,24] or active nematics [3,25,26]. However, for other systems where the time scale of inertial motion becomes relevant an extended underdamped description of the dynamics becomes necessary.
In this article, we propose a continuum model for classical underdamped active matter by including activity to a dynamical density functional theory for underdamped passive systems [27]. Then we use the model to investigate an example system where we capture the consequences of introducing translational inertia by comparing the found states to those of the corresponding overdamped limit system. The interaction rules of particle orientations we specify on the continuum level are derived from microscopic interactions and discussed in detail in [21,22]. Consequently, we nd in our system several of the states reported there. However, as discussed later on, differences arise due to the varying dynamics. Speci cally, we note that in our approach the actual velocity of a particle that is the cause of inertia effects does not have to coincide with the orientation of a particle that determines the direction of the intrinsic self-propulsion.
The article is structured as follows: in section 2 we derive the model and specify the system. Simulation results on the latter are discussed in section 3 and concluded in section 4.

Derivation
In this section, we outline the derivation of the general model which follows [27] concerning the calculation steps and coincides with it when active components are neglected. These additional active terms do not change the principal calculation but give additional contributions in each step. The inclusion of activity follows similar approaches used for overdamped active systems [10,28]. Afterward, we introduce the speci c system which is investigated later on.
We consider a general system which consists of N identical polar particles of mass m with position and momentum coordinates r N = {r 1 , r 2 , . . . , r N } and p N = {p 1 , p 2 , . . . , p N }. All particles self-propel with an active force f 0 into the direction of their orientation given by the unit vectorŝ u N = {û 1 ,û 2 , . . . ,û N }. For convenience, we restrict the derivation to effectively two dimensional systems where an angle ϕ i determines the orientationû i = (cos(ϕ i ), sin(ϕ i )) ⊤ . A particle's phase space coordinates are summarized as (2) (r N ) summarizes an external potential force F (1) (r i ) = −∇ r i V ext (r i ) and particle interactions F (2) (2) (r i , r j ) with the pair interaction potential V (2) (r i , r j ). Analogously to these translational forces particles change their orientation due to the torque G i (r N , ϕ N ) = G (1) (r i , ϕ i ) + G (2) (r N , ϕ N ) with a one particle contribution G (1) (r i , ϕ i ) and a pairwise interaction torque G (2) (r N , ϕ N ). The corresponding Langevin equations for the underdamped motion of the particle ensemble then read where γ = α/m is the damping constant with friction α and D, D R are translational and rotational diffusion constants.
and analogously for the components of Γ i . Note that the damping γ and diffusion D are not necessarily related via the Stokes-Einstein relation since they might describe, e.g., the interaction with a substrate instead of a thermal bath. We also want to point out that here we only consider translational inertia but that we neglect inertia effects in angular direction, i.e. we consider a case where the inertia time scale of the rotation is negligible while the inertia time scale of the translational motion still is considered. In case of particles that have to be rotated in order to change the orientation, this is the case if the mass of the particle is mainly located close to the center of the particle or in the limit of point particles as in [21,29]. These microscopic equations of motion for the stochastic state variables are converted into the corresponding Fokker-Planck equation of the N body phase space probability density f (N ) (x N , t) which gives us the probability of nding the system in a con guration x N = {r N , p N ,û N } at time t [30]. The resulting dynamical equation reads This N body problem is simpli ed by de ning the n body reduced phase space distribution functions where N − n of the N bodies' state variables are integrated out. Equivalently, integrating the N body Fokker-Planck equation over N − 1 sets of particle variables yields the dynamical equation for the one body distribution − dr 2 dp 2 dϕ 2 ∂ ϕ 1 G (2) f (2) where we assume that the N body density and its rst derivatives decay to zero for all r i , p i → ∞ and are periodic in ϕ i . The latter assumption is also used for the torques G 1 , G 2 . The reduced one body phase space density f (1) (x 1 , t) gives the probability of nding a particle at the position r 1 with momentum p 1 and orientationû 1 at time t. Note that due to particle interactions, (4) still depends on the two body density f (2) For the next step, we de ne the mean eld quantities number density ρ, momentum current j and orientation current ρP with the mean orientation eld P as ρ P(r 1 , t) = dp 1 dϕ 1û1 f (1) (x 1 , t).
Our goal is now to reduce the equation for f (1) to equations for the mean elds and simplify the problem further from there. Equation (4) may be integrated over p 1 and ϕ 1 which yields the continuity equation for the number density Similarly, an equation for j is found by multiplying (4) with p 1 /m and then integrating over p 1 and ϕ 1 to obtain Here, the two body number density ρ (2) (r 1 , r 2 , t) = dp 1 dp 2 dϕ 1 dϕ 2 f (2) (8) is used in the interaction integral. Now, we effectively make the same two approximations as in [27] to proceed further. The only difference is the dependence of f (1) on ϕ 1 . First, we approximate the interaction integral containing the two body number density via which is exact only in equilibrium. Therefore, the commonly used approximation is to use this expression also in the nonequilibrium case [31]. F exc is the excess free energy which contains energy contributions due to particle interactions. The second approximation is to assume that the ϕ 1 dependencies in f (1) decouple from dependencies on p 1 − p 1 where p 1 = m v(r 1 , t) with the local average particle velocity v that is analogous to the local mean orientation P. To be speci c, we assume that p 1 approximately is Gaussian distributed [32] around p 1 in a way that is independent of ϕ 1 . Therefore, we write the one body density as where Φ(r 1 , ϕ 1 , t) contains the orientational dependence and kT gives the width of the Gaussian distribution and in equilibrium systems would correspond to the equilibrium thermal energy. With this form of the one body density it follows from the de nition of the mean eld densities that Applying now the approximations (9) and (10) to the current equation in (7) results in a dynamical equation for v which only depends on the other mean elds. This calculation is done and described in [27] and therefore not rewritten here. All terms originating from the additional particle orientation vanish up to a coupling term in the momentum equation which does not change the principal calculation but is also found again in the resulting mean velocity equation Diffusion, interactions with other particles and external elds are summarized here under the Helmholtz free energy functional [32] The rst term is the ideal gas free energy with the thermal wave length Λ = h/(2πmkT) 1/2 where h is the Planck constant and the so called excess free energy F exc [ρ] denotes whatever is left after the ideal gas contribution and the contribution due to an external potential V ext (r) have been subtracted from the free energy F [ρ]. Next, we nd a dynamical equation for the mean orientation P by rst multiplying equation (4) withû 1 and then integrating over p 1 and ϕ 1 . The terms of the resulting dynamical equation of P are calculated separately in the following. First, we rewrite the left hand side of equation (4) by using (6) and the expression for the current in equation (11) to obtain Afterward, we simplify the rst term on the right hand side with the approximation (10) to Like the pair interaction in the current dynamics (7) we now approximate the interaction term for the orientation via an equilibrium excess functional dr 2 dp 2 dϕ 2 G (2) (r 1 , leading to the expression for the orientational interaction − dr 2 dp 1 dp 2 dϕ 1 dϕ 2û1 To arrive at the second line, the momentum integrations on the two body density are carried out, partial integration on ϕ 1 is used, and the approximation (16) is inserted. In the last step, the excess functional ] is expressed via the mean orientation P by applying the chain rule for functional differentiation Here, P[Φ] is given by inserting the approximation (5) into the de nition of P in equation (5) which yields All other terms in the resulting dynamical equation of P ever vanish or are dealt with via partial integration. A resulting nonvanishing term is the average one particle torque which describes a particle's behavior of orienting independently of other particles' orientations. Together with equations (6) and (12) we can now write down the nal form of our model for underdamped active systems (1) with the convective derivative D Dt = ∂ ∂t + (v · ∇). The selfpropulsion velocity v 0 = f 0 /α corresponds to the steady state velocity of a particle in an environment with friction constant α and accelerated by the active force f 0 . Although closure relations for particle interactions do not need to be formulated in a functional form this conceptual connection to equilibrium physics proves to be a useful and instructive ansatz in the context of active matter models [3,28]. Equation (21) represents the starting point for many possible mean eld approaches for various active systems with underdamped translational motion. In the next section, we will study one speci c example, where translational and orientational interactions are given by simple free energy approaches.

Example system
By orienting on minimal ingredients for phenomenological models on the microscopic scale [9,23,24,[33][34][35] we now use equation (21) to specify a concrete system by choosing Our model allows density variations which occur preferably in dilute systems. We consider repulsive interactions with an effective nite compressibility of the system described by F where c is the compressibility parameter. The interaction of orientations given by F P has the form of a Swift-Hohenberg functional [36] with the preferred wave number for structure formation q 0 and additional free parameters a, λ, β.
Its physical meaning becomes more evident when inserting equation (22) into (21) resulting in The κ = a − D R and β terms were introduced before by Toner and Tu [37] for the description of ocking in active systems. For κ < 0, the system evolves into the isotropic state v = 0. We restrict to the non-trivial case κ > 0 where the local alignment strength of orientations a outweighs rotational diffusion D R , resulting in a preferred net velocity amplitude κ/β. The second term in the orientation dynamics is interpreted as an anti-alignment interaction with strength λ at a preferred particle distance l 0 /2 = π/q 0 [21,22]. This type of interaction supports the formation of dynamical patterns that are aligned at short distances, but usually do not favor a globally aligned state due to the preference of the anti-alignment on longer distances. Additional one body torques G (1) are not considered. A detailed discussion of the physical background and consequences of the orientational interactions can be found in [21,22]. However, in these works the particles move with an arti cially xed velocity, while the velocity in our approach can vary in amplitude and direction with respect to the orientation of self-propulsion. The orientational interaction in equation (23) corresponds to the one used in [21,22] where results of particle-based simulations are presented. Inserting these orientational interaction rules into our underdamped dynamics equation (1) one arrives at a particle-based model ready for simulations to compare to our corresponding mean eld example in equation (23). We note that in [38], we also use equation (21) to investigate another system with different free energy interactions.
We are especially interested in the consequences of introducing translational inertia into active systems. Therefore, as an important reference case, we determine the overdamped limit of (23) by taking a low-mass limit. The convective derivative of the velocity equation can then be neglected leading to a quasi-stationary velocity dynamics. Additionally, we neglect convection of P which also results from the inclusion of translational inertia. With this, we simplify the dynamics in the overdamped limit to We rescale and numerically implement equations (23) and (24) as described in appendix A. This reduces the set of dimensionless and independent parameters to v 0 , λ for the overdamped and to m, v 0 , λ for the underdamped model. From now on we solely refer to the dimensionless form of these parameters.

Results
We investigate the role of translational inertia by comparing predicted states along different regimes of the anti-alignment strength λ in the under as well as the overdamped model. The in uence of the active drive v 0 and mass m on the state diagram is discussed afterward. We distinguish qualitatively different states from changes of observables for density uctuation and velocity alignment. First, we use the variance of the space averaged density and determine its time average for a time interval T that is chosen long enough such that the average does not depend on T. We use the mean value ∆ρ to quantify density uctuations in space while the corresponding variance is a measure for temporal uctuations. Second, we measure velocity alignment with the space averaged polar orientational order of the normalized velocity eld where ||·|| denotes the Cartesian norm. A value near one indicates global orientational order of the velocity eld while zero indicates the absence of global ordering. In the overdamped model the velocity eld for p v is given from the density dynamics in equation (24) which has the form of a continuity equation. Thus, we have We measure average density uctuation ∆ρ and velocity order p v for varying λ. From the shown results in gure 1 three states can be distinguished within the underdamped model which are discussed later in more detail. First, for high antialignment strengths the average density uctuation is close to the reference value ∆ρ 0 and the temporal variance is negligible. The system shows no global ordering of velocities. In this regime we observe alternating high and low density lanes along which particles move in opposite directions which we refer to as laning [18][19][20]39]. This is the only state which is equally predicted by the overdamped model. Second, near the critical value λ = 0.2, ∆ρ and its temporal variance spontaneously increase which is accompanied by an onset of global orientational velocity ordering in p v . We refer to this transition state as turbulent due to its non-steady character. And third, in the low λ regime ∆ρ continuously decreases again relative to the turbulent state and even below the value of the laning state. The temporal variance stays high. In this regime global ordering of velocities is observed with p v close to one leading to the identi cation of this state as collective motion. We now rst discuss the laning state in the overdamped case in order to compare with the underdamped model afterward.

Laning in the overdamped model
The interaction of particle orientations given by F P in (22) has two contributions. The Toner-Tu local alignment and the anti-alignment at distance l 0 /2. In gure 2(a) we observe that for a high anti-alignment strength λ the system favors a periodic modulation of orientations. Such emergent structures due to competing interactions are for example also observed in Ising models with ferromagnetic and anti-ferromagnetic couplings on different length scales [40]. The resulting orientation eld accumulates density in bands with the same periodicity. This process is balanced by the compressibility of the system which counteracts along the arising density gradients. In the shown steady state the density and orientation elds may be approximated by harmonic modulations around their mean along one spatial direction. From inserting these expressions in the steady state force balance condition we predict the amplitude of the density variation and the corresponding which is in accordance with our numerically found global stripe pattern and its spatial uctuation in gure 1(a). Along the maxima and minima of the density bands the active drive is not balanced by the compressibility of the system and therefore induces particle uxes organized in lanes. Due to the preferred anti-alignment particles in neighboring high and low density lanes propel in opposite directions with the same selfpropulsion velocity v 0 which explains the observed absence of global velocity ordering in gure 1(b). However, we emphasize that the density difference between the opposite propulsion directions leads to a net current in the direction of movement within high density lanes. For high enough active drive we expect all particles might accumulate in the high density lanes and move in the same direction, thereby maximizing the current, which is however not explicitly tested. Such unidirectional laning states are observed in particle simulations with orientational alignment [18]. There, the maximum distance of local alignment coincides with the resulting lane distance. In this reference the formation of lanes is explained as an overreaction of the alignment interaction. This differs from the laning mechanism observed here since alignment of particles happens only locally and an anti-alignment rule dominates at further distances thereby determining the lane spacing. In other active systems laning structures form due to a combination of crowding effects and ocking [41,42]. When the anti-alignment strength λ is lowered the orientation eld becomes more likely to locally form vortices with diameter l 0 /2. Seeing the orientational interactions as the derivative of the vectorial phase eld crystal model (PFC) functional F P the increasing occurrence of vortices between laning domains can be seen as a transition from the stripe to the crystal phase. The difference to a typically used one component PFC interaction [43][44][45] is the coupling of the two vector components of P. So instead of a clear transition from laning (stripe phase) to a regular lattice of vortices (crystal phase) we observe in gure 2(b) that the system is typically stuck in metastable states of laning domains which are separated by defects in the form of vortices in the orientation eld. Those defects increase the functional free energy of the λ term relative to the global laning state as anti-parallel alignment is only locally given.
For the lowest tested λ values the number of vortices in the orientation eld steadily increases until clear laning domains are absent. Instead, the density pattern varies locally depending on the metastable state of the orientation eld. Additionally, we nd several vortices in the orientation eld which arrange in local lattices of alternating clockwise and anticlockwise motion like illustrated in gure 2(c). Different to the hexagonal lattice symmetry arising for one component PFC functionals we solely observe square lattices which enable a frustration free arrangement of the vortices' orientation of circular motion. We nally note that local square lattices are typically distorted by their surrounding leading to inward and outward spiraling vortices in the orientation eld instead of perfectly circular ones. Inward spiraling orientations accumulate density while outward spiraling ones spread density. In the steady state the resulting density variations between contrary rotating vortices are balanced by the compressibility of the system resulting in square lattices of perfectly circular vortices in the velocity eld.

States in the underdamped model
Laning. For high λ values, the laning state found within the overdamped model is equally predicted in the underdamped model. We nd the same value of spatial uctuation ∆ρ 0 and absence of global velocity ordering p v = 0 (see gure 1). We nd vortical defects in the orientation eld for anti-alignment strengths λ 1. However, in contrast to the overdamped limit they, are absent when λ is lowered as exemplarily shown in gure 3(a). We explain this difference with the convection of orientation. A vortex in the orientation eld induces a corresponding circulating ux. Due to their inertia particles radially leave the vortex which gets distorted due to this convection. The anti-alignment interaction counteracts this distortion similar to a centripetal force which holds objects on circular orbits. Consequently, if λ is lowered enough the formation of vortices in the orientation eld is inhibited. In the steady state the grain boundaries between different laning domains then resolve by branching and linkage of lanes with equivalent orientation at the boundary. In this sense, inertial convection is actually benecial for the formation of a global laning state at intermediate anti-alignment strengths since it heals out vortical defects in the orientation eld.
Turbulence. Lowering λ into the turbulent regime suddenly changes the dynamics considerably since now convective ows not only destabilize vortical defects in the orientation eld but also the laning structures due to the weakened anti-alignment. As can be seen from the orientation eld in gure 3(b) the laning structure evolves only locally. The system does no longer reach a steady state but instead large spatio-temporal uctuations are self-sustained over time which is re ected by an increase in density uctuation ∆ρ and its temporal variance. The movie in the supplemental materials gives an impression of the uctuating density eld. The nonsteady character of the turbulent state results from the continuous interplay of local lane formation and their destabilization due to convective ows. The sudden increase of density uctuations coincides with the onset of global velocity ordering in gure 1. Since convective ows destabilize lane formation local alignment is no longer restricted to single lanes resulting in a small but global drift velocity which is re ected in a moderate velocity order p v . Our nding that inertial convection destabilizes the laning state is veri ed by switching off the convective terms in the dynamical equations of velocity and orientation within the underdamped model. Then, we again nd local laning structures and vortical defects like in the overdamped model where otherwise turbulence would arise.
Collective motion. In the turbulent transition state, antialignment is just weak enough so that laning is unstable but still strong enough to inhibit global velocity ordering. This changes when λ is further lowered since then local alignment of orientations increasingly dominates over the weak anti-alignment resulting in the emergence of global orientational velocity order as measured in gure 1(b) and shown in gure 3(c). The homogeneous ow eld transports particle density without accumulating it too much leading to small density uctuations in gure 1(a) comparable to the laning state or even smaller for low enough λ. Since anti-alignment continually perturbs the homogeneous velocity eld the temporal variance of ∆ρ stays on a high level. We suggest that convective ows stabilize the global collective motion as they do not in uence a homogeneous state of aligned orientations but mix any arising misaligned clusters with their surrounding due to convective transport of orientation. Hence, the formation of a larger misaligned cluster is suppressed by convection.

Remarks
Even for a small mass value in the underdamped model the predicted states in gure 3 differ from the ones observed in the overdamped model in gure 2. This seemingly contradictive observation is explained with the convective term in the underdamped orientation dynamics. We neglect this term in the overdamped model to ensure consistency with the corresponding overdamped microscopic equations of motion, which would not yield any convective terms on the continuum level when doing a derivation as in section 2.1. Therefore, we expect the results within the underdamped model to approach those of the overdamped limit model when in the low mass regime also the time scale of convection is negligible compared to those of interactions. This, for example, happens when we decrease the active drive v 0 in the turbulent regime of the underdamped model. Then, the values of the order parameters in gure 1 gradually decrease until the laning state emerges again since typical particle velocities and therefore the impact of the convective terms is reduced. We explored the transition from under-to overdamped dynamics for a different system in another work [38]. There, we neglect convection of the orientation in the underdamped system and nd that its results approach those of the overdamped limit model for m → 0. Furthermore, we note that the value of the particle mass m as a third free parameter is not relevant for the state diagram of the underdamped model. When choosing the mass two orders of magnitude higher as in gure 3, we nd the same qualitative states and the same values for the observables in gure 1 in the long time limit. However, the mass parameter might have an impact on the relaxation dynamics to the steady state, as we have shown for another system [38].
The dynamical equations of our overdamped model (24) and (28) are comparable to those in [21,22]. Since we employ the same interaction rules, we nd a laning state comparable to the one in [22] for predominant anti-alignment. More importantly, for lower anti-alignment strength a transition to collective motion with a turbulent transition regime is observed in [21,22] reminiscent to the ndings within our underdamped model. The authors also explain this with destabilizing convective ows arising from the local alignment. Different to [21,22], in our underdamped model vortex array patterns like in gure 2(c) are absent and instead the laning state is also observed for intermediate λ values. We explain the absence of vortex patterns in our underdamped model due to convection in radial direction away from a vortex center which is enhanced in our model since the actual velocity follows the active propulsion with an inertial delay. The transient patterns in the turbulent regime of [21,22] are described as local vortices or lane like jets in the orientation eld while we solely observe the latter patterns as seen in gure 3(b). This re ects that the transient structures in the turbulent regime correspond to the stable states observed for higher anti-alignment strength. Another mechanism favoring lane structures in our underdamped model over circular vortices arises due do the implementation of repulsion. In [21,22], particles move with constant velocity and actively rotate their orientation to avoid high density regions while in the present model this repulsion, modeled by the compressibility, is a passive one in the sense that it lowers a particle's velocity if it moves to higher densities instead of turning its orientation.

Conclusion
In this work, we have derived the dynamics of a continuum model for underdamped active matter equation (21) based on a dynamical density functional theory for passive systems [27]. Further, we applied the model to an example system by specifying the interaction rules equation (22) on the continuum level which were previously derived in [21] and shown to describe local alignment and distant anti-alignment of orientations.
Finally, we emphasize here that, especially in fully overdamped active systems, any effectively convective effects necessarily arise from interactions or noise while in the present case their physical origin is translational inertia. We stress this difference by comparing the destabilization of lanes in the present underdamped model with active turbulence in overdamped systems. For example, in active nematics the ordered state generally is unstable for non-vanishing activity due to hydrodynamic interactions [3]. There is an effective feedback loop between the nematic relaxing to its aligned state and the thereby induced ows in the background uid which again destabilize this state. In our underdamped system, a similar feedback is present in the turbulent state. The role of the background uid is replaced by inertial convection of the particles themselves which continually destabilizes lane formation. In polar active systems exhibiting active turbulence one term describing hydrodynamic interactions is even structurally equivalent to convection [23,35,46]. Increasing the coupling parameter of this term drives the structured state unstable, leading to the emergence of active turbulence [34].

Appendix A. Rescaling and implementation
We rescale equations (23) and (24) in order to extract physically relevant parameters and for numerical implementation. Independent quantities are rescaled to their dimensionless form according to the rules r → q 0 r, t → κ −1 t, And for the overdamped model in equation (24) we have For numerical implementation, we discretize all elds on a rectangular grid. The simulation box has lengths of multiples of l 0 and employs periodic boundary conditions. A pseudospectral algorithm is used for numerical time iteration. For the evaluation of the convective terms in real space a third order upwind scheme is applied [47]. Time stepping is implemented via a semi-implicit Euler discretization with xed step size. All simulations are started from homogeneous initial conditions. The average density is always set to ρ = 1.