Hidden velocity ordering in dense suspensions of self-propelled disks

Recent studies of the phase diagram for spherical, purely repulsive, active particles established the existence of a transition from a liquid-like to a solid-like phase analogous to the one observed in colloidal systems at thermal equilibrium, in particular in two dimensions an intermediate hexatic phase is observed. Here, we present evidence that the active dense phases (solid, hexatic and liquid) exhibit interesting dynamical anomalies. First, we unveil the growth - with density and activity - of ordered domains where the particles' velocities align in parallel or vortex-like patches. Second, when activity is strong, the spatial distribution of kinetic energy becomes heterogeneous with high energy regions correlated to defects of the crystalline structure. This spatial heterogeneity is accompanied by temporal intermittency, with sudden peaks in the time-series of kinetic energy. The observed dynamical anomalies are not present in a dense equilibrium system and cannot be detected by considering only the structural properties of the system.


I. INTRODUCTION
The dynamics of colloidal particles at high densities has been largely explored, theoretically, numerically and experimentally, in the last decades [1,2]. At thermodynamic equilibrium, the Mermin-Wagner theorem rules out the existence of a crystalline phase in two dimensional systems, characterized by long-range translational order [3,4]. As shown by Halperin, Nelson [5], and Young [6], in two dimensions the melting transition proceeds via two continuous Berezinskii-Kosterlitz-Thouless [7,8] transitions driven by topological defects, i.e., a hexatic-liquid transition, with a quasi-long range orientational order, and a solid-hexatic transition, characterized by quasi-long range translational order and long-range orientational order [6,9]. Density and temperature are the control parameters to move from liquid to hexatic and to solid-like aggregation phases. This scenario has been verified employing dense suspensions of equilibrium colloids [10,11].
Recently, the study of two dimensional systems of selfpropelled particles at high packing fractions has attracted the attention of the active matter community [12,13], since it may offer interesting engineering applications, for instance, in the design of new materials [14]. Although most of the experimental studies have so far focused on the low densities regime, some novel experiments have investigated Janus particles also in the case of very dense suspensions [15]. Another interesting class of high-density non-equilibrium systems is represented by driven granular media where mono-disperse polar grains under shaking [16] display persistent motion. It is also worth to mention some recent studies of artificial mi- * lorenzo.caprini@gssi.it croswimmers at such densities [17], revealing an intriguing experimental scenario for non-equilibrium aggregation. Specimens of active matter systems at very high densities are very interesting also in the biological realm. Typical examples are tissues composed of highly packed eukaryotic cells. Particular attention has been devoted to the dynamics of confluent monolayers [18,19], which slows down as density increases. More recently, an amorphous "solidification" process has been investigated during the process of vertebrate body axis elongation [20], where cells become solid-like.
Theoretical approaches in the statistical physics of active matter focus on simplified models of self-propelled particles, the Active Brownian Particle (ABP) being one of the most studied. Despite the existence of a vast literature concerning the regime of moderate packing fractions, ABP dynamics in the high-density regime is by far less explored. For instance, active crystallization is studied in Ref. [21], where a a shift of the liquid-solid transition line towards larger densities with respect to the Brownian counterpart is revealed. Besides, this transition is accompanied by a true non-equilibrium phenomenon: liquid and solid phases are separated by a region where the suspension is globally ordered but bubbles of "liquid" still persist. A more detailed analysis revealed the occurrence of traveling crystals [22,23], accompanied by the transition to rhombic, square and even lamellar patterns. This phenomenology has been recently confirmed by experiments realized with vibrating granular disks [24].
The construction of the phase diagram of the ABP model [25][26][27], follows the idea that its behavior [28] for small active forces resembles that of passive Brownian particles with the occurrence of "gas-like", "liquid-like", "hexatic-like" and "solid-like" phases, with a shift of the transition lines towards larger densities when activity increases [25,28]. On the contrary, for large self-propulsion (but moderate densities) an unexpected phenomenon oc- curs: the system phase separates even in the absence of attractive interactions [26,[29][30][31][32][33][34][35][36] (the so-called Motility Induced Phase Separation (MIPS) [37][38][39][40]). Differently, at high densities but far from equilibrium, "standard" crystallization seems to occur [21,25]. The general picture suggested by these studies is that the highdensity equilibrium scenario extends, qualitatively identical, to active systems, with the only difference that selfpropulsion may destabilize the ordered phases or induce a phase-separation.
Some other studies suggest that active dense phases in the MIPS region display a richer picture with respect to passive phases in the coexistence region. For instance, a quite different behavior can be found by analyzing the pressure [41,42] and the interfacial tension between the gas and the cluster phase [43,44]. Recently, the study of the particles' velocities has revealed unexpected features which are certainly absent in equilibrium fluids, e.g. the different kinetic temperatures inside and outside a cluster [45], and the spontaneous alignment of velocities in the phase-separated regime [46].

A. Summary of results
Our main findings can be summarised by the introduction of dynamical information about the particles' velocities into the structural phase diagram, as shown in Fig. 1. The phase diagram concerns the two-dimensional high-density regimes (both homogeneous and phaseseparated), with two main control parameters: the per-sistence time, τ , of the active force (which is proportional to the Péclet number and inversely proportional to the rotational diffusivity) and the packing fraction, φ. We recall the definition of φ = N σ 2 π/4L 2 , being N/L 2 the number density and σ the particles' diameter. Our "augmented" phase diagram challenges the widespread idea that the structural properties alone are enough to describe the dense phases of self-propelled particles, and suggests that a richer picture is obtained by including velocity correlations which, in turn, represent an exquisitely off-equilibrium feature of active systems.
Panel (a) of Fig. 1 portrays the phase diagram as a function φ and τ , which reproduces [25], with three homogeneous phases, i.e. the solid phase (S), the hexatic phase (H) and the liquid phase (L), and a nonhomogeneous regime with MIPS-like phase coexistence (P-S), see Sec. II B for details. Our first finding is that the alignment of particles' velocities discovered in [46] in the P-S regime is observed also in the homogeneous liquid, hexatic and solid phases. Particles' velocities form patterns arranging in aligned or vortex-like domains, as shown by the arrows in panels (c) and (d) of Fig. 1, even if the orientation of self-propulsion has negligible order (see color coding in the same panels). The size of aligned domains -quantified by R, defined later, encoded by colors in panel (a) -grows as τ and φ increase, as discussed in detail in Sec. III A. In the homogeneous liquid configurations, the size of the aligned velocity domains is rather small as a result of the absence of translational order, at variance with the solid (denser) case where the sizes reach large values. In the non-homogeneous configurations, only the growth of τ induces the increase of the correlation length as a result of density saturation. We recall that order in the velocity field is absent in any passive Brownian suspension, even at high densities: it is, in fact, a pure non-equilibrium feature due to the presence of propulsion forces in the active dynamics. Interestingly, a similar effect is also observed in fluidized granular materials, but it is caused by the presence of dissipative interactions [47][48][49]. Fig. 1 (b) focuses on the top part (largest densities) of the phase diagram: the color encodes different information here, i.e. the correlation length of the spatial velocity correlation which takes into account also kinetic energy (square modulus of velocity) and not only the orientation of the velocity vectors as in the case of R in panel (a). We observe that increases in the solid phase and saturates in the hexatic or liquid phases, because of the absence of both translational and/or orientational order (it increases again in the phase-separated regime as explained in detail in Sec. III C). In the same panel, a dashed black line delimits a region where heterogeneous spatial distributions and temporal intermittent behaviors of the kinetic energy are observed: interestingly, these anomalies in the kinetic energy field are correlated to the structural defects of the crystal arrangement (see details in Sec. III B).
The article is structured as follows: in Sec. II, we introduce the ABP model for interacting self-propelled particles, summarizing the structural properties of the system. In Sec. III, we present a detailed study of all the dynamical anomalies in the velocity orientation, velocity vector and kinetic energy fields, correlating these anomalies to the different structural properties of the system. A theoretical approach is also presented in Sec. III C that allows us to predict the features of the spatial correlation functions of the velocity field. Section IV is devoted to conclusions and perspectives.

II. THE SYSTEM OF INTERACTING SELF-PROPELLED PARTICLES
We study a system of N interacting ABP disks in two dimensions moving in a fluid at high viscosity (low-Reynolds number). We neglect both hydrodynamic interactions among the particles and inertial terms [14]. The center of mass of each disk, x i , evolves according to the following stochastic differential equation: where γ is the drag coefficient of the fluid and the effect of the thermal noise due to the solvent is assumed to be much smaller than the effect due to the random active force f a i [14]. The term F i represents the force contribution due to steric interactions, such that, Following several studies in the literature [29,46], we choose U (r) as a truncated and shifted Lennard-Jones potential: The constant σ represents the nominal particle diameter while is the typical energy scale of the interaction. For numerical convenience, both these parameters are set to one in the simulations. Even in very packed configurations, each particle can only interact with its first neighbors due to the truncated potential. The self-propulsion force, f a i = γv 0 n i , evolves with ABP dynamics: v 0 is the modulus of the speed induced by the self-propulsion for a force-free particle, n i is a unit vector of components (cos θ i , sin θ i ). The orientational angle, θ i , performs angular diffusive motion described by: where ξ i is a white noise with unit variance and zero average. The constant D r represents the rotational diffusion coefficient and its inverse defines the typical relaxation time, τ = 1/D r , of the active force [50]. We remark that no explicit aligning force is included in the present model at variance with Vicsek-like models where particles' velocities are forced to align to the mean orientation of surrounding particles' velocities [51,52]. Thus, in contrast with Vicsek-like models, the dynamics (1) does not produce any polarization of the directors n i . We also avoid employing any form of self-alignment between the particle velocity and the self-propulsion force responsible for orientation-velocity ordering as recently proposed in [53,54].
A. The effective velocity dynamics of the particle In order to obtain theoretical predictions and interpret the results, following [46], it is useful to switch from the set of variables {x i , f a i } to the transformed variables {x i , v i } eliminating the self-propulsions in favor of the particles' velocities, v i =ẋ i . We underline that the vectors f a i and v i are not aligned, because of the interaction force F i . This is true, in particular, at high densities. The transformed dynamics reads: where both v i and x i belong to the plane xy and r ij = x i − x j . Each term Γ ij is a 2 × 2 matrix, whose elements are: where Latin indices identifying the particles assume values i, j = 1, N while the Greek indices stand for spatial components (x, y). The last addend in Eq. (5), k i , is a noise term which reads: ξ i is the stochastic vector with components (0, 0, ξ i ) normal to plane of motion. At variance with the dynamics of f a i , the modulus of v i is not constant because of the term ∝ ξ i × F i . The result can be easily generalized to the case of finite thermal noise, as shown in Appendix A, using the strategy of [55]. The deterministic part of Eq. (5) represents the dynamics of an underdamped passive particle under the action of a space-dependent friction. Indeed, the first term in the right-hand side of Eq. (5) can be split into a contribution −γΓ ii v i , representing a generalized Stokes force acting on the i-th particle plus a second contribution, −γ N j =i Γ ij (r ij )v j , which depends on particles' relative positions and velocities and vanishes in a passive system.
The ABP equation (5) for the velocity strongly resembles the analogous equation of another schematic model of self-propelled particles, namely the Active Ornstein-Uhlenbeck model (AOUP) [56][57][58][59][60][61][62][63][64][65]. Such a connection has recently been established in some studies [66,67]: the only difference between ABP and AOUP arises from the noise term k i . In the latter case, k i is a simple white noise acting on the velocity as an effective thermal bath. On the other hand, in ABP, k i is perpendicular to n i , the orientation of the active force, and is a multiplicative noise since its amplitude depends both on v i and F i through the unit vector n i . Also notice that (since ξ i has unit variance and n is a unit vector) the variances of the ABP and AOUP noises coincide, being v 2 0 /τ . We underline that, upon fixing v 0 , the noise variance decreases in the large persistence regime.

B. Known results and positional order
In the present Section, we illustrate the phase diagram for the values of the control parameters φ and τ (at fixed v 0 ) explored in this work. Our results are obtained by means of numerical solutions of Eqs. (1) in a square domain of size L under periodic boundary conditions. In the considered interval of packing fractions (φ ∈ [0.78, 1.1]), a suspension of passive Brownian particles exhibits liquid, hexatic and solid phase depending on the temperature. Our ABP phase diagram is consistent with the findings of Di Gregorio et al. [25]. As shown by these authors, as τ is increased the liquid-hexatic and the hexatic-solid boundaries shift to larger values of φ and the hexatic region of the phase diagram is enlarged with respect to the passive equilibrium picture.
When the persistence is smaller than the time-scale associated with the potential, i.e. when τ (∇ · F(x)/γ) −1 , (beingx the average distance between neighboring particles, which is fixed by the density in any homogeneous configurations) we expect the same behav-ior as passive Brownian particles [68]: particles are homogeneously distributed in the box and arranged in the solid, hexatic or liquid phases, as shown in [25], depending on the interplay between φ and τ . In this regime, the active force acts as a thermal bath with effective diffusivity ∼ v 2 0 τ . Thus, the growth of τ , at fixed v 0 , can be mapped in the increase of the effective diffusivity in the corresponding Brownian system. Depending on the interplay between φ and the effective temperature, the system explores liquid, hexatic or solid phases. The structural properties are detected monitoring the behavior of the orientational order parameter [26,69], ψ 6 (x i ), defined as ψ 6 (x i ) = j e 6iαij /N i , where α ij is the angle -with respect to x axis -of the segment joining the i-th and the j-th particle and the sum is restricted to the first neighbors of the particle i, namely N i . In particular, we focus on the correlation function, , to distinguish between different phases. While g 6 (r) is roughly constant with the distance in the solid phase, it decays as an inverse power-law in the hexatic phase. Differently, in the liquid phase, g 6 (r) shows an exponential decay. Examples of the different decays of g 6 (r) in the three structural phases are shown in Fig. 2 (a).
In addition, we also measure the pair correlation function, g(r). At high density, the arrangement of particles is close to the hexagonal lattice, so that the peaks of the g(r) are placed at positionsx, √ 3x, 2x, √ 7x, ... and so on, wherex is the typical distance between neighboring particles, fixed by the density in any homogeneous configurations. It is worthy to note thatx is quite smaller than σ, meaning that particles climb on the interacting potential due to the large densities. Each particle only interacts with its six neighbors due to the truncated nature of U (r). In analogy with an equilibrium system, one can roughly identify the liquid phase with the parameter region where the second peak of g(r) is not split. We see that for all values of the packing fraction, the peaks of the g(r) decrease as τ is increased, in agreement with the interpretation of the self-propulsion in terms of effective diffusivity: fluidization occurs. The g(r) is shown in panels (g) and (h) of Fig. 2 for several values of τ and two densities φ = 0.94, 0.82, respectively: at high density (g), the system maintains a split second peak, while at moderate density (h), the system shows a transition -increasing τ -towards a single second peak where a liquid-like structure occurs [68]. Increasing τ beyond some threshold value (that depends on φ) denoted as a red dotted line in panel (g), the curves representing g(r) saturate, meaning that the internal structure is no more influenced by the value of τ . On the other hand, the shape of g(r) changes towards less fluid configurations where the system shows phase-separation or inhomogeneities. This is not a surprise since, in those cases, particles in the clusters attain a more compact configuration.
The boundary line between the homogeneous and inhomogeneous (phase-separated) regimes is obtained by monitoring the distribution of the local packing fraction,  Fig. 2 for two different configurations at the same densities and different τ . When the system is spatially homogeneous P (φ) displays small tails and a peak located atφ = φ, while in the inhomogeneous region it displays a long tail for smallφ < φ and a shift of the main peak forφ > φ. The line of the transition from homogeneous (L, H, S) to phase-separated (P-H) phases is tracked in Fig. 1 (a), in correspondence of the first point showing such a shift.
Finally, the fraction of defects vs τ is measured for the denser configurations spanning solid and hexatic phases, as illustrated in Fig. 2 (f). A defect is detected by counting the number of neighbors of a particle inside a circular radius of size σ: if the number of neighbors is different from six we mark this point as a defect. The measure is stopped when the system becomes inhomogeneous, in the proximity of colored stars. We observe that the solidhexatic transition takes place where the fraction of defects reaches ∼ 5%. Fig. 1 (c) and (d) are snapshots of the system representing particles' positions and velocities for a large value of the persistence: they do not show MIPS, since at φ = 1.1 the density remains homogeneous in the considered range of τ . In the case of interacting systems, the velocities of the particles, v i , represented by black arrows in Fig.1 (c)-(d), differ from f a i and, in spite of the absence of any alignment interactions, align and selforganize in large oriented domains. On the contrary, the self-propulsion f a i remains randomly oriented. The alignment of velocity orientations corresponds to the collective movements of large domains of particles. Such domains rearrange continuously in time and, sometimes, collapse into vortex structures, at variance with the wellknown traveling bands occurring in the Vicsek-like models [52,[70][71][72].

III. ORDER IN THE VELOCITIES
Hereafter, such a velocity order is studied quantitatively in terms of spatial alignment velocity correlations and suitable order parameters, useful to estimate the size of the domains. We find that there are different aspects in the velocity ordering phenomenology: i) order in the orientation of the velocity vectors ii) order in the full velocity vectors, accounting also for the occurrence of large regions with the same kinetic energy.

A. Velocity orientation
For Vicsek-like models, the global alignment of the particles, also known as polarization, is commonly measured by the following order parameters [51,70,73,74]: where i is the imaginary unit and Θ k is the velocity orientation of the k-th particle. This observable is almost zero for particles without any alignment, typically at low numerical densities and high noise, and returns one, for perfectly aligned particles, e.g. for large values of φ [70]. On the contrary in most systems of swimming active particles, typically evolved by means of overdamped equations, velocity vectors are ignored and the self-propulsion orientation is the only information used to characterize polarization: for instance, in the case of spherical (apolar) ABP particles, the parameter ϕ in Eq. (8) -with Θ k replaced by θ k of Eq. (3) -is close to zero. In this model, since self-propulsions do not coincide with velocities, it is more suitable to consider the orientation of the velocity vector,ẋ k , in Eq. (8), i.e. replacing Θ k with the angle formed by the velocity of the particle with respect to the x axis. However, due to the presence of several domains with different orientations, there is no global velocity-polarization. In principle, for very large persistence we could observe a large oriented domain spanning the whole box, but such a finite size effect occurs only when v 0 τ L, i.e. when the persistence length exceeds the size of the box [22,23]. We do not consider such a case in this manuscript.
A more appropriate indicator, which -even in the absence of a global polarization -gives information about the local alignment of the velocities and its dependence on physical parameters, is the spatial correlation function of the velocity-orientation, Q i (r) [46]. This observable measures the velocity alignment between the particle i and its neighboring particles located in the circular crown of thicknessr = σ and mean radius r = kr, being k an integer positive number, and reads: where the sum runs over the particles within the circular crown defined by the value of k and N k corresponds to the number of particles contained in it. The term d ij is the angular distance between the two angles of the velocities of particles i and j, namely β i and β j , calculated We average over all particles by defining Q(r) = i Q i (r)/N which has the property of being 1 and −1 for perfectly aligned and anti-aligned particles, respectively, and 0 in the absence of any form of alignment. In panels (c-d) of Fig.3, we report Q(r) for different values of the persistence time, τ , and for two different densities. As expected Q(r) is a decreasing function of r. For the smallest values of τ , the alignment is appreciable only in the first shells, a finding consistent with the scenario where the self-propulsion only acts as an effective thermal diffusion. Instead, as τ increases, Q(r) takes on larger values in the first shell and decays slower and slower towards zero, with a typical decay length which roughly represents the average size of one domain: larger values of τ and density (almost always) produce both the increasing of Q(r) in the first shell (k = 1) and the slower decay of the whole function. This observation is, also, qualitatively confirmed by three snapshot configurations obtained for increasing values of τ from panel (e) to panel (f) and (g) of Fig. 3. There, the colors encode the velocity orientations, showing the growth of the average size of the individual velocity domains with τ . The velocity ordering is well captured by the order parameter R, obtained by integrating over the whole box the correlation Q(r): Such a parameter is a measure of the domain size and is studied in Fig. 3 (a), as τ varies for several values of φ, evaluating both the homogeneous and the nonhomogeneous configurations. We recall that the nonhomogeneous regimes correspond to the emergence of empty regions or phase-separation, signaled by the presence of a non-single-peak in the packing fraction distribution (see Fig. 2 (d-e)). For the purpose of evaluating the impact of density inhomogeneity upon R, we show the main peak position,φ in Fig. 3 (b), for different values of τ and average packing fraction φ. We remark thatφ ≈ φ up to some critical value of τ , then it increases. Such a critical value grows as the average packing φ is increased. At the larger packing fraction studied φ = 1.1 the system remains homogeneous for all the explored values of τ . It is worthy to note that the values ofφ (the densities in the denser portion of the system) become independent from φ for τ sufficiently large. As shown in panel (a), R increases with τ . For φ = 1.1, i.e. when the system is spatially homogeneous for all the values of τ , the growth is steady. This proves that the growth of R occurs even in the absence of any local density inhomogeneity since it is not associated with some local change of density. Instead, for smaller values of φ, a first slow monotonic increase is followed by a sharp one occurring at a value of τ for which the homogeneous liquid-like or hexatic phases break down in favor of an inhomogeneous phase. The comparison between panels (a) and (b) also suggests that R has a strong dependence on φ.
The occurrence of such a velocity order is an evidence of the non-equilibrium nature of the dense active phases.
Even if the structural (positional) information suggests an analogy with the liquid, hexatic or solid phases of a passive -equilibrium -system, there is not an equilibrium counterpart of the velocity ordering phenomenon.
In general, it is believed that the occurrence of local velocity alignment is a consequence of the breaking of some microscopic isotropy (as occurs in the case of rod or elongated particles [75][76][77][78]) or the introduction of explicit alignment interactions. The present ABP model subject to random independent active driving leads to velocities' alignment -in the high-density regimes -even for spherical particles. This phenomenology simply arises from the interplay between self-propulsion and steric inter-particle repulsion.

B. Kinetic energy and intermittency
Besides the local order of velocity orientations, spatial correlations also manifest in the speed, v = |v|, and, thus, in the kinetic energy of the particles, ∝ v 2 . Panels (a), (b) and (c) of Fig. 4 show the map of v 2 / v 2 for three snapshot configurations obtained varying τ , for φ = 1.1, i.e. at a density value such that the active system attains a solid-like state for every τ . For small τ , kinetic energies display uncorrelated spatial fluctuations (with Gaussian statistics, not shown here). As τ grows, structures characterized by similar (or correlated) values of v 2 appear, with alternation of fast and slow regions (each identified by a given color).
We also highlight an interesting connection between the kinetic energy spatial distribution and the structural properties of the system. In panels (d), (e) and (f) of Fig. 4, we plot the observable |ψ 6 (x i )| -which is the field pertaining to the crystalline orientational order -relative to the three configurations of panels (a), (b) and (c). The comparison between the maps of |ψ 6 (x i )| and v 2 reveals that the regions with large kinetic energies develop close to the defects of the crystalline structure. A similar scenario occurs for a smaller value of φ, namely φ = 1.02. In this case, the v 2 -map is shown in panels (l), (m) and (n), while panels (o), (p) and (q) report the |ψ 6 |-map. For this choice of φ, the three values of τ distinguish between different aggregation phases: phase-separated, hexatic and solid (from the left to the right). The solid phase for the smaller value of τ is qualitatively indistinguishable from the denser case (compare the panels (c) and (n)). Instead, for the intermediate value of τ (panel m)) the occurrence of the hexatic phase is responsible for a larger number of defects and, thus, a larger number of mobile particles, as clearly shown from the comparison between panels (m) and (p). Finally, in the phase-separated configuration, the fastest regions are mostly concentrated near the boundary of the empty region (panel (l)).
To have another perspective, it is instructive to consider the time behavior of the kinetic energy, v 2 , calculated averaging over a box of size l such thatr L. This observable, as a function of time, is reported for differ- ent values of τ in Fig. 4 (g) and (r) for φ = 1.1, 1.02, respectively. In these two cases, the scenario is similar: for the smaller values of τ the kinetic energy displays symmetric and rapidly uncorrelated fluctuations around the mean value, v 2 . This is coherent with an effective equilibrium picture which is expected when τ → 0. For large values of τ , sparse anomalous peaks manifest, corresponding to rare fluctuations, which move away from their average by several standard deviations. Such peaks become higher and more isolated when τ increases. The observed behavior resembles the temporal intermittency observed in turbulence [79]. A more detailed analysis of such an issue will be presented in a future work, while, in this manuscript, we only consider the essential features of this phenomenon, in particular the role played by defects in the solid phase: particles near the defects attain -in fact -large kinetic energy. To corroborate our observation, we compare the fluctuations of the orientational order parameter, Ψ 6 = N l i=1 ψ 6,j /N l , and those of kinetic energy, v 2 , both averaged over a box of size l with N l particles. In particular, Fig. 4 (h) shows the time-trajectories of 1 − Ψ 6 and v 2 revealing a fair correlation between the occurrence of spikes for both these observables.
With the aim of introducing also the information about intermittency in the general picture, we have drawn a dashed black line in the phase diagram of Fig. 1: the line identifies the intermittency region and is tracked at the first values of τ and φ for which the peaks of the v 2 trajectory overcomes 3 times the standard deviation from its average value.

C. Vectorial velocity field
In the previous Sections, we have seen that the spatial correlation of velocity orientations grows with τ , even in the presence of crystalline defects or large voids (such as those in the phase-separated regimes). Speed (velocity modulus) is more sensitive to the presence of defects and creates patterns with sparse strong fluctuations. A natural question arises: what happens to the spatial correlation of the full velocity vectors (which incorporate both orientation and modulus)?
A quantitative measure of the ordering of the velocities can be obtained by measuring the spatial correlation function in the continuous limit, v i → v(r). v 2 is the variance of the velocity distribution calculated over the whole box. Our analysis limited to the case of homogeneous density under some assumptions is able to predict the form of C(r) (as illustrated in Appendix B). The equation of motion (5) is approximated by the AOUP dynamics, re- placing the multiplicative noise by a two-dimensional additive noise. By this method, it is possible to predict the spatial velocity correlation function. At variance with a similar calculation reported in [46], here we assume that particles are free to oscillate around their equilibrium positions, i.e. they form a hexagonal crystal structure with oscillating sites. Under these simple hypotheses, an expression for C(r) can be derived using the equation of evolution of the velocities. In particular, C(r) displays an exponential-like behavior: where is the correlation length: Thus, a strong potential and/or a large value of τ increases the value of . Also increasing the average packing fraction (i.e. decreasing the lattice constantr) leads to a growth of through the U (r) dependence on this quantity. Expression (13) coincides with the result obtained in [46], even if here it is derived under the less restrictive hypothesis of a vibrating (not rigid) lattice. Moreover, at variance with [46], here we are also considering very high average densities with homogeneous (not phase-separated) configurations, allowing us to directly check the scaling of C(r) with τ . On the contrary, when phase-separation occurs, the packing fractionφ of the dense regions (and therefore the effective value ofr) grows with τ , even at fixed average density φ.
To check the predictions (12) and (13), we study the spatial velocity correlation, C(r), for several values of τ and φ. C(r) for the denser case (φ = 1.1), which corresponds to the solid phase for the whole range of τ numerically explored, is reported in Fig. 5 (c) and reveals a good agreement with theory. The correlation length, , is reported in Fig. 5 (a) (green diamonds) and shows a monotonic growth in fairly good agreement with Eq. (13): under these high-density conditions, the defects are not so statistically relevant and do not interfere with the velocity order. A zoom of Fig. 5, shown in panel (b) and accompanied by a fit of the numerical data, displays a good agreement with formula (13). A similar analysis reveals discrepancies between theory and simulations when applied to lower values of φ: the growth of ceases at some value of τ which depends on φ. We mark with a colored circle the first value of τ where has reached a plateau, and call it τ * (φ). Interestingly, this is close (up to numerical errors) to the value of τ where the solid-hexatic transition takes place (i.e. the fraction of crystalline defects roughly overcomes a given threshold), Here, for completeness, we report the correlation function of Ψ 6 , namely g 6 = Ψ 6 (0)Ψ * 6 (r) , in Fig. 5 (f)-(i). This function shows the well-known transition from the solid to the hexatic phases, roughly at τ * , where g 6 goes from a nearly constant behavior to a power-law decay [69]. For comparison, we also show C(r) as a function of φ for two different values of τ , in panels (d) and (e), where each phase is signed by S, H (solid, hexatic) in the caption. In the hexatic phase, C(r) maintains the exponential-like shape predicted in Eq. (12) but decays abruptly much faster in the proximity of the solid-hexatic transition. As a consequence, the agreement between Eq. (13) and data, shown in Fig. 5 (b), holds up to τ * , while for τ > τ * the presence of defects invalidates the prediction (13) since remains nearly constant or decay very slowly. It is also remarkable that a further increase of τ produces a steep growth of . This is likely caused by phase-separation and the increase of local density in the clustered regions. Even in this case, we still expect that Eq. (13) holds even if the functionx(τ ) is unknown.
In conclusion, we have solid arguments to state that a large number of defects, occurring in the hexatic or liquid phase, is responsible for the saturation of . Indeed, in the proximity of a defect, regions with large kinetic energies are present, as seen in the previous section, and, as a consequence, the velocities are less correlated. Interestingly, the size R of the orientational domains is not so affected by the lack of orientational order, always revealing a monotonic growth with τ independently of the structural phase.

IV. DISCUSSION AND PERSPECTIVES
We have studied systems of self-propelled particles at high packing fractions displaying structural properties which resemble the equilibrium liquid, hexatic and solid phases, exploring both the small and the large persistence regime.
While at small values of the persistence time, τ , many observables behave as in equilibrium, the phenomenology in the high persistence regime is much richer and displays unexpected phenomena. We conclude that an improved active phase diagram benefits from the introduction of new order parameters, related to velocity correlations which have not a Brownian counterpart.
In particular, velocities exhibit intriguing patterns, forming aligned domains or vortex-like arrangements. We propose a suitable order parameter which quantifies the size of these domains, deduced from the spatial correlations of the velocity orientations. Such a parameter grows with both packing fraction, φ, and τ in the homogeneous phases, but it becomes independent from φ in the phase-separated regimes. We also observe the occurrence of large regions whose kinetic energy is highly correlated: these regions become larger when τ and φ grow. Besides, high energetic regions form in the proximity of orientational defects of the solid phase, much more visible in the hexatic phase where defects are diffuse. This is accompanied by a pronounced time intermittency phenomenon, apparently well-correlated to the fluctuations of the orientational order parameter. Correlations of the full velocity vector are also useful to get insights about these dynamical features, in particular, they can be successfully compared to a mesoscopic theory developed under the assumption of homogeneous density. Deviations from the theory appear together with the emergence of a large fraction of defects and the breakdown of homogeneity.
Our observations call for experimental verifications. Promising real platforms to confirm such interesting phenomenologies are Janus particles [15,80] or vibrated polar granular disks [17].
In the present work, we have been interested in understanding the dense phase in monodisperse active systems. Binary mixtures are often employed for gaining insight into active glass phases [61,[81][82][83][84][85][86][87]. In particular, it has been shown that the spatial velocity correlation function is an input ingredient for developing a self-consistent Mode-Coupling Theory of Active Matter [83]. Our findings prove that, as a general result, the statistical properties of the velocity field have to be taken into account for providing a complete description of active systems. As a future direction, it would be interesting to understand how the velocity alignment patterns we observed modify glassy transition.
The phenomenology of domains with aligned velocities could resemble the scenario of traveling crystals, observed in Refs. [22,23] in numerical simulations. Recently, traveling crystals have been observed in experiments using suspensions of micro-disks subjected to vertical vibrations [24]. In such studies, the whole hexagonal pattern moves coherently in space, even though each selfpropulsion vector points randomly. Our phenomenology is quite different since far particles (which belong to different domains) move independently and, thus, the whole crystal gets stuck. Instead, the movement of some clusters gives rise to the formation of defects. We find that the whole crystal travels coherently as in Ref. [22] only if the size of the box is smaller than the persistence length of the active motion.

Velocity correlation function
We, now, consider the real-space velocity correlation function: .
By replacing the lattice sum by a double dimensional integral and defining r = |x − x |, we have v x · v x ≈ 1 2π 2D τ a 2 2 K 0 (r/ ) (B15) where K 0 (r/ ) is the zero-order modified Bessel function of the second kind which has the following asymptotic behavior when r/ 1: which defines the correlation length in the harmonic hexagonal lattice in agreement with the result (13).

Bond angle order and ψ6-field in the harmonic crystal
We define the angle, α x , between the local crystallographic axes and the axes of the ideal lattice [88]: where we used the continuum representation. In Fourier space we have:α q = i 2 (q xûqy − q yûqx ) while theα q correlation reads α qα−q = Dγ 4 where σ 2 = Dγ/K. The real-space α x -correlation function is given by where we have used the expansion for small q and the asymptotic behavior of K 0 (r/l). Now, we consider the correlation function of ψ 6x = e i6αx ψ 6x ψ * 6x = e i6αx e −i6α x . Using the form of the α-correlation we find: Consequently, the correlator of ψ 6 does not vanishes at infinity, i.e. the order is maintained since: This is the expected result since the harmonic lattice always maintains the sixfold coordination number and no disclinations can be created. We notice that for this model the velocity correlation and the α x -correlation have the same long-range behavior.