Active Matter at high density: velocity distribution and kinetic temperature
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p Active Matter at high density: velocity distribution and kinetic temperature
Lorenzo Caprini ∗ and Umberto Marini Bettolo Marconi Universit´a di Camerino, Dipartimento di Fisica,Via Madonna delle Carceri, I-62032 Camerino, Italy
We consider the solid or hexatic non-equilibrium phases of an interacting two-dimensional sys-tem of Active Brownian Particles at high density and investigate numerically and theoretically theproperties of the velocity distribution function and the associated kinetic temperature. We obtainapproximate analytical predictions for the shape of the velocity distribution and find a transitionfrom a Mexican-hat-like to a Gaussian-like distribution as the persistence time of the active forcechanges from the small to the large persistence regime. Through a detailed numerical and theoreti-cal analysis of the single-particle velocity variance, we report an exact analytical expression for thekinetic temperature of dense spherical self-propelled particles that holds also in the non-equilibriumregimes with large persistence times and discuss its range of validity.
I. INTRODUCTION
Active Matter at high density is becoming a subject ofgreat interest since it plays a crucial role to understand abroad range of biological systems [1–4], such as cell mono-layers and living tissues. Experimental observations real-ized with cell monolayers reveal large-scale collective mo-tion, like swirls and velocity alignment [5–9]. In the spiritof minimal modeling, these systems have been recentlymodeled using high-density interacting Active BrownianParticles (ABP) [5, 10], thus, modeling the complex cellsand cell-substrate interactions through steric interactionand self-propulsion. Depending on their density, inter-acting systems of ABP display a variegate phenomenol-ogy. In particular, at moderate packing fractions, a non-equilibrium phase-coexistence known as Motility-InducePhase-Separation (MIPS) [11–15] occurs even in the ab-sence of attractive interactions [16–25]. Depending onthe active force and the packing fraction, ABP can alsoattain homogeneous configurations, such as active liquid,hexatic and solid phases [26–32]. With respect to equi-librium systems of Brownian colloids, the active liquid-hexatic and hexatic-solid transitions are shifted towardslarger values of the packing fractions and the hexaticphase occurs for a broad range of parameters [30, 33–35]. Moreover, the dense phases of Interacting ABP dis-play a plethora of dynamical phenomena making themquite different as compared to passive dense phases. Inparticular, the particle velocities spontaneously form or-dered domains even in the absence of explicit alignmentinteractions in phase-separated configurations [16] and inactive liquid, hexatic and solid phases [35] and give rise tofascinating intermittency phenomena [35, 36]. The spon-taneous alignment mechanism makes the ABP modelssuitable to describe the behavior of cell monolayers.For systems of interacting ABP, several authors,searching for an extension of equilibrium thermodynamicconcepts introduced an effective “temperature” in thestudy of non-equilibrium systems of self-propelled par- ∗ [email protected] ticles. Following ideas from glassy systems, several au-thors introduce the dynamical effective temperature bythe ratio between the mean-square-displacement and thetime-integrated linear response function due to small per-turbation, in the context of active disks [37–46], dumb-ells [47, 48], polymers [49], looking both at active ho-mogeneous (liquid, hexatic and solid) configurations andphase-separated regimes. In the homogeneous case, theeffective temperature increases as the propulsive speedincreases and decreases as the packing fraction grows. In-stead in the inhomogeneous case, i.e. when MIPS occurs,the net distinction between the populations of the two co-existent phases, e.g. slow particles in the dense clustersand fast particles in the disordered phase, allows us to in-troduce two distinct effective temperatures [46]. Mandalet. al. [50] focused on the kinetic temperature, i.e. thevariance of the velocity distribution, for underdampedself-propelled particles, outlining that the temperaturesin the two coexisting phases of MIPS are different. Analternative definition of the active temperature has beenalso proposed in the context of stochastic thermodynam-ics to generalize the Clausius relation to active systems, aprogram requiring the introduction of a space-dependenttemperature that depends on the potential itself [51–53].In this work, we shall not discuss the concept of tem-perature in non-equilibrium active systems [54] intendedas an observable satisfying well defined thermodynamicrelations, an issue still matter of debate, but focus onthe kinetic temperature of self-propelled particles. Wefind its exact analytical expression as a function of themodel parameters for dense homogeneous configurationsin the non-equilibrium active solid and hexatic phases.We also obtain the ABP single-particle velocity distri-bution in these highly packed configurations as the per-sistence time of the self-propulsion varies. The shape ofthis distribution obtained via numerical simulations iscompared with theoretical predictions both in the smalland the large persistence regimes. We find a crossoverbetween the two regimes which manifests itself in a qual-itative change in the shape of the velocity distribution.The manuscript is organized as follows: in Sec. II,we introduce the model, while, in Sec. III, we reportthe velocity dynamics representing the starting point ofour theoretical approach. Numerical and theoretical re-sults of single-particle velocity distributions are shown inSec. IV, while the analysis of the first moments of thedistribution and the discussion about the kinetic tem-perature are reported in Sec. V. Finally, we report somediscussions and conclusions in the final section. II. INTERACTING SELF-PROPELLEDPARTICLES
We consider a two-dimensional system of N self-propelled disks, described by the Active Brownian Parti-cles (ABP) model, where inertial and hydrodynamic ef-fects are neglected. The position, x i , of each disk evolvesby the following stochastic differential equation: γ ˙ x i = F i + f ai . (1)The constant γ is the solvent friction while we do nottake into account the thermal diffusivity since for severalexperimental active particle systems [1] is some ordersof magnitude smaller than the diffusivity associated withthe self-propulsion force, f ai . According to the popularABP model, f ai is a time-dependent force given by theequation: f ai = γv n i , (2)where v is the constant modulus of the swim velocityinduced by f ai and n i is the orientation vector of com-ponents (cos θ i , sin θ i ) evolving through a stochastic pro-cess. In particular, the orientational angle, θ i , performsangular diffusion: ˙ θ i = p D r ξ i , (3)where ξ i is a white noise with unit variance and zeroaverage and D r is the rotational diffusion coefficient. Weremark that the inverse of D r defines the correlation-timeof the active force, namely τ = 1 /D r [55], which will beassumed as a control parameter in the numerical studyperformed in this manuscript.The term F i represents the repulsive force betweenparticles due to steric interactions. In particular, F i = −∇ i U tot where the potential, U tot , can be expressed as U tot = P i
4, to the value 1 .
1, where the system attainsactive solid or hexatic configurations without showingdensity inhomogeneities [35]. In particular, the hexatic-solid transition is controlled by τ and occurs approxi-matively at τ = 0 .
1. Under these conditions, we studythe single-particle velocity distribution varying τ and itsmoments. We can distinguish between two regimes [35]:i) the small persistence regime where τ < U ′′ (¯ r ) /γ andii) the large persistence regime where τ > U ′′ (¯ r ) /γ , be-ing ¯ r the average distance between neighboring particlesthat is fixed by the density in any homogeneous con-figurations. In the case i), the self-propulsion f ai is thefastest degree of freedom: in this regime, the persistencetime, τ , is smaller than the typical time of the potential U ′′ (¯ r ) /γ , so that the behavior of ABP resembles that ofpassive Brownian particles and the x i just display oscil-lations around their equilibrium positions. Consideringthe structural properties of the system, this regime is in-distinguishable from the passive solid-state. In case ii),the evolution of f ai plays a relevant role and affects thedynamics of x i , manifesting itself in several dynamicalanomalies [16, 35] due to the intrinsic non-equilibriumnature of active models. III. THE VELOCITY DYNAMICS
As already reported in [16, 35, 57], the study of thevelocity dynamics reveals the existence of hidden collec-tive behavior of self-propelled particles at high densityin the regime of large persistence times. Nevertheless,many single-particle properties, such as the velocity dis-tribution and its moments, have not been yet explored.Following [16], we eliminate f ai in favor of v i = ˙ x i ,i.e. the velocity of the particle, which does not coincidewith the swim velocity, v n i , since the modulus of v i isnot fixed and its orientation is not parallel to n i . Thisstatement is true when particles interact and, thus, athigh densities, in particular. Transforming the dynamicsfrom the variables ( x i , f ai ) to the new variables ( x i , v i )(without any approximations), the equations of motionread: ˙ x i = v i (5a) τ γ ˙ v i = − γ N X j =1 Γ ij ( x i − x j ) v j + F i + τ γ k i (5b)where each Γ ij is two-dimensional matrix with compo-nents Γ αβij ( r ij ) = δ ij δ αβ + τγ ∇ iα ∇ jβ U ( | r ij | ) . (6)Greek indices are used to denote the spatial components α, β = x, y while Latin indices identify the particle num-ber i, j = 1 , ..., N . Finally, the term k i is a noise vectorthat reads: k i = v r τ ξ i × γ v i + ∇ i U tot γv , (7)where ξ i is a vector with components (0 , , ξ i ) and nor-mal to the plane of motion, ( x, y, k i isa multiplicative noise depending both on v i and F i andis perpendicular to n i , i.e. the orientation of the activeforce. Its amplitude scales simply as ∼ v p /τ since n i is a unit vector.The dynamics (5b) resembles the evolution of under-damped passive particles which are out from equilib-rium because of the occurrence of space-dependent fric-tion forces (e.g. the diagonal terms of the matrix Γ )and effective forces depending on positions and veloci-ties of neighboring particles (e.g. the non-diagonal termsof Γ ). Eq. (5b) resembles the dynamics of the ActiveOrnstein-Uhlenbeck particle (AOUP) [58–66], an alter-native model used to study the behavior of self-propelledparticles. Upon a suitable mapping of the self-propulsionparameters [67, 68], the difference between AOUP andABP dynamics is represented by the noise term k i , whichin the former is a white noise vector with independentcomponents [58, 59]. IV. PROBABILITY DISTRIBUTIONFUNCTION OF THE VELOCITY
We numerically study the distribution of the velocityin the steady-state to evaluate the effect of the persis-tent time, τ . As illustrated in Fig. 1 (a) and (b), thedistinction between large and small persistence regimesproduces different shapes in the probability distributionfunction of the velocity, p ( v x , v y ). In the small- τ con-figurations shown in panel (a) [case i)], p ( v x , v y ) has apronounced non-Gaussian shape: the probability of find-ing a particle with v ≈ ≈ v . Instead, in the large- τ configurations reportedin panel (b) [case ii)], p ( v x , v y ) presents a Gaussian-likeshape quite similar to the case of passive Brownian par-ticles. Intuitively, in case i), the self-propulsion changesrapidly without producing any appreciable change in theparticle positions giving rise only to very small fluctua-tions. The net steric force exerted by the neighboringparticles on a tagged particle almost cancels out and haslittle effect on the particle velocity, that in practice onlyexperiences the influence of the active force. This ex-plains why the distribution of v is very similar to theone of f a . On the contrary, in the large τ regime, this isno longer true. The direction of the active force beforeappreciably changing need an interval ∼ τ much largerthan the relaxation time associated with the interparticlepotential, τ p = U ′′ (¯ r ) /γ . Now, the resultant of the activeand steric forces nearly vanishes and (see Eq.(1)), as aconsequence, the average velocity is almost zero and thesingle-particle kinetic energy decreases. The transition from regime i) to regime ii) is quanti-tatively evaluated in Fig. 1 (c)-(d) where the marginalprobability distribution of the velocity along one compo-nent, p ( v x ) = R dv y p ( v x , v y ), is studied for different val-ues of τ . For smaller values of τ , p ( v x ) displays two sym-metric peaks near v x ≈ v . When τ grows, the peaks shifttowards smaller values of v x and their heights decreasewith respect to the p ( v x ) value around the origin. For τ & × − , the two peaks merge and the bimodality ofthe distribution is suppressed, while for further τ -values asingle pronounced peak placed at v x ≈ p ( v x ) and a Gaussian distribution is notvery pronounced as revealed in Fig. 1 (d), as shown by thecomparison with the best Gaussian fit. To provide an-other perspective, we study the probability distributionof the velocity modulus, p ( v = | v | ), in Fig. 1 (e)-(f). Inthe small- τ regime (panel (f)), the distribution is peakedaround v = v displaying a quite symmetric shape fairlydescribed by a Gaussian centered in v = v : p ( v ) ≈ N v exp h − α v − v ) i (8)where N is a normalization factor and α a parameterthat satisfies, h v i = h v i + 3 α/
2. We observe that p ( v )becomes narrow as τ grows, in the small- τ regime, but,for further values of τ , the peak of the distribution shiftstowards smaller values. After a crossover regime, occur-ring for intermediate values of τ , the distribution p ( v )approaches a Gaussian-like shape such that p ( v ) ≈ v exp ( − v /β ) , (9)where 1 /β = h v i , consistently with the observations ofpanel (d). Interestingly, in panel (f), we show that the v -distribution collapses for v → vτ / for a large range of τ between (10 − , − ). Thus, 1 /β , which plays the roleof an effective temperature, decreases as τ is enlarged.We remark that the scaling of p ( v ) ceases to hold for val-ues τ & − when the homogeneous active solid phasebreaks down. These observations will be clarified in thenext theoretical sections. A. Theoretical predictions
From the set of stochastic equations (5), we derive theFokker-Planck equation for the probability distributionfunction, p = p ( { x } , { v } ) (where the symbol {·} has beenintroduced to denote all the space-components of the N particles): ∂∂t p = − v i · ∇ x i p + 1 τ (cid:18) I ij + τγ ∇ x i ∇ x j U (cid:19) ∇ v i · ( v j p )+ ∇ x i Uτ γ · ∇ v i p + v τ ∇ v i ∇ v j ( D ij p ) , (10) FIG. 1. Probability distribution of the velocity. Panels (a) and (b): map of the two-dimensional probability distributionfunction, p ( v x , v y ) for two different values of τ = 10 − (panel (a)) and τ = 10 − (panel (b)). Panels (c) and (d): marginalprobability distribution function, p ( v x ), for several values of τ (colored lines). The dashed black lines in panel (d) are obtainedby numerical fits, obtained with Gaussian distributions. Panels (e) and (f) report the probability distribution of the velocitymodulus, p ( | v | = v ), for different values of τ . In particular, in panel (f), we show p ( v ) rescaled by τ / . Panels (c), (e), and(d), (f) share the captions. Numerical simulations are obtained with v = 50, ǫ = σ = 1 = γ = 1. where I ij is the identity matrix and each element D ij isa 2 × D ij = δ ij n y i − n x i n y i − n x i n y i n x i
00 0 0 . We remark that each D ij is a non-diagonal matrix as aconsequence of the complex noise structure in Eq. (5) andof the fact that n x ( y ) is a function of the particle velocityand position through Eq. (1). We point out that Eq. (10)has the same form as the AOUP Fokker-Planck equationexcept for the diffusion-like term (i.e. the term containing D ij in Eq. (10)). However, in the AOUP equation, thenon-diagonal matrix D ij is replaced by a diagonal one,˜ D ij , with components˜ D ij = δ ij . We observe that ˜ D ij can be obtained from D ij just byreplacing n x ( y ) and n x n y by their averages, i.e. h n x ( y ) i =1 / h n x n y i = 0, respectively (with the addition of anextra factor 2 needed for consistency between the param-eters of the two models [67]). We remark that, even in thesimplified AOUP case, the solutions of Eq. (10) for τ > τ γ around a Gaussian distribution [59, 69]. On the contrary, in the case of interacting ABP, there are neitherasymptotic nor approximated results for the probabilitydistribution function of the velocity.Being the general solution of Eq. (10) unknown, wewill employ suitable approximations supported by nu-merical observations. In Fig. 2, we report h| F |i and h| v |i as a function of τ . In the small-persistence regime, v ≈ h| v |i ≫ h| F |i while in the large-persistence regimethe opposite relation holds, namely v ≈ h| F |i ≫ h| v |i ,confirming the physical explanation mentioned before.Such observation will be crucial in the following to deriveapproximate analytical solutions of p ( v ).
1. Large persistence regime
Using the observation v ≈ h| F |i ≫ h| v |i , holding inthe large persistence regime, we have n i ≈ ∇ i U , and thediffusive matrix can be approximated as follows: D ij ≈ δ ij v γ ( ∇ y i U ) −∇ y i U ∇ x i U −∇ y i U ∇ x i U ( ∇ x i U )
00 0 0 . In the active solid phases, where the defects of the crys-talline arrangement are negligible, the sum of the forcesexerted by neighboring particles cancel out and we canassume ∇ x i ( y i ) U = 0. In other words, our approxima-tion consists in replacing ∇ x i ( y i ) U with h∇ x i ( y i ) U i = 0.This is true only in the high-density regime where we can FIG. 2. h| v |i (yellow data) and h| F |i (green data) as a func-tion of τ . The black dashed line is drawn in correspondenceof v that is chosen as v = 50 in the numerical simulations.The other parameters are ǫ = σ = γ = 1. approximate the gradient of the potential expanding thedistance between neighboring particles around ¯ r . Since h∇ x i U ∇ y i U i = h∇ x i U ih∇ y i U i = 0 and h∇ x i U ∇ x i U i = v , we obtain that D ij ≈ ˜ D ij , proving that, at very highdensity in the large persistence regime, the ABP velocitydynamics is well approximated by the AOUP dynamics.Even with this simplification, an exact solution is notknown, and we shall employ an approximation for themany-body velocity distribution [58]: p ( { v }|{ x } ) ∝ exp − v X ij v i · Γ ij · v j . (11) p ( { v }|{ x } ) is a multivariate Gaussian coupling the wholeset of velocities through the space dependent matrix Γ ij .We remark that the prediction (11) is not the exact so-lution of the Fokker-Planck equation associated with theAOUP interacting dynamics, but it is a suitable approx-imation which works in the large persistent regime. Thisprediction has been tested in several cases, even underthe action of external potentials [52, 70].In the active solid phase, the matrix Γ ij , which de-pends on particles’ relative positions, simplifies due tothe hexagonal structure and to the short-range nature ofthe interaction potential. Thus, the conditional proba-bility distribution of v i (i.e. knowing the velocity of theother particles) is given by Eq. (11) where the sum isrestricted to the first six neighbors of the target parti-cle. Integrating out all the velocity degrees of freedomexcept v i , we still obtain a Gaussian distribution withzero average confirming the shape reported in Fig. (1) inthe small persistence regime: p ( v ) ∝ exp (cid:18) − β v (cid:19) (12)where β is the variance of the distribution or the inverseof the kinetic temperature. Its exact expression as a func- tion on the active force parameters will be derived inSec. V.
2. Small persistence-regime
In the small persistence regime, the AOUP approxima-tion is no longer valid, as numerically shown in Fig. (1).Indeed, according to the AOUP model, the shape of p ( v x , v y ) should always be Gaussian (with asymptoticcorrections) at variance with our numerical results ob-tained with ABP simulations. As shown in Fig. 2, wesimplify the noise matrix assuming v ≈ h| v |i ≫ h| F |i ,obtaining D ij ≈ δ ij v v y − v y v x − v y v x v x
00 0 0 . (13)The Fokker-Planck equation (10) with the matrix (13)(in the small τ limit) turns to be ∂∂t p ≈ − v i · ∇ x i p + 1 τ ∇ v i · ( v j p )+ ∇ x i Uτ γ · ∇ v i p + v τ ∇ v i ∇ v j [ D ij p ] . (14)Neglecting the term ∝ ∇ x i U because the forces almostcancel out in the solid phase (and their modulus is smallerthan the velocity modulus, as shown in Fig. 2), we caneasily check that Eq. (14) admits a solution of the form: p ( v x , v y ) ∝ exp h − α | v | − v ) i . (15)The shape of Eq. (15) corresponds to the Cartesian ver-sion of the velocity distribution shape numerically ob-served, i.e. Eq. (8). V. THE KINETIC TEMPERATURE
In Fig. 3, we show the first two moments of the ve-locity modulus distribution as a function of τ , namely h| v |i and h v i . The latter coincides by definition withthe kinetic temperature of a system of ABP and, forthis reason, we will denote h v i simply as “kinetic tem-perature” in the rest of the paper. For τ . − , thesystem is in the small persistence regime [case i)] andboth h| v |i and h v i are roughly constant with τ , beingapproximatively h| v |i ≈ v and h v i ≈ v . We recallthat, in this regime, the interparticle forces almost bal-ance and the velocity displays the same statistical prop-erties of the self-propulsion in such a way that the ki-netic temperature does not display any τ -dependence.Upon increasing τ , the values of h| v |i and h v i mono-tonically decrease reaching very small values. The morepersistent is the particle motion, the slower it becomesand, as a consequence, the kinetic temperature decreasesmonotonically with τ . After a crossover regime occur-ring for 10 − . τ . − , a clear power-law scalingwith τ appears for 10 − . τ . − in both the mo-ments. In particular, we have h| v |i ∼ v ( γτ ) − / and h v i ∼ v ( γτ ) − / , as clearly shown in Fig. 3. The valid-ity of these scalings ceases approximatively at τ = 10 − ,i.e. near the solid-hexatic transition. Starting from thisvalue of τ , both h| v |i and h v i decrease slower than ∼ τ − / and ∼ τ − / as the persistence time is increasedwithout showing any clear power-law scaling with τ .In what follows, we develop an exact, analytical predic-tion valid in active the solid-state for the kinetic tempera-ture which will explain the scaling with τ numerically ob-served shedding light also on the role of the other parame-ters. Indeed, the periodicity of the almost-solid structureand, in particular, its hexagonal order suggests switchingin the Fourier space to perform calculations [35, 57]. Asreported in Appendix A, the velocity correlation in theFourier space reads: h ˆ v q · ˆ v − q i = v τγ ω q (16)where ˆ v q is the Fourier transform of the velocity vector v and q = ( q x , q y ) is a vector of the reciprocal Bravaislattice. The factor ω q has the following form: ω q = − K h cos( q x ¯ r ) + 2 cos (cid:16) q x ¯ r (cid:17) cos (cid:16) √ q y ¯ r (cid:17) − i (17)where the dimensional constant K reads2 K = U ′′ (¯ r ) − U ′ (¯ r )¯ r . (18)and ¯ r is the average distance between neighoring parti-cles. To obtain the variance of the velocity, we need to togo back to the real space calculating the inverse Fouriertrasform in the origin: h v i = v N X q τγ ω q ) . (19)Taking the continuum limit in the q -sum and accountingfor the periodicity of the lattice, we restrict the integralto the first Brillouin zone: h v i ≈ v I (cid:20) τγ (cid:21) , (20)where I (cid:20) τγ (cid:21) = ¯ r |B| Z B d q τγ ω q ) , (21)and |B| is the area of the Brillouin region associated withthe hexagonal lattice. We remark that it is not possibleto approximate ω q for small q truncating at the quadraticorder since the integral diverges at q = 0. To the best of FIG. 3. h v i (red data) and h| v |i (blue data) as a functionof τ . The colored dashed lines are plotted as a eye-guides,while the solid red line represents the theoretical prediction,Eq. (20). Numerical simulations are realized with v = 50and ǫ = σ = γ = 1. our knowledge, Eq. (20) is the first analytical expressionfor the kinetic temperature of interacting self-propelledparticles that does not require fitting parameters. Weobserve that our expression increases quadratically with v in agreement with previous results [46] while the de-pendence on packing fraction and persistence time is con-tained in the integral I h τγ i . As shown in Fig. 3 (see thecomparison between red points and the solid red line),Eq. (20) is in fair agreement with numerical data whenthe system attains solid configurations for τ . − . Onone hand, I [0] ≈ τ , while, on theother hand, the numerical integration of I confirms boththe crossover regime and the scaling ∼ ( τ /γ ) − / in thelarge persistence regime. For τ & − , Eq. (20) underes-timates the values of h v i with respect to numerical databecause, for these values of τ , the structure of the systemis no longer a solid without defects. Thus, a fundamen-tal hypothesis behind the derivation of the prediction isviolated and, thus, Eq. (20) is no longer valid. In partic-ular, it has been already shown that, in the proximity ofdefects, active particles have kinetic energies much largerthan the ones in the absence of defects, as occurs in activesolid configurations [35]. This is a clue to understandingwhy the decrease of h v i with τ in active hexatic phasesis slower than the decrease for active solids. We remarkthat the integral I contains also the dependence on thepacking fraction through the constant K in Eq.(17). In-deed, K is mainly determined by the second derivativeof the potential calculated at ¯ r that is uniquely fixed bythe packing fraction in any homogeneous configurations.Thus, the growth of φ induces the increase of the ki-netic temperature through the non-linear derivatives ofthe function U (¯ x ). The explicit dependence on the po-tential shape is in agreement with previous studies basedon temperature definitions derived in simpler cases, i.e. FIG. 4. h v x i / h v x i (blue data) and h v x v y i / h v x i (red data)as a function of τ . The colored solid lines are eye-guideswhile the dashed black lines are marked in correspondence of h v x i / h v x i = 3 and h v x v y i / h v x i = 1, i.e. at the expected val-ues for a velocity following the Gaussian statistics. Numericalsimulations are realized with v = 50 and ǫ = σ = γ = 1. a one-dimensional particle confined through an externalpotential [51, 52]. A. Higher order moments and non-Gaussianity
In spite of the fact that the equilibrium-like Gaussianprediction is a good approximation of the velocity dis-tribution, at least in the large persistence regime, evenfrom Fig. 1 (d) clear deviations from the Gaussian the-ory are evident in its tails. To get a quantitative analysisof the non-Gaussianity, we report the behavior of thehigher-order moments in Fig.4. In particular, we showthe kurtosis of the v x -distribution, namely h v x i / h v x i .This observable is rather small (around ∼
1) in the small- τ regime as a result of the non-Gaussianity of the distri-bution (Fig. 1 (c)). On the contrary, in the Gaussian-like regime, for τ & − , the kurtosis shows just smalldepartures from the Gaussian prediction correspondingto h v x i / h v x i = 3. In particular, our numerical obser-vations reveal that h v x i / h v x i ≥
3, meaning that thetails of the distribution are a little fatter than the Gaus-sian prediction. Finally, for larger values of τ , i.e. for τ > − at the solid-hexatic transition point, the kur-tosis abruptly increases and the system departs from theGaussian-like regime. As already mentioned, this is con-sistent with the occurrence of intermittency phenomenain the hexatic phase [35] that manifest themselves alsoin high and non-Gaussian peaks in the time-trajectoryof the single-particle kinetic energy. As a further con-firmation, a similar scenario, consistent with the obser-vation regarding the kurtosis, occurs for the observable, h v x v y i / h v x i . In particular, h v x v y i / h v x i . − ≤ τ ≤ − (as expected in Gaussian regimes) and & VI. DISCUSSION AND CONCLUSIONS
In this paper, we have studied the properties of the ve-locity of highly packed systems of self-propelled particles(active hexatic and solid phases) to understand the influ-ence of the activity. A transition from Mexican-hat-likevelocity distribution (i.e. peaked in the proximity of a cir-cular crown with a radius larger than zero) to Gaussian-like velocity distribution is observed going from the smallpersistence to the large persistence time regime. Analyz-ing the velocity dynamics, we derive suitable approxima-tions to predict the functional form of the probability dis-tribution function of the velocities in both these regimes.Concerning the active solid, we derive by a Fourier-spacemethod a theoretical expression for the variance of thevelocity distribution, giving the kinetic temperature ofABP in this phase. Thus, on one hand, we have derivednew approximate analytical results concerning the veloc-ity distribution of ABP particles holding both near andfar from equilibrium and, on the other hand, we haveprovided the analytical expression for the active kinetictemperature in the solid phase.At least in homogeneous solid configurations, our an-alytical expression for the kinetic temperature increasesquadratically with the speed of the self-propelled parti-cles (that is proportional to the Peclet number). Thisquadratic scaling has been also observed by means of dif-ferent definitions of temperature, such as the active effec-tive temperature for homogeneous configurations [47, 49].In those cases, our kinetic temperature displays a mono-tonic decrease as a function of the packing fraction asalso for the effective temperature [46]. Our results showthat the kinetic temperature contains a strong depen-dence on the shape of the interacting potential in agree-ment with another temperature definition obtained in thecase of single-particle confined through external poten-tials [51, 52]. Thus, the concavity of the potential playsa fundamental role not only for confined non-interactingactive particles [71] but also for interacting systems, be-ing relevant to determine the velocity variance (and, thus,the kinetic temperature) in both cases. While the ki-netic temperature does not show a dependence on τ inthe small persistence regime, a power-law decay with τ isnumerically observed and theoretically predicted in thelarge persistence regime. We remark that our predictionsare valid in the active solid-state while does not workwhere orientational and/or positional orders are broken(i.e. active hexatic and liquid state, respectively). Inthose cases, the occurrence of large non-Gaussianity andintermittency phenomena in the time-trajectory of thekinetic energy are consistent with the failure of the the-oretical predictions [35]. Acknowledgements
LC thanks I. Petrelli, M. Cencini and A. Puglisi forfruitful discussions. LC and UMBM acknowledge sup-port from the MIUR PRIN 2017 project 201798CZLJ.
Appendix A: The kinetic temperature of activeparticles in the solid-state
To develop a prediction for the kinetic temperatureof self-propelled particles in the active solid state, weshall employ two approximations to simplify the dynam-ics, Eq. (1). i) The dynamics of each component of f ai is replaced by independent Ornstein-Uhlenbeck processeswith equivalent persistence time, τ = 1 /D r , and variance v in such a way that h| f a |i = v consistently with theABP model. ii) Each particle oscillates around a node ofa hexagonal lattice so that the total interparticle poten-tial is approximated as the sum of quadratic terms. Withthese two assumptions, the original dynamics, Eq. (1),becomes: ˙ x i ( t ) = f ai ( t ) − n.n X j ∇ i U ( | x j − x i | ) γ (A1) τ ˙ f ai ( t ) = − f ai ( t ) + v √ τ ξ i ( t ) , (A2)where the sum involves the nearest neighbors of the lat-tice node i , and the symbol ∇ i is the gradient with re-spect to x i . Introducing the displacement u i of the parti-cle i with respect to its equilibrium position, x i , namely u i = x i − x i , (A3) we obtain ˙ u i ( t ) = f ai ( t ) + Kγ n.n X j ( u j − u i ) (A4) τ ˙ f ai ( t ) = − f ai ( t ) + v √ τ ξ i ( t ) , (A5)being K the strength of the potential in the harmonicapproximation, i.e. U ≈ K ( u j − u i ) , that explicitlyreads: 2 K = (cid:18) U ′′ (¯ r ) + U ′ (¯ r )¯ r (cid:19) , where ¯ r is the lattice constant. Because of the linearity ofthe system, it is useful to switch to normal coordinates,in the Fourier space representation:ˆ u q = 1 N X i u i e − i q · x i (A6)ˆ η q = 1 N X i η i e − i q · x i , (A7)where ˆ u q and ˆ η q are the Fourier transform of u and f a ,respectively. The dynamics in the Fourier Space reads: ddt ˆ u q ( t ) = − ω q γ ˆ u q ( t ) + ˆ η q (A8) τ ddt ˆ η q ( t ) = − ˆ η q + v √ τ ˆ ξ q , (A9)where ω q = − K h cos( q x ¯ r ) + 2 cos (cid:16) q x ¯ r (cid:17) cos (cid:16) √ q y ¯ r (cid:17) − i where q = ( q x , q y ) are vectors of the reciprocal Bravaislattice. Defining ˆ v q as the Fourier transform of the ve-locity v , that satisfies ˆ v q = ddt ˆ u q , we can easily calculatethe steady-state equal time correlations in the Fourierspace, that is: h ˆ v q · ˆ v − q i = 2 v τγ ω q . (A10)Eq. (A10) is the final expression for the spatial veloc-ity correlation in the Fourier space and corresponds toEq. (16). [1] C. Bechinger, R. Di Leonardo, H. L¨owen, C. Reichhardt,G. Volpe, and G. Volpe, Reviews of Modern Physics ,045006 (2016).[2] M. Marchetti, J. Joanny, S. Ramaswamy, T. Liverpool,J. Prost, M. Rao, and R. A. Simha, Reviews of ModernPhysics , 1143 (2013).[3] ´E. Fodor and M. C. Marchetti, Physica A: StatisticalMechanics and its Applications , 106 (2018). [4] G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nar-dini, F. Peruani, H. L¨owen, R. Golestanian, U. B. Kaupp,L. Alvarez, et al., Journal of Physics: Condensed Matter , 193001 (2020).[5] S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek,and E. Bertin, Nature communications , 1 (2020).[6] N. Sep´ulveda, L. Petitjean, O. Cochet, E. Grasland-Mongrain, P. Silberzan, and V. Hakim, PLoS ComputBiol , e1002944 (2013). [7] S. Garcia, E. Hannezo, J. Elgeti, J.-F. Joanny, P. Sil-berzan, and N. S. Gov, Proceedings of the NationalAcademy of Sciences , 15314 (2015).[8] M. Basan, J. Elgeti, E. Hannezo, W.-J. Rappel, andH. Levine, Proceedings of the National Academy of Sci-ences , 2452 (2013).[9] J. M. Nava-Sede˜no, A. Voß-B¨ohme, H. Hatzikirou,A. Deutsch, and F. Peruani, Philosophical Transactionsof the Royal Society B , 20190378 (2020).[10] D. Sarkar, G. Gompper, and J. Elgeti, arXiv preprintarXiv:2006.04519 (2020).[11] Y. Fily and M. C. Marchetti, Physical Review Letters , 235702 (2012).[12] I. Buttinoni, J. Bialk´e, F. K¨ummel, H. L¨owen,C. Bechinger, and T. Speck, Physical Review Letters , 238301 (2013).[13] M. E. Cates and J. Tailleur, Annu. Rev. Condens. MatterPhys. , 219 (2015).[14] G. Gonnella, D. Marenduzzo, A. Suma, and A. Tiriboc-chi, Comptes Rendus Physique , 316 (2015).[15] Z. Ma, M. Yang, and R. Ni, arXiv preprintarXiv:2004.02376 (2020).[16] L. Caprini, U. M. B. Marconi, and A. Puglisi, PhysicalReview Letters , 078001 (2020).[17] J. Stenhammar, D. Marenduzzo, R. J. Allen, and M. E.Cates, Soft matter , 1489 (2014).[18] J. T. Siebert, J. Letz, T. Speck, and P. Virnau, SoftMatter , 1020 (2017).[19] F. Ginot, I. Theurkauff, F. Detcheverry, C. Ybert,and C. Cottin-Bizonne, Nature Communications , 696(2018).[20] A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri, andJ. Tailleur, Physical Review E , 020602 (2018).[21] M. N. van der Linden, L. C. Alexander, D. G. Aarts, andO. Dauchot, Physical review letters , 098001 (2019).[22] S. A. Mallory and A. Cacciuto, Journal of the AmericanChemical Society , 2500 (2019).[23] J. Grauer, H. L¨owen, A. Beer, and B. Liebchen, Scientificreports , 1 (2020).[24] F. Jose, S. K. Anand, and S. P. Singh, arXiv preprintarXiv:2004.01996 (2020).[25] P. Chiarantoni, F. Cagnetta, F. Corberi, G. Gonnella,and A. Suma, Journal of Physics A: Mathematical andTheoretical (2020).[26] J. Bialk´e, T. Speck, and H. L¨owen, Physical review let-ters , 168301 (2012).[27] A. M. Menzel and H. L¨owen, Physical Review Letters , 055702 (2013).[28] G. Briand, M. Schindler, and O. Dauchot, Physical re-view letters , 208001 (2018).[29] J. U. Klamser, S. C. Kapfer, and W. Krauth, Naturecommunications , 1 (2018).[30] P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo,G. Gonnella, and I. Pagonabarraga, Physical ReviewLetters , 098003 (2018).[31] S. De Karmakar and R. Ganesh, Physical Review E ,032121 (2020).[32] S. Paliwal and M. Dijkstra, Physical Review Research ,012013 (2020).[33] L. F. Cugliandolo, P. Digregorio, G. Gonnella, andA. Suma, Physical Review Letters , 268002 (2017).[34] P. Digregorio, D. Levis, L. F. Cugliandolo, G. Gonnella,and I. Pagonabarraga, arXiv preprint arXiv:1911.06366(2019). [35] L. Caprini, U. M. B. Marconi, C. Maggi, M. Paoluzzi,and A. Puglisi, Physical Review Research , 023321(2020).[36] R. Mandal, P. J. Bhuyan, P. Chaudhuri, C. Dasgupta,and M. Rao, Nature communications , 1 (2020).[37] J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet,Physical Review Letters , 088304 (2010).[38] L. Berthier and J. Kurchan, Nature Physics , 310(2013).[39] D. Levis and L. Berthier, EPL (Europhysics Letters) ,60006 (2015).[40] Z. Preisler and M. Dijkstra, Soft matter , 6043 (2016).[41] S. K. Nandi and N. S. Gov, Soft Matter , 7609 (2017).[42] G. Szamel, EPL (Europhysics Letters) , 50010 (2017).[43] S. K. Nandi and N. Gov, The European Physical JournalE , 117 (2018).[44] L. F. Cugliandolo, G. Gonnella, and I. Petrelli, Fluctu-ation and Noise Letters , 1940008 (2019).[45] S. Dal Cengio, D. Levis, and I. Pagonabarraga, PhysicalReview Letters , 238003 (2019).[46] I. Petrelli, L. F. Cugliandolo, G. Gonnella, and A. Suma,arXiv preprint arXiv:2005.02303 (2020).[47] A. Suma, G. Gonnella, G. Laghezza, A. Lamura,A. Mossa, and L. F. Cugliandolo, Physical Review E , 052130 (2014).[48] I. Petrelli, P. Digregorio, L. F. Cugliandolo, G. Gonnella,and A. Suma, The European Physical Journal E , 128(2018).[49] D. Loi, S. Mossa, and L. F. Cugliandolo, Soft Matter ,10193 (2011).[50] S. Mandal, B. Liebchen, and H. L¨owen, Physical ReviewLetters , 228001 (2019).[51] U. M. B. Marconi, A. Puglisi, and C. Maggi, Scientificreports , 46496 (2017).[52] L. Caprini, U. M. B. Marconi, and A. Puglisi, ScientificReports , 1386 (2019).[53] G. Szamel, Physical Review E , 012111 (2014).[54] A. Puglisi, A. Sarracino, and A. Vulpiani, Physics Re-ports , 1 (2017).[55] T. F. Farage, P. Krinninger, and J. M. Brader, PhysicalReview E , 042310 (2015).[56] G. S. Redner, M. F. Hagan, and A. Baskaran, PhysicalReview Letters , 055701 (2013).[57] L. Caprini and U. M. B. Marconi, arXiv preprintarXiv:2006.09551 (2020).[58] U. M. B. Marconi, N. Gnan, M. Paoluzzi, C. Maggi, andR. Di Leonardo, Scientific Reports , 23297 (2016).[59] ´E. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco,and F. van Wijland, Physical Review Letters , 038103(2016).[60] R. Wittmann, F. Smallenburg, and J. M. Brader, TheJournal of chemical physics , 174908 (2019).[61] L. Berthier, E. Flenner, and G. Szamel, The Journal ofChemical Physics , 200901 (2019).[62] L. Caprini and U. M. B. Marconi, Soft matter , 2627(2019).[63] E. Woillez, Y. Kafri, and N. S. Gov, Physical ReviewLetters , 118002 (2020).[64] C. Maggi, M. Paoluzzi, L. Angelani, and R. Di Leonardo,Scientific reports , 1 (2017).[65] Y. Fily, The Journal of Chemical Physics , 174906(2019). [66] L. Dabelow, S. Bo, and R. Eichhorn, Physical Review X , 021009 (2019).[67] L. Caprini, E. Hern´andez-Garc´ıa, C. L´opez, and U. M. B.Marconi, Scientific reports , 1 (2019).[68] S. Das, G. Gompper, and R. G. Winkler, New Journalof Physics , 015001 (2018). [69] D. Martin, J. O’Byrne, M. E. Cates, ´E. Fodor, C. Nar-dini, J. Tailleur, and F. van Wijland, arXiv preprintarXiv:2008.12972 (2020).[70] L. Caprini, U. Marini Bettolo Marconi, A. Puglisi, andA. Vulpiani, The Journal of Chemical Physics ,024902 (2019).[71] L. Caprini, U. Marini Bettolo Marconi, A. Puglisi, andA. Vulpiani, The Journal of chemical physics150