Hydrodynamics of simple active liquids: the emergence of velocity correlations

We derive the Hydrodynamics for a system of N active, spherical, underdamped particles, interacting through conservative forces. At the microscopic level, we represent the evolution of the particles in terms of the Kramers equation for the probability density distribution of their positions, velocities, and orientations, while at a mesoscopic level we switch to a coarse-grained description introducing an appropriate set of hydrodynamic fields given by the lower-order moments of the distribution. In addition to the usual density and polarization fields, the Hydrodynamics developed in this paper takes into account the velocity and kinetic temperature fields, which are crucial to understanding new aspects of the behavior of active liquids. By imposing a suitable closure of the hydrodynamic moment equations and truncation of the Born-Bogolubov-Green-Kirkwood-Yvon hierarchy, we obtain a closed set of mesoscopic balance equations. At this stage, we focus our interest on the small deviations of the hydrodynamic fields from their averages and apply the methods of the theory of linear hydrodynamic fluctuations. Our treatment sheds light on the peculiar properties of isotropic active liquids and their emergent dynamical collective phenomena, such as the spontaneous alignment of the particle velocities. We predict the existence within the liquid phase of spatial equal-time Ornstein-Zernike-like velocity correlations both for the longitudinal and the transverse modes. At variance with active solids, in active liquids, the correlation length of the transverse velocity fluctuations is sensibly shorter than the length of the longitudinal fluctuations. In particular, the latter depends on the sound speed and increases with the persistence time, while the former displays a weaker dependence on these parameters. Finally, within the same framework, we derive the dynamical structure factors ...


I. INTRODUCTION
In the last decade, significant progress has been made in the study of the collective behavior of active (or also self-propelled) particles, which comprise bacteria, cell assemblies, active colloidal suspensions, vibrated granular particles, autonomous micromotors, bird flocks, etc. [1][2][3]. Understanding how to control and modify their unusual properties is of capital importance in many practical applications, and could revolutionize wide-ranging fields from medicine to robotics. Since one of the possible practical applications of active particles could be self-assembly it would be important to learn how to use them to obtain emergent materials and substances which have reliable, expected, and predictable properties. From a thermodynamic viewpoint, active particles are systems out of equilibrium since they consume energy from the environment or internal chemical processes and generate mechanical persistent motion [4]. In statistical mechanics, this everlasting energy flow corresponds to a violation of the detailed balance condition. Many properties of Active matter are peculiar and absent in passive systems subject only to random thermal fluctuations. Even * lorenzo.caprini@gssi.it in the absence of explicit attractive forces, active particles may exhibit novel types of self-organization: they undergo a type of phase separation known as mobility induced phase separation (MIPS) [5][6][7][8][9], crowd in the proximity of surfaces [10], and form "living crystals" [11,12] that are mobile, break apart and reform again.
Besides the appearance of spontaneous density inhomogeneities, recent experimental, theoretical and numerical investigations have shown the existence of spontaneous equal-time velocity correlations in active matter systems, an emergent collective phenomenon. This new property has been experimentally observed in different contexts, such as cell monolayers [13][14][15][16] and bacterial colonies [17][18][19][20]. Spontaneous velocity correlations represent a peculiar property of active matter systems and, thus, understanding the underlying physical mechanism could be an important advancement. In most cases, these correlations have been explained by invoking velocityaligning interactions which are the basic ingredient to produce the flocking transition in Vicsek-like models and Toner-Tu Hydrodynamics. However, in some systems of self-propelled particles, it is possible to observe fascinating velocity/activity patterns even in the absence of this kind of interactions [21]. Henkes et al. and Caprini et al. [14,22,23] have shown that spatial velocity correlations do appear also in systems of spherical particles without alignment forces: this phenomenon is sim-ply induced by the interplay between persistent active forces and steric repulsion. The first research group compared experimental results observed in systems of highdensity cell monolayers with a phenomenological theory [14], while the second group studied a suspension of two-dimensional active Brownian particles (ABP) under high-packing conditions [22] and put forward a microscopic theory predicting the exponential decay of the spatial velocity correlations. The characteristic coherence length was obtained in terms of the model parameters with [24] or without the effect of the inertia [23]. Recently, Szamel and Flenner [25] investigated active liquids and found numerical evidence of spontaneous velocity correlations. Such a study leads to the conclusion that the emergence of self-organized patterns in the velocity field is a general property of active matter that may occur even in the absence of direct aligning interactions.
The existing theoretical treatments explain this phenomenon employing a microscopic approach where the evolution of the positions and velocities of each particle is explicitly considered. However, a systematic treatment based on coarse-grained collective variables is still lacking. This paper aims to bridge the microscopic and mesoscopic levels by developing a hydrodynamic theory of active liquids able not only to reproduce the observed onset of longitudinal/transverse velocity correlations but also to predict new phenomena such as the slow relaxation of velocity fluctuations at large scales.
The main actors and the logic line of our work are summarised in Fig. 1. Specifically, we consider a system of underdamped interacting active particles characterized by an active force following Active Brownian (ABP) or Active Ornstein-Uhlenbeck (AOUP) dynamics (for recent reviews, see Refs. [26,27] and Ref. [28], respectively). These are minimal models for active matter, which combine mutual conservative interactions (e.g. repulsion for volume exclusion) and Brownian directed motion, but neglect hydrodynamic interactions or other kinds of explicit alignment forces. We first obtain a microscopic description in terms of the stochastic Kramers-Fokker-Planck (KFP) evolution equation of the N -particle distribution function. Then, we connect the microscopic to the macroscopic by defining suitable hydrodynamic variables as averages (over the N -particle distribution) of local single-particle observables. Finally, we derive from the KFP equation the relevant set of hydrodynamic balance equations which represent the first members of an infinite hierarchy involving the velocity and orientation moments of the distribution.
While such a procedure is a standard tool in nonequilibrium statistical mechanics in the framework of passive systems [29], its implementation in the case of active matter offers new interesting challenges [30]: due to the presence of the orientational degrees of freedom associated with the active force, the procedure not only entrains a larger number of fields, hence a different truncation of the hierarchy of equations for the moments, but also requires different approximations to treat the effec-tive force which appears as a result of the interplay of active and repulsive forces. The truncation of the hierarchy can be obtained by various procedures, either by heuristic methods or as in the case of the ABP model more systematically by the use of the Chapman-Enskog method as shown by Steffenoni et al. [31]. This leads to consider the dynamics of at least an extra field, the polarization, in addition to the usual fields of standard (passive) Hydrodynamics. Besides the truncation of the hierarchy, to obtain a closed set of equations, it is necessary to find a suitable representation of the interactions in terms of the hydrodynamic fields. In passive systems, it is well known, as shown in the pioneering theoretical treatment of Irving and Kirkwood (IK) [32], that (purely repulsive) interactions give rise to stress and energy flux contributions and viscous dissipation terms, in the equations for the momentum and energy densities, respectively. For a system of passive particles, the average internal force density can be expressed as the divergence of a stress tensor, providing a connection between statistical mechanics and continuum mechanics. In active systems, the peculiar interplay between repulsive forces and directed motion induces trapping of the particles, reducing their effective motility. Such a mechanism has been represented as an effective interaction that has no counterpart in passive fluids. The understanding of this complex force, which has been the object of systematic investigation in the last two decades [8,10,33,34], is fundamental to obtain the correct hydrodynamic description of the active fluids.
We recall the existence of a rich literature where effective hydrodynamic equations are employed to understand collective phenomena in active systems. There are two main differences with respect to our treatment. First, most of these studies focus on models containing different ingredients, such as explicit alignment interactions (such as in Vicsek-like models) or anisotropies of the particles' shape. Second, the large majority of hydrodynamic studies concerning the ABP neglects the velocity field, typically involving only density and polarization: these fields account for the main macroscopic observed symmetry breaking phenomena, i.e. motility-inducedphase-separation and flocking transition. In many of such theories, the effect of the active force is represented by an effective "active" pressure term [35] or an effective density-dependent local mobility [8].
Neglecting the velocity fluctuations comes from the analogy with equilibrium models and the general argument that, in systems of passive particles dispersed in viscous fluids, such velocity fluctuations decay much faster than the density or polarization fluctuations, typically on a time-scale determined by the viscous drag. From the present study, we learn that this assumption should be reconsidered, as velocity and polarization fluctuations relax on the same slow time-scale. However, from the theories, where the velocity fluctuations are neglected, we learn a few lessons which are crucial to close also our equations. For instance, in some steps of our derivation,

Slowdown of velocity relaxation
Intermediate scattering functions At small k: both decay like polarization (Sections VIII) FIG. 1. A summary of the paper. The panels of the figure illustrate the main logical steps of this work, highlighting the most important observables considered and symbols introduced. They show the essential ingredients of the model (top frames), the hydrodynamic description in terms of coarse-grain fields, the closure approximations (middle frames), and the main results of the paper (last two frames in the third row).
we took inspiration from the recent work of Speck and coworkers [36,37]: they considered a system of overdamped two-dimensional dry active disks (in the absence of hydrodynamic interactions) and by imposing the so-called force closure derived a self-consistent set of equations for the number and polarization densities [38,39]. Their approach captures the interplay between self-propulsion and repulsive interactions through an effective active force proportional to a density-dependent active speed whose form is similar to the one proposed by Cates and Tailleur in their seminal works [40,41] (see here for a review [8]).
Our theory is suitable to describe isotropic active liquids and sheds light on their peculiar properties compared to those of usual passive liquids. A fluctuating linear hydrodynamic approach around the homogenous state [42,43] allows us to derive closed expressions for the spatial correlations of the velocity field. They originate from dynamical mechanisms that have no counterpart in equilibrium liquids, where the velocity field plays a marginal role only determining the relaxation towards the equilibrium. Our theory is coherent with the scenario numerically observed in particle-based simulations both in solid and liquid configurations [14,22,23,25]. A particular outcome of the theory concerns the correlation lengths: we discover that the longitudinal velocity fluctuations are characterized by a correlation length which increases with the persistence time of the active force and is mainly determined by the sound speed. The transverse velocity fluctuations also display some degree of coherence, but the associated correlation length is much shorter, essentially because it depends on the shear, but not on the bulk modulus. We also investigate the spectrum of the hydrodynamic matrix as a function of the wavelength of the perturbation and, finally, we determine the dynamical structure factors and the intermediate scattering functions: from the latter we predict a pronounced slowdown of the velocity field, due to the activity. To the best of our knowledge, the latter point is a central result that has not been highlighted in the literature. At small values of the wavenumber, that is in a range of length-scales relevant for macroscopic behavior, the fluctuation modes of the velocity field of a system of overdamped passive particles (such as colloidal suspensions) decay with a rate dictated by the viscous damping (normalized by mass). As demonstrated here, active liquids display a different behavior: the coupling with the polarization mode induces an important reduction of the velocity field decay rate which is mainly determined by the inverse of the persistence time. At the end of this introduction, we anticipate a simplified discussion of such a general slowdown mechanism, in order to bring to light its essential features.
The paper is organized as follows: in Sec. II, we introduce the stochastic model employed to describe a system of interacting active particles and develop the Born-Bogolubov-Green-Kirkwood-Yvon (BBGKY) hierarchy in the active case, while, in Sec. III, the hydrodynamic equations are derived, providing suitable closures for high-orders moments and the BBGKY hierarchy. Secs. IV and V introduce a linearization procedure around the isotropic steady-state suitable to describe homogeneous active liquids, while Sec.VI employs the method of fluctuating Hydrodynamics to determine, in Sec. VII, the spatial velocity correlation of transverse and longitudinal velocity modes. Finally, Sec. VIII is dedicated to the study of the time-dependent properties of active liquids, showing results for the dynamical structure factor, the intermediate scattering function, and the eigenvalues of the hydrodynamic matrix as a function of the wavelength, both for transverse and longitudinal velocity modes.
A. Summary of results: transverse and longitudinal velocity correlation and slowdown of the velocity field Many hydrodynamic theories for overdamped active particles neglect the velocity field since this is expected to undergo a fast decay likewise in systems of passive colloids in viscous solvents. However, various arguments and analytical results from active particle models suggest that the self-propulsion induces an effective nonequilibrium memory term (effective inertia and spacedependent mobility) with a typical time determined by the persistence of the active force. This description explains the slowdown of the dynamics due to the interactions but seems to be in contrast with the fast relaxation hypothesis of the velocity field assumed in some active hydrodynamic theories. To settle the question, in this paper, we derive coupled macroscopic equations for density n(x, t), velocity u(r, t), polarization p(r, t) and kinetic temperature T (r, t) fields (see Sec. III for definitions and balance equations), starting from a general microscopic active model with only conservative isotropic interactions. In the second part of the paper, the theory is applied to homogeneous active liquids, resorting to linearized fluctuating Hydrodynamics. As illustrated in Fig. 1), a first conclusion concerns the existence of spontaneous order in the velocity field, already found in active solids, but here bearing new features: in particular, the correlation length of longitudinal modes, ξ , is larger than that of transverse modes, ξ ⊥ : and is mainly determined by the sound speed, v s (see Eqs. (82) and (77) of Sec. VII). As a second important result, we show the slowdown of the velocity field: at a large spatial scale and in the overdamped regime, the time decay of both velocity modes (longitudinal and transverse) is of the order of the persistence time, τ , instead of being of the order of the inverse viscosity, 1/γ τ (see Sec. VIII) as in passive suspensions. In the rest of this subsection, we discuss the essential mechanism underlying this slowdown effect. Quite naturally, the active force induces a coupling between the macroscopic velocity and polarization fields. This is clearly shown in Sec. VIII, where dynamical correlations (the so-called intermediate scattering functions) are worked out in full generality both for transverse and longitudinal velocity modes. To give a flavor of the main underlying mechanism, we anticipate the case of zerowavenumber (k → 0) transverse modes. In linear fluctuating Hydrodynamics, the transverse velocity and polarization modes at large scale, V (t) = u ⊥ (k = 0, t) and P (t) = p ⊥ (k = 0, t) respectively, obey the following Langevin equation: where v 0 is the self-propulsion velocity, n 0 is the average number density of the active liquid, (D V , D P ) are effective diffusion coefficients (dependent on the model parameters), and (η V , η P ) are independent white noises with zero average and unitary variance. As illustrated in the rest of the paper, at finite k > 0 the entries of the dynamical matrix and the noise amplitudes depend on k through the transport coefficients of the theory. When v 0 = 0, a liquid of spherical passive colloids is recovered, the polarization becomes irrelevant and the velocity fluctuates around V = 0 with its unique relaxation rate γ. When v 0 > 0, the eigenvalues of the dynamical matrix remain the same (as the determinant and the trace of such a 2 × 2 matrix are independent of v 0 ), therefore one would be tempted to conclude that v 0 does not affect the relaxation of velocity fluctuations. Such a conclusion would be wrong as shown by the exact calculation of the correlation function: where A and B depend on the various parameters, see below, Eq. (92). In general, when γ > 1/τ (thus, in the overdamped regime) the coefficients satisfy A B, therefore the first -rapidly decreasing -exponential becomes negligible in a time of order 1/γ and the decay is dominated by the slower relaxation of the second exponential.
The mechanism summarised here gets more complicated in the case of the velocity longitudinal modes, which in addition to the polarization are coupled to the density fluctuations, but the result is similar: the dominant relaxation rate is 1/τ also for these modes. The nontrivial effect of non-equilibrium coupling between modes, even within a linear modelization, leading to slow relaxations and memory effects, has been studied in the past in the context of fluctuation-dissipation relations [44][45][46][47] and granular materials [48][49][50].

II. THE MODEL
We consider an assembly of identical active particles of mass, m, mutually interacting through a short-range pairwise potential and immersed in a viscous solvent in a d-dimensional space. The solvent is quiescent and not affected by the particles' motion and is assimilated to a heat bath, acting as a source/sink of energy and momentum. Each particle, identified by an index n, is subjected to an external time-independent force, f ex n = f ex (r n ), and a non-gradient force, f a n , simply known as active force. At such a level of description, we assume that f a n can be represented by a stochastic process with memory without specifying the detailed biological/physical mechanisms responsible for the active motion. The active force is expressed as: where v 0 is the typical swim speed induced (in the absence of interactions or external forces) by the active force and γ is the solvent viscosity. The termê n is a d-dimensional unit vector, representing the orientation of the active force. It evolves in time according to an unbiased random process with exponential memory having a typical correlation time, τ . The memory gives to the particle's trajectories the persistence character that is one of the salient features of active motion. In the framework of the continous stochastic processes, the most popular dynamics ofê n are provided by the Active Brownian Particles (ABP) [23,36,[51][52][53][54] and the active Ornstein-Uhlenbeck particle (AOUP) models [55][56][57][58][59][60][61]. The latter has been often used in the theoretical studies (being simpler than the ABP) since it can reproduce the same phenomenology observed with ABP simulations, namely MIPS [62,63] or accumulation near boundaries [64,65]. Unlike Vicsek models or variants, it does not assume any explicit alignment interaction among the active forces of different particles. In the ABP, the orientation is subject to the constraintê 2 = 1 and evolves according to the law: where D r is the rotational diffusion coefficient and ζ n (t) is a white noise vector such that ζ α n (t)ζ β m (s) = 2δ nm δ αβ δ(t − s). Here, latin and greek indices refer to particle numbers and spatial components, respectively. The multiplicative noise in Eq. (4) is treated according to the Stratonovich convention. In the AOUP case, eacĥ e n evolves according to an Ornstein-Uhlenbeck process with correlation time, τ : In the dynamics (5) different Cartesian components of e n are uncorrelated and take on values ranging in the interval [−∞, ∞] and satisfy ê 2 = d. Such a choice guarantees the correspondence with the ABP together with the relation (d−1)D r = 1/τ and a suitable rescaling of the orientationê →ê/ √ d [53,66]. The equations of motion for underdamped active particles of mass, m, with position r n and velocity v n , read [24,67]: where ξ n is a white noise vector with unit variance. Both the friction and ξ n originate from the solvent so that their amplitudes are related by the fluctuation-dissipation relation: where T 0 is the solvent temperature while, again, latin and greek indices refer to particle numbers and spatial components, respectively. The force deterministic term, F n , is due to the interactions with the other particles and is expressed as being U a generic pairwise repulsive potential accounting for volume exclusion effects.

A. Many-particle Fokker-Planck description
It is straightforward to derive the associated Kramers-Fokker-Planck [68] evolution equation for the 9N dimensional phase-space probability density distribution f N = f N ({ê, r, v}, t), where {r, v} indicates a 2 × d × N dimensional phase space point and {ê n } the d × N space of the particle "orientations": where repeated indices are summed. Here, the operator L a accounts for the dynamics of the active force and it is the only term containing the explicit dependence on the model considered (AOUP or ABP). In particular, the operator L a has the following form: where the difference between AOUP and ABP is contained in the form of the d-dimensional matrix D n for the n-th particle. Suppressing the particle index n, for convenience of notation, in the case of AOUP, we have D AOU P = I, that is the identity matrix, while, in the case of two-dimensional ABP, the matrix has the following spatial components [69]: The generalization to the three-dimensional case is straightforward and leads to similar results. We remark that L a is formed by two terms: i) a "deterministic" drag term (proportional to the first derivative with respect tô e) common for ABP and AOUP and ii) a "diffusive" term (proportional to the second derivative with respect toê) that contains the only difference between ABP and AOUP. Moreover, we remark the following property that will be useful later: where the average is performed over all the variables of the system. Therefore, to obtain the mapping between ABP and AOUP it is enough to rescale the ABP orientational vectorê with the constant factor √ d, that is equivalent to map the ABP swim velocity onto v 0 → √ d v 0 . Bearing in mind this minor difference between AOUP and ABP, i.e. rescalingn in the ABP case, the matrix D satisfies the following tensorial relation: In the following, we discuss a theoretical approach that applies to both models since the differences between ABP and AOUP do not lead to major differences in the hydrodynamic description.

B. Reduced description
In order to proceed further as in the passive case, we define the reduced single-particle (marginal) distribution By integrating out the (ê, v, r) coordinates of (N − 1) particles in Eq. (7), we obtain the following non linearequation where the interaction term, Ω = Ω(ê, r, v, t), is defined as: and involves the two-particle distribution function, f 2 = f 2 (r,ê, v; r ,ê , v , t), which is obtained from f N by integrating out (N − 2) particle coordinates. We remark that the relaxation of the system towards the steadystate occurs even in the absence of particle-particle couplings, due to the combined action of the solvent forces and self-propulsion. On the other hand, the interaction Ω not only contributes to the relaxation process via viscous effects, heat transport, and polarization diffusion but is also responsible for the terms describing steric repulsion, effective attraction, and deviation of the local kinetic temperature from the heat bath temperature. Clearly, when the self-propulsion, v 0 , vanishes, the distribution function factorizes into a translational and an orientational part and Eq. (7) can be reduced to the Kramers equation describing the inertial passive colloidal particles in equilibrium with a thermal bath.

III. HYDRODYNAMIC BALANCE EQUATIONS
To handle Eq. (11), which cannot be solved in general even when Ω = 0, we derive the evolution equations for a finite set of moments of the single-particle distribution function (see also [31,70]). The equations for the moments are obtained by multiplying Eq. (11) by a suitable number of products of the velocity, v and ofê and integrating with respect to v,ê the evolution equation for the distribution. In contrast with the hydrodynamical treatment of passive fluids, where only (2 + d) variables are considered, namely, the number and momentum densities, and the kinetic temperature, the state space required to describe active systems is larger and it is necessary to include the local polarization vector (that is the polarization of the active force).
Even in the absence of interactions, the resulting set of balance equations forms part of an infinite hierarchy. As in the case of passive colloids, a truncation scheme is necessary. Moreover, the presence of interactions in Eq. (11), brings in an additional difficulty because the Ωterm contains, through f 2 , the two-particle correlations. Thus, one is faced with two problems: a) break the moments' hierarchy and b) obtain an expression of the interaction in terms of the moments of the single-particle distribution and some known pair correlations.
To develop an active hydrodynamic theoretical description, we introduce the following local fields (depending both on the position r and time t): i) The local density, n = n(r, t): ii) The local velocity, u = u(r, t), iii) The local kinetic temperature, T = T (r, t), where the Boltzmann constant is set equal to 1.
iv) Following the current literature, the local polarization [10] vector, p = p(r, t) is: The fields i), ii) and iii) are those usually accounted for by the Hydrodynamics of passive particles while the field iv) has an active origin. In order to derive an adequate closure of the hydrodynamic equations, it is necessary to write the balance equation for at least two additional tensorial fields constructed combining the components of orientation and velocity vectors: and where the subscript · r indicates that the ensemble averages r still depend on the position of the fluid element.
Here and in the following, the symbol ⊗ stands for the tensorial product, so that for instanceê ⊗ v is the matrix of elementsê jvi . Instead, the symbol : will be used to denote the dyadic product between vectors or matrices. We remark that the tensor ê ⊗ê r is strictly related to the quadrupole tensor employed in the treatment of overdamped active particle systems [35,71], whereas ê⊗v r is a tensor whose trace is proportional to the work performed by the active force.
A. Balance equation for number, momentum, kinetic temperature densities We derive a set of balance equations for number and momentum densities, and kinetic temperature by projecting Eq. (11) onto the Hilbert space spanned by the functions 1, v, m(v−u) 2 /2, respectively, as usual for passive colloids. This procedure immediately leads to the continuity equation by integrating Eq. (11) over velocity and orientation degrees of freedom (in what follows, this procedure will be simply called "integration"): that expresses the density conservation. Here, we have used that dv dê Ω = 0, because the interaction conserves the number of particles. Projecting onto the v-element, i.e. multiplying Eq. (11) by mv and integrating, leads to the momentum balance equation: while the equation for the kinetic temperature field is obtained by multiplying Eq. (11) by m(v − u) 2 /2. In that case, after not shown integrations and manipulations, usual for passive liquids, we obtain: Eqs. (19), (20) and (21) are the minimal set of equations necessary to describe the Hydrodynamics of passive colloids and are expressed in terms of the so-called kinetic contributions and interaction terms. The term P (K) = P (K) (r, t) represents the kinetic contribution to the pressure tensor and reads: while q (K) = q (K) (r, t) the one to the heat flux vector defined as: Finally, Eqs. (20) and (21) contain the interactions terms (simply called B (·) ) that are those involving the two-body distribution through Ω. These terms account for the interactions among particles and are defined as: and In Eq. (20) the second, third, and fourth terms in the l.h.s. represent a momentum flow. In the r.h.s., n f ex and γv 0 p are two sources of momentum due to the presence of the external field and the active force, respectively. Finally, −mγnu represents the average solvent drag force density which opposes the motion in the direction of the velocity, while term −mB (v) is the internal force density.
It represents the rate of change of the momentum of the active particles induced by their mutual interactions. In passive systems, the Irving-Kirkwood theory [32] identifies such a force with the divergence of the pressure tensor, P (the negative of the stress tensor) according to B (v) , which vanishes in the bulk, i.e. under homogeneous conditions, is a well-studied quantity in the case of passive particles. In the case of spherically repulsive interparticle interaction potential and when u = 0, the term B (v) describes the repulsion that is solely due to the inhomogeneous density distribution. When the fluid velocity is inhomogeneous, in addition to the hydrostatic contribution to the pressure, B (v) contains also a viscous contribution [72,73]. The Irving-Kirkwood procedure [32] gives an explicit expression to derive this term. When all terms on the right hand side of Eq. (20) are set to zero, i.e. when the friction, the active force, and the external force are suppressed, such an equation contains the same physics as the Navier-Stokes equation (NSE) describing a compressible viscous fluid in motion. In fact, in the limit of small gradients, the pressure contributions can be expressed in terms of the static pressure and the dynamical stress tensor via the strain rates. The dissipation contained in the NSE is only due to the internal friction of the fluid, while here two additional mechanisms of injection and dissipation of energy are present. The first one, typical of colloidal suspensions, corresponds to the interaction with the solvent, causing momentum suppression and energy reinjection. The second mechanism is the active force which injects momentum through the term γv 0 p but randomizes orientation sinceê follows a stochastic evolution (ABP or AOUP models). The kinetic temperature field T describes how the variance of the velocity distribution depends on the interactions and the active force. The second, third and fourth term in the l.h.s. of Eq. (21) are analogous to those of a passive fluid and represent compressional work and heat transport.
In particular, P (K) and q (K) are the kinetic contributions to the pressure tensor and to the heat flux, respectively. Similarly, the B (vv) terms, representing the rate of change of kinetic energy induced by the interactions, can separated into a collisional heat flux contribution and a term accounting for the viscous and compressional work [29,74]: Such an equation defines the divergence of the field q (C) , the "collisional component" of the heat flux, which in the high-density regime becomes dominant over q (K) . In the energy equation Eq. (21), the effect of the active force is encapsulated in the last parenthesis. This term is proportional to v 0 and has an energetic interpretation: it is nothing but the average work density (performed by the active force) and acts as an additional source of energy, usually larger than T 0 , i.e. the contribution of the solvent.

B. Balance equation for the polarization field and breaking of the moment hierachy
At variance with the Hydrodynamics of passive particles, the momentum balance equation (20) and the temperature equation (21) involve two new fields, the polarization, p and ê · v r , respectively. To fix them we enlarge the set of hydrodynamic equations by projecting the distribution f overê andê ⊗ v and obtain two additional relations. First, multiplying Eq. (11) byê and integrating, we obtain balance equation for the polarization: where we have used that dvdê Ωê = 0. Equation (28) explicitly depends on the two-body cross-correlation between particle velocity and active force v 0 v ⊗ê r and does not depend on the choice of the active force, being the same both for ABP and AOUP models. As shown by Cates & Tailleur, the equation for the polarization brings in the key ingredient responsible for MIPS. Equation (28) reveals that the local polarization changes for two reasons: the advection associated with the incoming flux of particles carrying the active force with them, and a sink term due to the reorientation of the polarization after a typical time τ . The form of Eq. (28) suggests that on a time scale t < τ , p remains approximately constant and the polarization can be expressed as the divergence of the tensor n v ⊗ê r . To determine the latter quantity, we need to consider the polarization flux equation, which is obtained by multiplying Eq. (11) by v ⊗ê and integrating: where the interaction contribution B (êv) = B (êv) (r, t) is defined by the tensor This term represents the highly non-trivial contribution of the interaction operator Ω to theê i v j -moment equation. It describes the combined effect of self-propulsion and steric repulsion and, in recent papers, it has been termed "indirect interaction" [10]. Again, Eq. (29) does not depend on the ABP or AOUP choice. To close the moments hierarchy, we consider the equation for, ê ⊗ê r , simply by multiplying Eq. (11) byê ⊗ê and integrating: By factorizing the averages of the products of three operators in Eq. (31), neglecting the third order cumulant and using the continuity equation, we find the following relation between steady averages: By neglecting the gradient terms in Eq. (32), we get the simpler approximation for ê ⊗ê r : In the AOUP model, Eq. (33) is exact, while in the case of ABP, we have used a further approximation valid in the homogeneous phases, D r ≈ D . Now, Eq. (9) (which holds also for ABP, after a suitable variable rescaling) leads to the result. Further expressions involving gradient terms in the steady-state of ê⊗ê r have been derived for systems of overdamped ABP [75]. These more refined treatments of ê ⊗ê r could be used in our theory, but, here, we have employed the simplest closure since we will consider only homogeneous rather than phase-separated configurations.
To estimate the third-order moment in Eq. (29), we employ a Gaussian ansatz analogous to the one employed to obtain Eq. (32): +n v ⊗ u ⊗ê r + nu ⊗ v ⊗ê r so that the following relation holds: Plugging the approximation (33) in the steady-state version of Eq. (29) and using Eq. (35), we get the following relation for ê ⊗ v r : which is our closure of the moment-hierarchy.
Summarizing the results of this section, we have obtained a set of balance equations for the first moments of the distribution function: where the symbols P and q are a compact notation for P (K) + P (C) and q (K) + q (C) , respectively, and the term ê ⊗ v and ê · v can be obtained by Eq. (37).
To proceed further, we need to develop a suitable closure procedure to express the pressure, heat flux and B (ê,v) as a function of the hydrodynamic moments and their gradients. This task will require the estimate of the various collisional contributions.

C. Closure of the BBGKY hierarchy for active particles
Since Ω, given by Eq. (12), contains the two-particle distribution function, f 2 , it is necessary to have an approximation to determine how it can be expressed in terms of the one-particle and two-particle properties of the system such as the configurational correlation function. For passive systems, there are several approximations concerning f 2 ranging from the Boltzmann Stosszahl-ansatz [76] to the Enskog-like theories [77], to the Kirkwood superposition approximation [78]. They are not exact but, in many cases of interest, lead to fairly good predictions both for hydrodynamic properties. In active systems, the dependence of the distribution functions on the orientation vectorê poses an extra challenge to the effort of finding an accurate approximation for f 2 . However, it is possible to resort to phenomenological arguments to determine the form of B (êv) . The seminal idea was conceived by Tailleur and Cates [8] and elaborated by Speck [39] and more recently by de Pirey et al. [79] within a model of active hard-spheres in infinitely many dimensions. According to Ref. [39], in a homogeneous system of ABP, the conditional probability of finding a second particle a distance r from a first particle fixed at the origin and having a fixed orientationê is axisymmetric with respect toê and gives rise to a force that opposes direct motion. Such an effect due to the interplay between the self-propulsion and the repulsive force can be represented as a force parallel toê but having a renormalized propagation speed: its value changes from v 0 , typical of isolated active particle, to a densitydependent value [80]: where n c is a constant with the dimension of a density that is a function of the model parameters. How does this phenomenological theory translate into the formalism of this paper? Is it possible to choose a particular form of f 2 which is consistent with Eq. (38)? We give the following argument: let us consider first the single-particle distribution f and then proceed to guess the form of the f 2 distribution function. Supposing that the first four hydrodynamic moments are known, we can construct a trial single-particle distribution having these moments by the following ansatz: where φ M = φ M (v, r, t) is the local Maxwellian distribution that is defined by and ψ 0 = ψ 0 (ê) reads: We remark that, in Eq. (39), we have neglected deviations from the local Maxwellian velocity distribution apart from those stemming from the polarization. This assumption can be questioned in the case of phase separation, while it is reasonable in homogeneous configurations. However, still in the latter case, the largest effect of theê field is to distort the distribution function from its Maxwellian form. Now, we turn to the two-particle distribution and, taking inspiration from Eq. (39), we approximate f 2 = f 2 (r,ê, v; r ,ê , v , t) by the following form: where the prime means that the observable is calculated at a phase-space point (r , v ,ê ) that differs from (r, v,ê). The fields ρ = ρ(r, t) and m = m(r, t) are defined as: m ≡ N nê n δ(r − r n (t)) while g 2 = g 2 (r, r , t) represents the configurational pair correlation function.
Once the form of g 2 (|r−r |) is fixed, using the approximation (42) and the specific form of Ω given by Eq. (12), it is possible to obtain closed expressions for the two collisional terms B (v) i and B (vv) . Such a task, normally performed in the study of passive fluids with repulsive interactions, is the subject of a vast literature [77,[81][82][83]. In the following, we shall approximate B (v) i and B (vv) by the corresponding quantities of an elastic hard-disk system. Instead, the evaluation of the B (êv) integral requires the knowledge of orientational correlations which are not normally considered in the studies of simple passive liquids and involves the correlation term proportional tô e · m(r, t)ρ(r , t) in the square parenthesis of Eq. (42). After defining the following marginal correlation: it is possible to rewrite Eq. (30) as where we have integrated over the velocity degrees of freedom. To go further, we need a prescription to evaluate the integral in Eq. (45) containing the pair correlation function G 2 , whose detailed form is unknown. As indicated by our parametrization, G 2 contains a spherically symmetric part g 2 plus a contribution that depends on the angles associated withê,ê and the vector r − r. By symmetry arguments one sees that only theê · m(r, t)ρ(r , t) term contributes to the integral in Eq. (45). To obtain an expression for B (êv) , we approximate the integral in Eq. (45) with the help of a semi-empirical formula that has been proved for active hard-spheres in infinitely many dimensions [79]: Such an equation represents the integral of the interparticle force with the correlation of a particle with polarizationê at r with a second particle at r [39,79]. The formula can be understood as follows: by balancing the active force with the repulsive force we obtain the typical scaling of the potential ∇U (|r − r |) ≈ mγv 0ê . On the other hand, the pair correlation G 2 is very peaked when the two following conditions are satisfied: a) |r−r |/σ ∼ 1 (where the two particles are at the distance of closest approach being σ their effective diameter) and b) (r − r) · (ê −ê) < 0 corresponding to the case of a pair of colliding particles. On the contrary, G 2 is depleted when two particles move apart.
To conclude, by inserting Eq. (46) into Eq. (45), we obtain an approximate but explicit representation of the collisional term: This result is consistent with the force closure employed and recently reviewed by Speck [39], in the simpler case of homogeneous configurations. Notice, that for a passive system, v 0 = 0, this term vanishes being proportional to the active force strength γmv 0 . Moreover, the approximated expression for B (êv) depends quadratically on the local density and is negative. Remarkably, B (êv) is a non-gradient term, resulting from the correlations between particles orientation and the pairwise repulsive force.
Replacing B (êv) given by Eq. (47) in Eq. (37) , we obtain the following self-consistent equation for the elements of the tensor ê ⊗ v r : Such a result has a simple physical interpretation [80]: at finite densities, collisions slow particles down, and reduce the propulsion speed with respect to the value, v 0 , of an isolated particle. At each collision, the combined effects of persistent motion and excluded volume lead to a temporary immobilization of a particle lasting a time τ C ≤ τ . Thus, N C collisions occurring during the propagation time τ at constant speed v 0 reduce the effective distance (i.e. the persistence length) travelled during the persistence time τ , from = v 0 τ to = v 0 (τ − N C τ C ). The number of collisions, N C , is roughly given by / M F , where the mean free path M F ∼ 1/nσ. Therefore, the effective density-dependent propulsion speed becomes v[n] = /τ = v 0 (1 − v 0 στ C n) in agreement with Eq. (38) upon fixing n c ≈ v 0 στ C .

Estimate of the collisional terms B (v) and B (vv)
We approximate the terms B (v) and B (vv) , and the related tensors q and P featuring in Eqs. (37b) and (37c), by the corresponding expressions of a passive suspensions having the same density and temperature. Thus, in the limit of small spatial variations, i.e. in weakly inhomogeneous active liquids, they are represented in terms of gradients of the fields n, T and u. The contributions of pressure tensor and heat flux vector stemming from the direct interactions among the particles, i.e. the ones not explicitly involving the self-propulsion, can be expressed using the standard macroscopic Navier expressions for the heat flux and momentum flux in terms of the hydrostatic pressure, the velocity gradients and the temperature gradient. This is because, at variance with previous approaches, the effect of the active force is not simply recast onto an additional stress contribution (leading to the well-known swim/active pressure [71,84,85]), but its dynamics are still explicitly considered on the same footing as those of velocity and kinetic temperature fields [86]. Hence, in the case of small velocity and temperature gradients, we write: where κ is the thermal conductivity, η the dynamic viscosity, ζ the so-called bulk viscosity and the superscript T denotes the transpose of vectors or matrices. The symbol P h represents the hydrostatic contribution to the pressure and, in equilibrium systems, it can be clearly identified in terms of an ensemble average. In the present case, it can be obtained in principle by evaluating the trace of P (K) + P (C) when the temperature and the velocity fields are uniform. For our scopes, we approximate such a quantity by using the hard-disks equation of state of Henderson [87]: where the effective temperature T * includes the effect of the active force and is defined later in Eq. (59) (see the uniform solution reported in Sec. IV). The remaining contributions to the pressure and heat flux involve the transport coefficients, κ, ζ, η. which were fixed according to the Enskog theory of hard-disks. In the following, it was practical to replace them with the shear kinematic viscosity, the longitudinal kinematic viscosity and the thermal diffusion coefficients defined as: After using the expressions (49) and (50)), the hydrodynamic equations turn out to be: where for compactness in Eq. (53b) we employed P given by (50) . In this way, after eliminating ê·v r and the elements of the tensor ê⊗v r with the help of Eq. (48), the hydrodynamic equations (53a)-(53d) are self-consistent and depend solely on the number, density, temperature, velocity and polarization densities fields. As in the case of the Hydrodynamics of passive liquids, Eqs. (53a)-(53d) are non-linear and apart from a few exceptions, their solutions require numerical tools. For this reason, to obtain some insight, we shall look at their linearized version around a configuration describing a uniform active liquid.

IV. UNIFORM STEADY SOLUTIONS OF THE TRANSPORT EQUATION
The solution of the coupled system of hydrodynamic equations can only be achieved numerically, however, a useful starting point, particularly relevant for our successive linear study, is to consider the stationary homogeneous state, where the system is subject only to a uniform driving force field having introduced the uniform velocity u to ease the notation. Linearizing around this state consists of restricting the validity of our theory to homogeneous active liquids, excluding MIPS from the theoretical analysis. The corresponding FPE for the probabilityf =f (ê, v, t) formally coincides with Eq. (7), upon replacing f ex with Eq. (54) and suppressing the spatial gradient terms in view of our hypothesis of spatial uniformity: whereΩ =Ω(ê, v, t) and L a is still given by Eq. (8).
Notice that the coupling term in the l.h.s. betweenê and v is not symmetric. This setup corresponds to a system of active particles (ABP or AOUP) subject to a constant force field which induces a translation of the center of mass of the fluid at a constant velocity u.
When Ω = 0, i.e. when the interaction is suppressed, the exact steady-state solution of Eq. (55) can be ob-tained in the case of AOUP and reads: (56) where N is a normalization factor and This distribution formally corresponds to that of a free active particle (See, here for a discussion [88]), and is characterized by the following expectation values: Although it is not possible to write the exact uniform solution under the same form as Eq. (56) in the case of interacting particles, one may insist on looking for it as a multivariate Gaussian distribution. This idea is supported by particle-based numerical studies of both ABP and AOUP at a high density which reveal the almost Gaussianity of the single-particle velocity distribution also in regimes of large persistence [69]. In practice, in the homogeneous high-density configurations, the influence of the active force could be recast onto a renormalization of the swim velocity, v 0 → v[n], affecting the work performed by the active force. In this way, we can include the interaction just by modifying the coefficients of the correlation matrix, i.e. Eq. (57d) and Eq. (57e) but not Eq. (57c) (let us observe that the latter is consistent with the closure employed for ê ⊗ê r ). This agrees with the assumption of Cates and Tailleur employed so far and means that Eqs. (57d) and (57e) should be replaced by: The expectation values derived in this section (Eqs. (57a), (57b), (57c) and Eqs. (58b) and (58b), taking f ex = 0) characterize the homogenous steadystate around which we linearize the hydrodynamic equations.

V. LINEARIZED HYDRODYNAMICS
The theory formulated by Landau and Lifshitz in 1957 allows studying the fluctuations of the hydrodynamic fields around their averages [89]. It consists of linearizing the equations for the fields about their bulk values and study their small deviations when stochastic fluxes to the stress tensor and heat flux are added. The amplitudes of the noise terms are determined by the temperature and the transport coefficients of the fluid. In the following, for simplicity, we do not consider the contribution from the stochastic heat flux, but instead, include a stochastic polarization flux to capture the fluctuations of the active force.
First, we resort to a linearization scheme to simplify the structure of Eq. Eq.(48) and Eqs. (53a)-(53d). The linearization is performed around the stationary homogeneous state, where the hydrodynamic fields (according to the results of Sec. IV) take the values n = n 0 , T = T * , u = 0 and p = 0, where we have defined the effective temperature, T * , for notational convenience: We stress that T * takes contributions from both the thermal agitation induced by the surrounding fluid bath and both the active force [90]. As discussed before, the latter includes a phenomenological density-dependence [8].
Since in active colloidal and bacterial suspensions, T 0 mv 2 0 [3] the effective temperature is usually dominated by the active contribution. However, it is useful to keep the thermal contribution, since it yields the correct passive limit when v 0 → 0. We remark that our choice of the steady-state and T * is not arbitrary but, in the case of isotropic homogeneous active fluids, follows from the assumptions made.
To proceed further, we linearize the term ê ⊗ v r in Eq. (48) around the uniform steady-state (Eqs. (57a), (57b), (57c) and (58b). Keeping only terms linear in the fluctuations, Eq. (48) is approximated as: The structure of the hydrodynamic equations suggests separating the contribution of transverse and longitudinal components both for velocity and polarization fields, namely u ⊥ , u and p ⊥ , p , to take advantage of the decomposition into curl-free and divergence-free components. In this way, the linearization about the steadystate (using also Eq. (60)) leads to a couple of equations for u ⊥ and p ⊥ , that are not affected by the other fields: The form of the equations for the longitudinal fields suggests defining the following quantities to simplify the calculations: Taking the divergence of Eqs. (53c) and (53d) and linearizing about the homogeneous state (using also Eq. (60)), in the case of a two-dimensional system (d = 2), we obtain the following set of four coupled equations for n, T , ψ, φ: where we have introduced the sound speed, v s , the thermal expansion coefficient α, and used the following thermodynamic identity: It is convenient to evaluate the two sets of Eqs. (61a)-(61b) and Eqs. (63a)-(63d) in the Fourier reciprocal space. The solutions depend on the phenomenological parameter n c , and the transport coefficients, such as the transverse and longitudinal viscosities ν ⊥ and ν , and the sound speed v 2 s . Additional parameters, such as the thermal expansion coefficient α and the thermal conductivity κ, appear because of the coupling between velocity and temperature fields. Since, to the best of our knowledge, there are no analytical expressions for the transport coefficients of active liquids, we resorted to a drastic but clear choice for their estimate. To fix all these quantities, we employed one of the theories of major success dealing with the hard-disk systems assimilating the steric effects of the interacting active particles to those of two-dimensional disks and treated the contribution of the self-propulsion as an additive effect. Although this choice is somehow arbitrary, it should give reasonable information about the relative importance of the transport coefficients in the expressions for the correlations. As shown in Appendix A, these parameters contain the activity only via T * (Eq. (59)), so that they show a strong dependence on v 0 , while τ and γ play a marginal role. For high density, we expect that the thermal expansion coefficient, α, remains small in such a way that, in Eq. (63c), the coupling term between velocity and temperature fields, that is ∼ αv 2 s is weak. The argument sketched here is described in detail in Appendix B where the explicit expression of α is reported. With this simplification which leads to more clear and transparent results, we assume that the temperature field effectively decouples from the remaining longitudinal modes so that the 4 × 4 linear system given by Eqs. (63a)-(63d), is reduced to a 3 × 3 problem, namely Eqs. (63a), (63c) and (63d) with α = 0.

VI. FLUCTUATING HYDRODYNAMICS
To solve the linear system of partial differential equations for the deviations δa = a − a of the hydrodynamic fields from their homogeneous value, a, we consider their Fourier transforms (denoted by hat-symbols): where we used the following compact notation: The transverse modes of velocity and polarization fields read:û wherek ⊥ is a unit vector such thatk ⊥ ·k = 0. Following the methods of the fluctuating Hydrodynamics, we study the fluctuations of the hydrodynamic fields by considering the stochastic equations obtained by adding suitable noise sources to Eqs. (61a)-(61b) and Eqs. (63a)-(63d). These equations in Fourier representation read: where the matrix M(k) is made up of two blocks and has the following representation: The termξ(k, t) is a noise field that models the contribution of both the stochastic thermal forces due to the interaction with the solvent and the fast degrees of freedom (viz. the higher moments) which are not accounted for in the hydrodynamic description. We fix the noise by requiring that, when the active force vanishes (v 0 = 0), the steady-state velocity fluctuations reduce to those corresponding to an equilibrium system. In other words, we assume that the fluctuation-dissipation theorem holds in the reference passive system. The derivation of this prin-ciple is described in detail in Appendix C both for the longitudinal and the transverse fields. According to this principle the noise has the following properties: whereξ has the same number of components asâ and the effective diffusion matrix D is diagonal with the following non-vanishing elements: The validity of the fluctuation-dissipation theorem only constrains the noise correlations in the equilibrium limit v 0 → 0. In principle, in the active case, the hydrodynamic noise could have different expressions for the diffusion matrix and even non-white time-correlations (for instance if a strong separation of scale is lacking, between fast and slow variables). A systematic coarse-graining procedure that derives rigorously the properties of the noise starting from the full Fokker-Planck equation of the system (Eq. (7)) is typically very difficult and certainly beyond our scope. Historical examples have been obtained in the case of molecular fluids [91][92][93], granular fluids [94], and lattice models [95,96].

VII. EQUAL-TIME CORRELATIONS
After fixing the diffusion matrix with elements D αβ , the equal-time (stationary) correlations in Fourier space, defined as may be determined by solving the Lyapunov equations associated with the dynamics (69): Hereafter, we restrict to the velocity-velocity spatial correlation evaluating separately transverse and longitudinal components as suggested by the block-structure of the matrix M.

A. Transverse equal-time correlations
Considering the lower-right block of the matrix (70), M ⊥ , we determine the Fourier transform of the equaltime velocity-velocity transverse correlations, C ⊥ (k), defined as: For the sake of conciseness, in the following we shall denote by ⊥ and the transverse and longitudinal elements of the tensors. The details of the calculations are reported in Appendix C and lead to the formula: where ξ ⊥ is the correlation length associated with the transverse mode of the spatial velocity correlation, which reads: C ⊥ (k) displays an Ornstein-Zernike form, viz. decays exponentially in real space, a property not having an equilibrium counterpart. Indeed, in the thermal limit, v 0 → 0 or τ → 0, the k-dependence in Eq. (75) disappears and the amplitude of C ⊥ (k) is simply T 0 /m, as expected at equilibrium. When γτ 1, inertial effects induce a similar suppression of the spatial ordering. The spatial velocity correlations become irrelevant when the so-called active temperature is smaller than the solvent temperature, i.e. v 2 0 τ γ T 0 /m (see Eq. (75)). The expression of ξ ⊥ contains two distinct contributions: i) the "thermal" one, which depends on T * and, thus, both on swim velocity and solvent temperature (usually v 2 0 T 0 /m in active colloids); ii) the second one proportional to the transverse viscosity ν ⊥ . According to the estimates of Appendix A, the dependence on τ appears only explicitly both in Eq. (75) and Eq. (76) because ν ⊥ and T * are τ -independent: In the overdamped limit, τ γ 1, ξ ⊥ does not vanish showing that the presence of a finite correlation length is not a consequence of the inertial dynamics. Specifically: where we used Eq. (59) to eliminate the temperature T * and neglected the small contribution T 0 . In the dense regime, we expect that the term proportional to ν ⊥ dominates because the viscosity grows with the packing faster than the first term. Equation (77) has an explicit decreasing dependence on γ: the larger the friction, the smaller ξ ⊥ . Interestingly, ξ ⊥ does not depend on τ (ν ⊥ is τ -independent), a prediction almost in agreement with recent observations made by Szamel et al. who numerically studied the transverse correlation length of active liquids in particle-based simulations (in the overdamped regime) [25]. It is instructive to evaluate Eq. (76) in the opposite inertial limit, τ γ 1, assuming that the active temperature remains larger than the thermal temperature. In this case, the expression for ξ ⊥ reads: In this limit, the term ∝ ν ⊥ in the expression of ξ ⊥ is still τ -independent while the first term scales as ∼ τ . Again, one could expect that the ν ⊥ term is dominant because of its dependence on the packing fraction although particlebased simulations have yet to be performed, in this case.

B. Longitudinal equal-time correlations
The longitudinal spatial velocity correlation, defined by is obtained by solving Eq. (73) considering the reduced 3 × 3 problem, employing the matrix M (k) given by the upper left block of the expression (70). Going back from ψ to u ⊥ , we obtain (see Appendix C): where ξ is the correlation length associated with the longitudinal modes of the velocity field, that reads: Likewise C ⊥ (k), also C (k) has an Ornstein-Zernike form corresponding to an exponential-like decay in real space and becomes k-independent both in the thermal equilibrium limit, v 0 → 0, and in the inertial limit, τ γ 1. The major difference between Eq. (75) and Eq. (80) relies in the dependence of the two correlation lengths ξ ⊥ and ξ on the control parameters. The expression for ξ contains four terms: the first three are similar to those of the expression for ξ ⊥ and depend on T * , v 2 0 and the longitudinal kinematic viscosity, ν ; the last term is determined by the sound speed, v s which has been estimated in Appendix A and that does not depend on τ .
In the overdamped limit, τ γ 1, the longitudinal correlation length assumes a simpler form: While all terms ( provided v[n] > 0) represent positive contributions to ξ , the v 2 s term becomes dominant because it has the fastest increase when the packing fraction becomes large. The coherence length, ξ , displays a strong dependence on the persistence time, scaling as ∼ √ τ , a result in agreement with the numerical results by Szamel et al. [25] based on particle simulations. Interestingly, this prediction is also in accord with the microscopic theory developed for the case of active solids where v 2 s is proportional to the second derivative of the interaction potential calculated at the lattice constant of the solid.
Finally, we discuss the expression for ξ in the underdamped regime, τ γ 1, always assuming that the active temperature is larger than the solvent temperature: lim γτ 1 Equation (83) shows that, in this case, ξ has a faster growth with τ than in the overdamped regime, where ξ ∼ √ τ . The linear scaling obtained in the inertial regime agrees with the numerical and theoretical results obtained by particle-simulations in inertial active solids [24] while further numerical studies are needed to check this result in inertial active liquids.
The above discussion shows that the main difference between liquid and solid (both in overdamped and underdamped regimes) is the presence of a much shorter velocity correlation length of the transverse modes with respect to the longitudinal length, ξ ξ ⊥ . Indeed, the solid supports traveling waves of both transverse and longitudinal type, while through the bulk of a fluid (liquid or gas) only longitudinal waves can propagate. If the medium is not rigid (fluids), the particles will slide past each other and will not generate a transverse wave but only a diffusive momentum propagation (the shear diffusion mode).

VIII. DYNAMICAL CORRELATIONS
The study of the time-dependent correlations reveals some new interesting aspects of the dynamics of the active system at hand. The dynamical structure factors, defined as are obtained by solving Eq. (69), in the frequency domain: with M(k, ω) = iωI + M(k). The vectorsξ(k, ω) and δã(k, ω) are the time-Fourier transform ofξ(k, t) and δâ(k, t), respectively. The latter is defined as: where ω is the frequency and a similar definition holds for the vector of noise,ξ(k, ω), characterized by zero average and the following correlations: with the matrix D introduced in Sec. VII. Multiplying Eq. (85) on the left by M −1 (k, ω) and on the right by δã T (−k, −ω) and averaging over the noise, we obtain the matrix of dynamical structure factors: where in the last equality we have used the Hermitian conjugate of Eq. (85) and the relation (87). Backtransforming from ω to t, we obtain the two-time correlation structure factors in the (k, t) representation, namely the matrix of intermediate scattering functions F(k, t).

A. Hydrodynamic spectrum
Before delving into the discussion of the dynamical structure factor, we consider the dependence of the eigenvalues of M as a function of the wavevector k and the control parameters. In simple liquids, each eigenvalue of the dynamical matrix is associated with a particular fluctuation such as a sound mode or a shear mode. It might be tempting to extend such correspondence to the active case, including two additional eigenvalues associated with the fields p and p ⊥ . However, as the analysis (see Sec.VIII C) of the intermediate scattering functions demonstrates, in active liquids these identifications become unclear, since the relaxation of each hydrodynamic fluctuation is determined by more than one eigenvalue.
We first obtain the eigenvalues associated with the transverse fluctuations from M ⊥ (the lower right 2 × 2 block of M in (70)). They are real and positive for any choice of the parameters: In passive liquids, the eigenvalue λ ⊥ γ (k) describes shear waves and is associated with the diffusion of the component of the momentum orthogonal to the direction of a propagating signal. At variance with simple (inviscid) liquids [78], this eigenvalue does not vanish as k → 0 due to the presence of the solvent drag force proportional to γ and grows as k 2 with a coefficient given by the shear kinematic viscosity, ν ⊥ , that in the present treatment depends on the active force through T * . However, the second eigenvalue, λ ⊥ τ (k), describes a transverse polarization fluctuation that has not a passive counterpart. As we shall see below, this eigenvalue is important to understanding how the velocity field relaxes. We remark that the k 2 -terms both in Eqs. (89a) and (89b) do not increase with τ in the overdamped regime but their amplitudes are mostly determined by v 0 . In the overdamped regime, we have the following inequality: λ ⊥ γ (k) λ ⊥ τ (k) (while in the inertial regime the opposite relation holds), a property that will influence the dynamic properties.
We turn, now, to the longitudinal modes by solving the eigenvalue problem for M (the upper left 3 × 3 block of M in (70)). When k = 0 and v 0 = 0 (passive limit), one can unambiguously identify a λ 0 -mode which describes a density fluctuation, and a λ γ -mode with a momentum fluctuation. In the active case, an additional polarization fluctuation will be described by a third eigenvalue λ τ without a passive counterpart which will be crucial also for the decay of velocity fluctuations. At first, we obtain perturbatively the longitudinal eigenvalues by expressing them as an expansion in powers of k around the k = 0values, namely 0, γ, 1/τ (for λ 0 , λ γ and λ τ respectively). Up to quadratic order, we find the following real and positive expressions: where after the symbol → we have written the dominant contributions in the overdamped regime, τ γ 1. The approximate expressions (90a), (90b) and (90c) are not valid for values of k above a certain threshold where the eigenvalues become complex. To explore these regimes, Fig. 2 shows the three eigenvalues as a function of k for a particular (but relevant) choice of the control parameters, corresponding to an overdamped system (τ γ 1) with a large active force such that, mv 2 0 T 0 . Bearing in mind that positive real parts provide stable modes, i.e. fluctuations that decay in time, while imaginary parts denote the presence of propagating waves, we identify four different regimes labeled with Roman numerals (we do not explore larger values of k because our theory only applies to small gradients). We employ three different colors to identify each eigenvalue and distinguish real and imaginary parts using solid and dashed lines, respectively.
For the smallest values of k (region I), the longitudinal eigenvalues are real and positive and are roughly described by the predictions (90a), (90b) and (90c). Starting from their values at k = 0 (0, γ, 1/τ ), they remain distinct: while the λ 0 -mode represents a diffusive mode (with diffusion constant ≈ v 2 s /γ), both the fluctuations associated with λ γ and λ τ decay in time at a finite rate even when k = 0. The eigenvalue λ γ decreases with k proportionally to −v 2 s /γ, while λ τ decreases as −v 0 v[n]τ 2 γk 2 (see the inset of Fig. 2). Remarkably, the active force affects the longitudinal eigenvalues through the T * dependence in the transport coefficients (as in the transverse case), but the expression for λ τ decreases faster for increasing τ at variance with λ ⊥ τ . This difference reflects the one occurring between ξ and ξ ⊥ in the overdamped case (see Eq. (82) and Eq. (77)).
Region II is different because the λ 0 and λ τ modes (red and green curves) mix and give rise to a pair of underdamped compression/polarization waves propagating in opposite directions (as they have opposite imaginary parts) but decaying at the same rate (see the inset of Fig. 2). As k increases, the real parts of λ 0 and λ τ increase too. The mode associated with λ γ instead maintains its identity and stays real and positive, showing a slow decrease. Region III is again characterized by three distinct real modes: λ 0 and λ γ monotonically decreasing, while λ τ increases until it crosses λ γ . Finally, in region IV, the two modes associated with λ γ and λ τ give birth to a pair of underdamped waves propagating in opposite directions and having a large adsorption rate, while the λ 0 -mode is again real.
It is perhaps useful to recall the scenario in the passive case, v 0 = 0. In this case the λ τ mode is positive, completely decoupled, increases quadratically with k and remains always real. The other two eigenvalues, sayλ 0 andλ γ are real both below a threshold k 1 ≈ γ/(2v s ) and above the value k 2 = 2vs ν 1 − ν /2v 2 s . Between k 1 and k 2 , the modes are complex and represent underdamped compression waves propagating in opposite directions. They result from the coupling between density and the momentum, like sound waves in ordinary liquids, and exist even in the presence of friction, γ > 0. On the other hand, in active systems, the polarizationmode (the one corresponding to λ τ = 1/τ at k = 0) may couple with the other two modes and sustain propagating waves which have a much lower damping rate than the corresponding waves in the reference passive system.

B. Transverse dynamical velocity-velocity structure function
The transverse velocity-velocity dynamical structure factor, S ⊥ (k, ω), can be easily obtained by solving Eq. (88) with M ⊥ (lower-right 2 × 2 block of M): This observable has a simple form being the sum of two Lorentzians each associated with one of the two transverse eigenvalues, λ ⊥ γ (k) and λ ⊥ τ (k) given by Eqs. (89b) and (89a), respectively. Multiplying by e iωt /2π and integrating with respect to ω, we obtain the time-dependent which decays as a linear combination of two timeexponentials, e −λ ⊥ γ |t| and e −λ ⊥ τ |t| , whose relative weight is mainly controlled by the ratio λ ⊥ τ /λ ⊥ γ (because T 0 is usually negligible). The leading term in Eq. (92) depends on the choice of the parameters: in the inertial active regime, τ γ 1, the first term is dominant since λ ⊥ γ λ ⊥ τ while, in the overdamped active regime, τ γ 1, only the exponential with rate λ ⊥ τ survives because λ ⊥ γ λ ⊥ τ . Instead, in passive liquids, where v 0 = 0, F ⊥ (k, t) ∼ e −λγ t . This is shown in Fig. 3 where F ⊥ (k, t) is reported for different values of k and the same choice of parameters as Fig. 2 corresponding to the overdamped regime (in the inset of Fig. 3 the corresponding S ⊥ (k, ω) is plotted as a function of ω). For active liquids (v 2 0 T 0 /m), our plot shows that the relaxation of F ⊥ (k, t) is mainly determined by λ τ , which is ∼ 1/τ in the small-k limit. At variance with passive liquids, a k-dependence is observed in the active case in agreement with the correlation length of the spatial velocity correlation discussed in Sec. VII. We conclude that the intermediate transverse velocity structure function decays more slowly (at a rate ∼ 1/τ ) than the corresponding quantity in the passive case (v 0 = 0) whose decay rate is ∼ γ, as illustrated in Fig. 3.

C. Longitudinal dynamical velocity-velocity structure function
The discussion concerning the longitudinal modes is algebraically more involved than the one regarding the transverse modes, but follows similar lines. By using Eq. (88) with M (the upper left 3 × 3 block of M in (70)), we can express the density-density structure factor, S nn (k, ω), in terms of the longitudinal eigenvalues λ 0 , λ τ , λ γ according to the following formula: where (λ α ) and (λ α ) represent the real and imaginary part, respectively, of the eigenvalue λ α with α = (0, γ, τ ). To gain better insight into the dynamics of the model, we inspect Eq. (93). In regions I and III of Fig. 2, S nn (k, ω), as a function of ω, is a linear combination of three Lorentzians. Instead, in regions II and IV, the dynamic structure factor may in principle develop three distinct peaks at ω = 0 and ω ≈ ±Im (λ C ) where λ C is one of the two complex conjugate eigenvalues in regions II and IV. The width of these peaks is approximately given by |Re(λ C )|. However, when the damping is large (τ γ 1) well-separated peaks are hardly observable because the real part of the eigenvalues is very large.
Let us go back to the main object of our investigation, namely the longitudinal velocity-velocity correla-tion function S (k, ω). This observable can be obtained by Eq. (93) through the general relation: Since this formula depends on a variety of parameters, we shall limit ourselves to illustrate the behavior of S (k, ω) for a special selection as shown in Fig. 4 (a) and (b) for active and passive particles, respectively. We vary ω in correspondence of several values of k and explore the regions represented in Fig. 2. In the active case, after a small-ω regime which depends on k, a maximum in ω is approached and, then each S (k, ω) shows a similar k-independent decay. As k increases, the height of the maximum as a function of ω decreases, and the peak The study of the intermediate scattering functions, F nn (k, t) and F (k, t), provides complementary information on the relaxation behavior of the fluctuations. Let us begin by F nn (k, t): the location in the complex frequency plane of the poles in Eq. (93) indicates that, in region I, this function can be expressed as a linear combination of three different temporal relaxations (each associated with one of the three Lorentzians mentioned above) characterized by exponential decays proportional to e −λ0(k)t , e −λγ (k)t , e −λτ (k)t . As illustrated in Appendix D, the full calculation is performed by time-Fourier transforming the lengthy formula for S nn (k, ω). Again, to proceed further, we restrict to the overdamped set-up analyzed so far. In Fig. 5, we display F nn (k, t) versus t for several values of k. We find that for λ 0 λ τ λ γ and small k, an approximate relation holds: The behavior (95) is derived in Appendix D and displayed in Fig. 5. The function F nn (k, t) is dominated by e −λ0(k)t whose characteristic time diverges diverges as k → 0: as expected, the fluctuations of a conserved density relaxes diffusively towards the steady-state. The relaxation term associated with the polarization fluctuation provides the first correction to the expression for F nn (k, t) with relative weight ∼ λ 0 /λ τ which becomes more relevant as k increases (see Fig. 5). On the contrary, such a contribution is absent in the passive case which is almost k-independent (not shown). By time-Fourier transforming S (k, ω), we compute the velocity-velocity intermediate scattering function, F (k, t), to shed light on the attenuation of the longitudinal velocity modes. Again, the form of S (k, ω) suggests that F (k, t) decays according to three exponential processes: e −λγ (k)t , e −λτ (k)t and e −λ0(k)t but having different weights with respect to those featuring in F nn (k, t) (see Appendix D.) In Fig. 6, F (k, t) is shown as a function of t for several values of k both in the passive (v 0 = 0) and active overdamped cases (such that λ 0 λ τ λ γ ). In these regimes, F (k, t) can be approximated as (see Appendix D): As in the transverse case, F (k, t) shows an initial fast decay mainly due to the e −λγ (k)t contribution. At variance with passive suspensions, where this exponential is dominant and F (k, t) ≈ F (k, 0)e −λγ t , in active suspensions this term has a small relative weight, ∝ λ τ /λ γ . The leading contribution to F (k, t) goes as e −λτ (k)t , but at variance with F ⊥ (k, t), one observes a stronger dependence on k: the larger the wavevector, the faster is the decay of F (k, t), an effect determined by the presence of the slowest exponential e −λ0t of relative weight λ 0 /λ τ and negative sign. As its amplitude vanishes when k → 0, its influence becomes completely negligible in the small-k region but is appreciable when k 10 −2 and is responsible for the negative values of F (k, t) for large times.
In conclusion, the analysis of the correlations using the (k, t) space shows that the slowest modes associated with the eigenvalue λ τ control the long-time behavior of the longitudinal velocity-velocity time correlation.

IX. CONCLUSIONS
In this paper, we have derived a simple hydrodynamic theory for interacting spherical active particles which applies to both ABP and AOUP underdamped systems, described by positions and velocities. We started from a microscopic level where the history of the system requires the knowledge of the full-space Fokker-Planck equation. Then, we coarse-grained the description switching to a hydrodynamic picture by taking the moments of the Fokker-Planck distribution function. While the application of this procedure is standard in many problems, in the case of the active models considered here, the choice of the appropriate hydrodynamic fields is not straightforward since the only conserved variable is the density and the assumption of strong scale-separation for the other modes (with respect to the fast degrees of freedom) cannot be rigorously guaranteed. Nevertheless, such a choice is crucial to observe the collective phenomena of interest. At variance with the majority of previous approaches, we considered the balance equations not only for the density and polarization but also for the momentum and temperature density fields. The interactions were accounted for by including the appropriate collisional terms which determine the so-called Irving-Kirkwood pressure and the transport coefficients of the model but produce also a decrease of motility of the particles, in agreement with previous works.
In the second part of the paper, we investigated how the above equations can be used to predict the fluctuations of the hydrodynamic fields about their average values by employing the methods of linear fluctuating Hydrodynamics. We showed that, by taking into account the momentum field, one can observe, even in the case of large drag coefficients, spatially extended velocity correlations reminiscent of those experimentally and numerically observed in high-density active matter systems. The interplay between persistent active driving and elastic response of the fluid to compression leads to correlations similar to those found in particle-based simulations and theoretical investigations of solid-like configurations of active particles [14,23,24].
The hydrodynamic theory predicts that the longitudinal velocity-velocity correlation function decays exponentially at large distances with a characteristic length that is mainly determined by the sound speed, and is an increasing function of the persistence time of the active force (i.e. the inverse of the rotational diffusion coefficient), but a decreasing function of the solvent drag co-efficient. Surprisingly, we found that the time evolution of the velocity modes at a large scale is dominated by the coupling with polarization and not by the viscous damping, therefore their decay is substantially slower with respect to that of passive liquids. The shear transverse correlation function displays a correlation length that is almost independent of the active force persistence and mainly determined by the transverse viscosity: it is shorter than the longitudinal one because only depends on the shear rate but not on the compressibility of the liquid. Thus, our theory represents an extension to the liquid realm of the analogous treatment of velocity correlations in active solids. Contrary to that case, the present study is not based on the equations of motion for the coordinates of the particles but considers the evolution of the collective hydrodynamic variables. This study shows that the velocity field is an important observable of active liquids, in contrast with the widespread approach where one only considers an effective equation for the density field and thus disregards the velocity. Despite the latter procedure has been successfully employed to understanding the onset of density inhomogeneity, it is legitimate to ask whether a wider picture as the one we presented could shed some new light on the active phase separation. Understanding the possible influence of the velocity ordering process on the structural properties of the system, such as the MIPS transition or the shift of liquid-hexatic and hexatic solid transitions of homogeneous phases [23,[97][98][99], still represents an open question.
where B is the bulk modulus. Explicitly we find: We also need the expressions for the dynamic shear viscosity: η = η 0 1 g 2 (y) + 2y + 1 + 8 π y 2 g 2 (y) (A4) and the bulk viscosity [100] ζ = 16 π y 2 g 2 (y)η 0 , with The formulas contain the pair correlation function at contact distance between two disks, g 2 (y), which in the Verlet-Levesque approximation can be estimated as: Finally, the kinematic viscosity is obtained by dividing η by mn 0 (assuming that the density is constant) Notice that ν 0 diverges as density n 0 → 0, because in the absence of collisions there is no diffusion of momentum, but the phenomenon becomes ballistic. Now, the expressions for ν ⊥ and ν can be obtained as: These explicit expressions of the transport coefficients have been employed in the numerical study. In order to calculate the intermediate scattering functions, we consider the following Fourier transforms obtained by Cauchy's integral formula: ω 2n e iωt (ω 2 + λ 2 0 )(ω 2 + λ 2 γ )(ω 2 + λ 2 τ ) .

Scattering functions in the overdamped regime
It is useful to extract the dominant terms in the expressions for F nn (k, t) and F (k, t) in the regime of parameters reported in the figures to get approximated but explicit expressions. Assuming the condition v 2 0 T 0 /m, we can neglect the term containing R 2 (t) in the expression for F nn (k, t) and the one containing R 4 (t) in the expression for F (k, t). In addition, in the overdamped regime τ γ 1, we have λ 0 λ τ λ γ so that we can simplify the expressions for R 0 (t) and R 2 (t) as follows: where K 0 = K 0 (k) and K 2 = K 2 (k) are two k-dependent amplitudes. In Eqs. (D3) and (D4), after expanding the relative weights of the different terms in powers of λ τ /λ γ 1 and λ 0 /λ τ 1, we have kept only the lead-ing contributions. Using these two expressions, we arrive at the relations leading to Eqs. (95) and (96) F nn (k, t) ≈ N n e −λ0|t| − λ 0 λ τ e −λτ |t| + λ 0 λ 2 τ λ 3 γ e −λγ |t| F (k, t) ≈ N e −λτ |t| − λ 0 λ τ e −λ0|t| − λ τ λ γ e −λγ |t| where N n and N are constants depending on k through K 0 and K 2 and contain the prefactors present in formulas (D1) and (D2). It is remarkable that the leading term in the expression for F (k, t) qualitatively agrees with the prediction obtained from a one-dimensional active solid [102].