Time-dependent properties of interacting active matter: dynamical behavior of one-dimensional systems of self-propelled particles

We study an interacting high-density one-dimensional system of self-propelled particles described by the Active Ornstein-Uhlenbeck particle (AOUP) model where, even in the absence of alignment interactions, velocity and energy domains spontaneously form in analogy with those already observed in two dimensions. Such domains are regions where the individual velocities are spatially correlated as a result of the interplay between self-propulsion and interactions. Their typical size is controlled by a characteristic correlation length. In the present work, we focus on a novel and lesser-known aspect of the model, namely its dynamical behavior. To this purpose, we investigate theoretically and numerically the time-dependent velocity autocorrelation and spatio-temporal velocity correlation functions. The study of these correlations provides a measure of the average life-time and, thus, the stability in time of the velocity domains.


I. INTRODUCTION
In the last period, there has been an upsurge of interest in the study of systems of self-propelled particles that convert energy into directed and persistent motion. Certain biological systems such as run-and-tumble bacteria or crawling cells, as well as non-biological systems such as self-driven colloids or artificial swimmers, commonly referred as active matter, can be described in terms of effective models able to capture their salient features [1][2][3]. Active particles display a very rich phenomenology, such as their accumulation at the boundaries [4][5][6][7] and near rigid obstacles [8][9][10][11][12][13] or a kind of non-equilibrium phasecoexistence, known as Motility Induced Phase Separation (MIPS) [14][15][16][17][18] occurring even in the absence of attractive [19][20][21][22][23][24][25][26][27][28] or depletion interactions [29]. Self-propelled particles are far-from-equilibrium systems, showing several dynamical anomalies which have not a Brownian counterpart [30,31]. A recently reported phenomenon is the formation of large domains characterised by the tendency of the particles towards a common alignement of their velocities. This fact is somehow surprising since the particles have spherical symmetry, interact via central potentials, and are not endowed with an alignment mechanism. These domains have been observed numerically in Ref. [32] in two-dimensional systems of repulsive self-propelled disks (Active Brownian Particles) both at moderate packing fraction in the phase-coexistence region and at large packing fraction in homogeneous active liquid, hexatic, and solid phases [33], where domains with aligned velocities can still be observed [34]. A nonequilibrium phase activity-density diagram has been introduced to represent both homogeneous and inhomogeneous regimes [34] and the structural properties of the system have been compared with the typical size of the aligned domains. The model reproduces several experimental results regarding confluent cell monolayers [35][36][37][38], whose velocity fields display alignment patterns quite similar to the corresponding predictions of the ABP model. Hence, even such a simple model can account for the phenomenology of active matter systems at high density observed in experiments. In particular, the minimal ingredients to induce the velocity-alignment are i) the excluded volume interaction and ii) the persistent self-propulsion, while the long-range attractive force and/or the Vicsek-like velocity interaction are not strictly necessary.
At present, notwithstanding the existing detailed information about the velocity domains, obtained by measuring the equal-time spatial velocity correlations, very little is known about their dynamical properties. The aim of this work, which will be carried out by numerical and analytical methods, is to characterise how these domains evolve, how stable they are, and what controls their lifetime. To achieve this goal, we focus on high-density systems of self-propelled particles in one-dimension. The motivation of the choice of this low dimensional system is threefold: a) it considerably reduces the time cost of the numerical simulations, b) it is possible to develop an analytical theory and c) there are situations where biological swimmers move in highly confining geometries behaving as almost one-dimensional systems. This is the case of highly confined bacteria [39], such as E.Coli in microdevices with the size of a single particle, or molecular motors [40]. Recently, experimental studies with trains of swimming water-droplets have been performed in microfluidic square channels with sections equal to the droplets' diameter [41,42]. In this case, the local velocity alignment between neighbouring droplets has been experimentally observed.
The major outcome of the present investigation is encapsulated in Fig. 1, where the life-time of the velocity domains, t * , defined as the typical decaying-time of the spatio-temporal velocity correlations, is reported as a function of several values of the persistence time, τ (see Sec. IV for further details). Data coloring reflects the value of the correlation length, , of the spatial velocity correlations, which quantifies the average size of each velocity domain (see Sec. IV). While ∝ √ τ , our study reveals that t * displays a linear increase with τ . The larger the persistence time of the self-propelled motion, the larger are the typical size and life-time of the velocity domains.
The article is structured as follows: after the introduction of the model reported in Sec. II, we ascertain the existence of velocity domains in one-dimensional dense systems of active particles. The steady-state properties, e.g. spatial velocity correlations and correlation lengths, are studied in Sec. III. The major novel insight of this work is reported in Sec. IV, where the velocity autocorrelation and the spatio-temporal velocity correlations are numerically and theoretically investigated. A final discussion is reported in the conclusive section. Appendixes contain not only lengthy calculations giving support to the theoretical results of the paper, but also deeper insights into the problem.

II. MODEL
We study dense systems of N interacting self-propelled particles at density ρ 0 , employing the Active Ornstein-Uhlenbeck (AOUP) model [43][44][45][46][47][48][49]. The AOUP is a versatile and popular model of active matter that can reproduce many aspects of the phenomenology of selfpropelled particles, including the accumulation near rigid boundaries [50][51][52][53][54] or obstacles and the motility induced phase separation (MIPS) [55]. To the best of our knowledge, AOUP is also one of the simplest ways to model selfpropelled particles moving in one-dimension in the presence of mutual interactions and/or external forces [56][57][58]. In this paper, particles are constrained to move on a line of length L and are subject to periodic boundary conditions. The particles' positions, x i , evolve with the fol-lowing stochastic equation: where γ is the drag coefficient and we have neglected the thermal noise due to the solvent since, for many active colloids and bacteria, the thermal diffusivity is usually rather smaller than the effective diffusivity due to the active force [1]. The term F i represents the steric interaction between particles and is given by where U is a truncated and shifted Lennard-Jones (LJ) potential, namely where θ is the Heaviside function and both the energy scale and length scales σ are set to one, for numerical convenience. Since the potential chosen has a limited range only the first neighbours mutually interact. The term f a i models the self-propulsion of the particle i and is described by an Ornstein-Uhlenbeck process according to the following AOUP dynamics: where ξ i is a white noise with zero average and unit variance. The parameter τ is the persistence time and v 0 the active velocity associated with the self-propulsion force. We remark that, in the model employed in this paper, the self-propulsion acts independently on each particle, at variance with Vicsek-like models [59][60][61][62] or more complex dynamics where explicit couplings between velocities and self-propulsions are postulated [63,64]. A convenient method to study the AOUP is achieved by switching from (x i , f a i ) variables [55,65,66] to position, x i , and velocity, v i =ẋ i variables. In one dimension, the equation of motion (1) can be recast as: where the matrix, Γ ij , in general depends on the spatial coordinates of different particles and contains the derivatives of the potential: The original active over-damped dynamics is effectively mapped onto a passive under-damped dynamics with a space dependent friction and additional forces which mutually couple the particles' velocities. When τ γ 1 the inertial term can be neglected while for τ /γ ∂ xi ∂ xj U 1 uniformly in x i and x j , the term Γ ij reduces to δ ij . In these cases, Eq. (4) corresponds to an equilibrium overdamped motion with long-time effective diffusion coefficient D = v 2 0 τ .

A. Numerical simulations
The numerical study presented in this paper has been conducted by setting v 0 = 50 and varying the persistence time of the active force, τ . Simulations are realised with lengths, L, much larger than the persistence length of the active force, in such a way that the condition L τ v 0 is satisfied for the whole range of τ considered. Such a regime guarantees that boundary conditions do not play a significant role and that the one-dimensional version of traveling crystals does not occur [67,68]. At variance with previous studies [69][70][71][72][73], we consider high-density regimes, in such a way that the one-dimensional system of active particles is compact enough and neither defects in the periodic arrangement of the particles nor clusters can easily form: the system attains homogeneous configurations for the whole set of parameters numerically explored in this work.

B. The one-dimensional harmonic active crystal
To develop a suitable analytical theory and interpret the numerical findings, we have considered an approximate treatment of (4), by replacing the full LJ potential by its Taylor expansion truncated at the second order: where the lengthx = L/N is the average inter-particle separation. Such an approximation works quite well thanks to two conditions: the large packing regime and the limited range of the LJ interaction. In the simulations, the formation of defects is practically absent and this makes the mapping onto the active harmonic crystal a successful strategy. Being the Langevin equation relative to the active crystal (4) linear and diagonalisable via Fourier modes we can determine all the stationary one-time and two-time correlation functions within this approximation.

III. FORMATION OF VELOCITY DOMAINS: SPATIAL VELOCITY CORRELATIONS
Since the 1D AOUP shares the same physics as its 2D companion, also in this case, at high-density we expect the formation of spatial domains where the velocities statistically point towards a common direction and have the same modulus. This is pictorially shown in Fig. 2 (a) where several one-dimensional instantaneous configurations of the system are reported for different values of τ . Particles are represented by vertical segments and are colored according to their velocities, in such a way that a spatial region with the same color can be identified with a velocity domain. When τ increases, the size of the domains grows, as revealed by the color gradients in Fig. 2 (a).
To understand the formation of the velocity domains, we note that the dynamics (4) can be approximated with The additional velocity term, w i , is the mean velocity of the adjacent particles: Further details about the derivation of Eq. (7) are reported in the Appendix A. The first term is an effective Stokes force while the second is assimilable to thermal noise, therefore it does not lead to any kind of velocity alignment. The term, F i , represents the collective field which constrains the particles to form and maintain a lattice structure of periodicityx. The last term forces the velocity of particle i to assume values equal to the average velocity of its neighbours: it resembles an interactionà la Vicsek. We observe that this effective alignment force, dominant if τ U (x)/γ 1, is a genuine non-equilibrium effect and leads to the spontaneous local velocity alignment. On the other hand, in the opposite regime, i.e. τ U (x)/γ 1, the alignment force is negligible and neighbouring velocities are uncoupled as in the passive Brownian case. The increase of the oscillation amplitude of a singleparticle is captured by measuring the variance of the single particle position which is theoretically calculated in the Appendix B in the case of highly dense configurations and almost perfect lattices. Under these assumptions, we obtain: where the velocity variance, v 2 , is given by In the right hand side of Eq. (8), the first term coincides with the variance of a passive system (τ = 0) whereas the second term gives a negative correction resulting from the formation of correlated velocity domains.
Remarkably, for small values of τ , such that τ U (x)/γ 1, the velocity variance is nearly constant and is approximately v 2 0 , while for larger values of τ , v 2 decreases as τ 1/2 . Also an increase of ρ 0 ∝ 1/x makes U (x) larger and leads to a decrease of v 2 . Both observations qualitatively agree with the results relative to 2D interacting particles in non-harmonic potentials where a similar trend was reported [66]. Moreover, the variance of the position linearly grows with τ , as shown in Fig.3, in agreement with the qualitative observations of Fig. 2 . Indeed, as explained in Appendix B, Eq. (8) holds only for systems with Nx/v 0 τ 1 and, thus, the growth of the positional variance remains monotonic with τ . Interestingly, for τ = 10 −1 and mostly τ = 1, i.e. when x, the positional fluctuations are so large as to force particles to synchronise their fluctuations, i.e. to correlate their velocities.

A. Spatial velocity and energy correlations
To characterise the size of the velocity domains, we study the spatial velocity correlation functions, v(x)v(0) , in the steady-state adapting the strategy of Ref. [34] to a one-dimensional system. In this simple one-dimensional case, v(x)v(0) can be analytically predicted in the active harmonic crystal approximation, as shown in Appendix C. We find that, in the stationary regime, the fluctuation amplitude of each velocity mode has the following shape: where the frequency ω q reads: and q takes values from −π to π in the limit N 1. Inverting analytically the average (10) to determine its real space representation is not so easy without additional approximations, although some formulae can be found in terms of trascendental functions (see the Appendix C). The function represented in Eq. (10), being peaked around q = 0, is approximated employing a small q-expansion. Fourier transforming back to real space, the resulting velocity correlation displays an exponential spatial decay. In the continuum limit we obtain: where is the correlation length given by: Fig. 4 (b) reports a set of spatial correlations functions corresponding to different values of τ . It shows the excellent agreement between the numerical data and the theoretical predictions (12). Based on such an exponential behaviour, we argue that represents a measure of  12) and (14). Simulations are realised with v0 = 50 and the interaction is given by Eq. (2).
the size of the velocity domains, since particles within a distance ≈ are correlated and, roughly speaking, share similar velocities. The average size of a typical velocity domain scales as ∼ τ 1/2 and increases with ρ 0 , due to the dependence of on the curvature of the potential (higher ρ 0 ∝ 1/x means smallerx and, thus, larger U (x)), as shown by (13).
The spatial correlation of the kinetic energy, E(x) = v 2 (x), shows a similar trend. The two-point average E(x)E(0) , in the steady state, is used to measure the size of domains sharing the same kinetic energy. This observable approaches a non-vanishing asymptotic value, given by E 2 2 . For this reason, we numerically evaluate the normalised cumulant E(x)E(0) c = [ E(x)E(0) − E 2 2 ]/ E 2 2 . This spatial correlation is analytically calculated in the Appendix C and shows the following exponential shape: where E is its correlation length which reads: Theoretical predictions fairly agree with data as revealed both in panel (a) and (c) of Fig. 4, reporting E and the energy correlations, respectively. Also in this case, E corresponds to the average size of energy domains, which is smaller than by a factor 2 −1/2 and maintains the same scaling with the parameters of the model. As a consequence, in the one-dimensional case, the spatial energy correlations do not contain further information concerning the spatial velocity correlations.

IV. DYNAMICAL PROPERTIES OF THE VELOCITY DOMAINS
In the previous section, we have investigated numerically and theoretically the spontaneous formation of velocity domains by studying the velocity correlation functions. We have seen that, despite the absence of any alignment interactions, the velocities of different particles develop a correlation which increases with the the persistence of the active force. In this section, we investigate the time-dependent properties of the system, considering, in particular, the velocity autocorrelation and the two-time spatial velocity correlation function to unveil the dynamics of the velocity domains and estimate their life-times and permanence.

A. Velocity autocorrelation function
In Fig. 5, we report the steady-state normalised velocity autocorrelation functions (VACF), v(t)v(0) / v 2 , as a function of t/τ for several values of τ , and explore both the large and the small persistence regimes where velocities are spatially uncorrelated. As shown in Fig. 5 (a), this observable decays within a typical time of order τ for the whole range of persistence times numerically explored. In the small τ regime, i.e. when τ U (x)/γ 1, the VACF displays the same exponential shape as the active-force autocorrelation, i.e. f a i (t)f a i (t ) = Dγ 2 τ e −|t−t |/τ , (see the yellow curves in Fig. 5). Such a behaviour is a direct consequence of Eq. (4) when τ U (x)/γ 1 since, in this limit, Γ ij ≈ δ ij . When τ increases, the shape of the VACF changes and the relaxation process becomes faster as measured with respect to rescaled time t/τ . Then, for further values of τ , in particular, for τ 10 −4 , the rescaled VACFs collapse onto the same curve as shown in Fig. 5. Quite surprisingly, the VACFs v(t)v(0) / v 2 assume negative values for t > τ , except in the small-τ regime, for τ U (x)/γ 1, where the decay is exponential as already discussed. In the former case, the autocorrelations decay very slowly (with a power-law behaviour) towards zero from negative values, as zoomed in Fig. 5 (b). The sign inversion of v(t)v(0) / v 2 is understood in terms of the presence of two competing mechanisms, the active force producing a negative contribution to the correlations and the restoring force of the individual oscillation modes driving the VCF towards its final value. Each mechanism is characterised by a different time scale, τ and γ/ω 2 q , respectively. Only if the active force acts on a time-scale not too small compared with γ/ω 2 q the above phenomenon can be observed, therefore it has not a passive Brownian counterpart and represents a pure nonequilibrium collective effect. Fig. 6 (a)-(c) displays the spatio-temporal velocity correlation function (VCF), v(x, t)v(0, 0) / v 2 , as a function of t/τ for several distances, x/σ, roughly from x = 0 to x ∼ . In particular, panel (a), (b) and (c) are obtained at different values of τ , corresponding to configurations ranging from those with negligible spatial correlations to those highly correlated. For the smallest value of τ , reported in panel (a), even nearest-neighbour pairs sitting at distance x/σ = 1 are practically uncorrelated in time, while for larger τ -values very far particles display a pronounced time-correlation (panels (b) and (c)). In the cases (b) and (c), the spatio-temporal correlations for "small" separation x/ closely resemble the corresponding VACF. On the contrary, for x ∼ , the VCF v(x, t)v(0, 0) / v 2 displays a sort of plateau for an initial time-window (always τ ) until it collapses onto the VACF curve for t τ .

B. Spatio-temporal velocity correlation functions
At the origin of such a plateau is the fact that the state of the particle placed at x = 0 changes under the influence of the active force only after a time of order τ . However, the information about such a change does not reach a second particle belonging to the same velocity domain (hence correlated with the first) and located n-sites away, before a time which increases with their separation has elapsed. Finally, for x , any sort of spatio-temporal correlation is absent because the two particles do not belong to the same domain.
Based on the observation of the collapse of the VCF data onto the VACF curve for x ∼ , we identify the average life-time, t * , of a domain with typical size as the time at which the normalised autocorrelation approaches the value 1/e being e the Neper number. This typical time t * shows a linear increase with τ above the dashed black line in Fig. 1 which marks the first value with nonvanishing spatial correlations. Below this line, even nearest neighbour particles are almost independent and velocity domains do not form. In this last case, since the domain contains only a single particle, t * is nothing but the VACF relaxation time. Such a quantity coincides with τ in the equilibrium-like regime where τ U (x)/γ 1 and shows a non-linear increase as a function of τ in the crossover regime before the linear for τ U (x)/γ 1.

C. Theoretical predictions
To predict the shape of the spatio-temporal VCF, we consider the correlation of th Fourier modes as in the case of the steady-state spatial velocity correlations. Solving the dynamics in the active harmonic crystal approximation, we obtain: where the factor ω 2 q is defined by Eq. (11). At variance with the spatial velocity correlations (Eq. (10)), here, the approximation for small q is no longer valid to predict the whole time-behaviour (see Appendix D for more details about the derivation and the approximations involved). The reason is that the wave-vector dependent VCF represented in Eq. (16) is made of two contributions each varying with its own relaxation time and both containing a divergence. Only by handling them together the two divergences cancel out. Both for VACF and spatio-temporal VCF, Figs. 5 and 6 show the excellent agreement between data obtained via numerical simulations and the numerical integration of Eq. (16) (normalised with v 2 ).

Short-time approximation
The integral in Eq. (16) can be evaluated numerically in a straightforward way but must be handled with care to extract analytical predictions about the temporal decay of the spatial autocorrelation function. In the Appendix D, we derive a suitable approximation, holding for t τ , which consists in expanding the exponentials for small (t/τ ) and resumming a class of terms. In this way, the small-time decay of the autocorrelation is predicted and reads: where the integral is given by .
(18) This class of integrals can be performed exactly in terms of modified Bessel functions of the first kind and integer order [74]. Prediction (17) fairly agrees with data in a time-window smaller than τ while the differences between small-time theory and simulations occur approximatively for t > 0.3 τ , as shown in Fig. 7 for two different values of τ . In particular, Eq. (17) cannot reproduce the negative values assumed by the autocorrelation for t ∼ τ . In addition, in the short-time regime t/τ < 1, and, in particular, from Eq. (17), it is possible to derive the leading correction of the departure of v(x, t)v(0, 0) from its equal-time value: where τ K = γ/2U (x) and n ≈ x/x. This prediction explains the plateau observed for large distances when t/τ is small. Indeed, the prefactor of the first-order timecorrection in Eq. (19) becomes smaller as x increases through the n-dependence. In agreement with the numerical observations, a longer time is needed before the VCF with large x deviates from v(x, 0)v(0, 0) .

Long-time behaviour
As already mentioned, in the small persistence time regime the VACF has a pure exponential decay and, thus, for t τ , approaches zero from positive values. As numerically revealed, the decay is more complex in the large persistence regime, i.e. τ U (x)/γ 1. While the VACF and the VCF are positive in the early stage of the relaxation process as shown by direct inspection of Eq. (16), they become negative at larger times. Each mode gives a positive contribution to the q-integral up to a wave-dependent crossover time, t * q = −τ ln(τ ω 2 q /γ)/(1 − τ ω 2 q /γ) > 0, defined as the instant when the function in Eq. (16) shows its first zero. Such a q-dependent time decreases as q grows and, thus, the final decay towards zero of the VACF (for t τ ) is controlled by the relaxation of the long wave-length modes since the short wave-lengths give negligible contributions to the integral. Roughly for t > τ the e −t/τ term in formula (16) plays a negligible role. Performing the qintegration of Eq. (16) we obtain the following long-time approximation: which shows that, as t → ∞, the VACF vanishes from negative values with an inverse power-law behaviour, explaining the behavior observed in Fig. 5. The prediction (20) is derived in the Appendix D and numerically checked in Fig. 7 for two different values of τ showing a good agreement with data. In addition, we remark that, in the small persistence regime, the relevant amplitudes in Eq. (16) remain positive, as expected for passive systems. Finally, it is interesting to realise that the nonmonotonic behaviour of the two-time correlation function discussed above is different from the monotonic behaviour of the linear response function, i.e. the response of the system to a a small impulsive perturbation on the positions of the self-propelled particles. The calculation, reported in Appendix E, clearly shows that response and correlation are not proportional.

V. CONCLUSIONS
The present work contains results concerning both the steady-state and time-dependent spatial properties of a system of interacting active particles. To the best of our knowledge, such an analysis has not been performed so far, in particular, for what concerns the characterisation of the two-time correlations in the non-equilibrium steady state of self-propelled particles.
The interest in such a study is twofold: the ensemble averages of the single time observables are constant in time but already contain interesting information about the spatial organisation of the system, namely display the presence of domains where the velocities of the particles are strongly correlated. The two-time observables reveal how the system responds to external stimuli or how a given spatial structure lasts in time. Our study employs a one-dimensional model of interacting self-propelled particles evolving under the Active Ornstein-Uhlenbeck dy-namics, as a test system. The numerical study is conducted by choosing a truncated and shifted Lennard-Jones (LJ) potential and large packing conditions. As a preliminary step, we have studied the stationary properties of the system to characterise the size of the velocity domains and focused on the spatial velocity correlations. Then, we investigated the dynamics of the velocity domains studying the autocorrelation and the spatio-temporal velocity correlation functions and determined the average life-time of a typical velocity domain. Our numerical study has been supported by theoretical arguments to derive the temporal dependence of the correlations. Thanks to the lack of an appreciable density of defects or voids we can faithfully describe the system as a one-dimensional lattice of coupled harmonic oscillators, whose elastic and spacing constants are derived from expanding quadratically the LJ potential near its equilibrium values. The resulting active harmonic lattice can be solved analytically in terms of Fourier normal modes giving explicit formulae for one-time and two-time observables.
The theoretical method employed in this paper bears some similarities with the one used in active polymer theory [75][76][77][78][79], i.e. the shape of the excitation spectrum to make an example. However, the observables, here considered, are quite different from those usually studied in the context of active polymers, described by one-dimensional chains in two or three-dimensional spaces. Those studies are mainly concerned with diffusive [80], or configurational properties [81][82][83], while the present work regards the dynamical properties of a one-dimensional single file system where the interplay between low dimensionality, steric interactions, and self-propulsion determines an offequilibrium, highly correlated motion not studied so far.

Appendix A: The harmonic active crystal
In this Appendix, we illustrate the approximations employed to develop the theoretical formulas that in the main text have been compared to the numerical data. We replace the full LJ potential by its Taylor expansion around the average particles' positions,x. The expansion is truncated at the second order and gives: Replacing U tot with U hc , the matrix Γ ij assumes a simple form and the dynamics (4) becomes: Remarkably, the force term, maintains the average distance between neighboring particles and is derived from U hc . Nevertheless, its shape is not particularly significant. Splitting the square brackets in Eq. (A2), we obtain the dynamics (7). Representing Eq. (A2), in Fourier Space, we obtain the equations of motion of the harmonic crystal: whereû q andv q are the Fourier transforms of the displacement and velocity, respectively, which explicitly are whileξ q is a noise term obtained from the fourier transform of ξ i . The reciprocal lattice wavevectors are q = 2πk/N with k = 1, N , while the frequency ω q is such that: We remark that, in the limit N 1, the sums can be approximated by integrals In the following, we shall maintain a discrete notation that identifies each particle with an integer index. The comparison between the discrete theoretical and the xdependent numerical observables is made by noting that, in the high-density regime considered, the particle's coordinates approximately satisfy the relation x = nx, making possible replacing v n by v(x).
Setting u 2 = (x n − nx) 2 (which is n-independent) and K = U (x), for the sake of simplicity, the resulting autocorrelation of the displacement, u, reads: .
The first integral in Eq. (B3) is divergent for N → ∞, i.e. when performed from 0 to π/2, at variance with the second one. However, for a finite but large sample with N 1, by discarding a small interval near the origin we obtain from Eq. (8) the following result Eq. (B4) makes sense just if u 2 > 0, meaning that does not hold for arbitrary values of τ . Indeed, the validity of Eq. (B4) is restricted to the regimes where the persistence length, v 0 τ , is smaller than L = Nx (i.e. the regimes numerically explored in this paper). A simple integration gives the formula for the velocity variance reported in Eq. (9). Interestingly, we can relate the velocity fluctuation to the positional fluctuation through the following relation: where the first term in the right hand side is the well known variance of the displacement fluctuations in the white noise limit, τ = 0 [84,85]. Thus, the single units, comprising the active harmonic crystal, have fluctuations smaller than those of the corresponding passive crystal. This observation agrees with the well-known results of single active harmonic oscillators, which display an enhanced rigidity with respect to their equilibrium counterparts.

Appendix C: Spatial correlations of the velocity and energies
The equal-time velocity pair correlation function in Eq. (12) is obtained by expanding ω 2 q up to quadratic order in q 2 and extending the Fourier integrals to the whole real axis. Such an approximation is valid when the discretization of the lattice becomes irrelevant. Let us write where g n (x) = 1 2π π −π dq cos(nq) x−cos(q) . For small values of the index n, g n can be expressed in terms of elementary functions, for instance g 1 (x) = x √ , whereas in general we have to write it in terms of Hypergeometric functions. By expanding the result in powers of , one verifies that the terms are very well reproduced by the formula (12). This procedure yields the same result as if we expand up to the second order in q the denominator in Eq. (C1) and perform an integration over the whole real axis.
By a simple extension of the previous method we may evaluate the energy correlation function. In our discrete notation we define E l = v 2 l and consider where we have used the Gaussian property of the distribution to factorize the operator average. Taking into account Eq. (12), we obtain the final expression: which coincides with Eq. (14).

Appendix D: Spatio-temporal correlation functions
For t > t , the Fourier Transform of the two-time displacement correlation function, in the steady-state, reads (D1) while the Fourier Transform of the velocity correlation function is given by Eq. (16). Thus, to come back to real space, we have to calculate the integral .
(D2) The apparent singularity of the denominator is eliminated by the concomitant vanishing of the numerator.

Short-time approximation in real space
Let us consider the following integral related to the velocity correlation by the identity v n (t)v 0 (0) = d 2 σ 2 n (t) dt 2 . After Taylor expanding the first and second exponential in powers of the rescaled time, t/τ , we rewrite it as Up to this order, we may approximate the above result by the following expression: dq cos(qn) (Γ(q)) 2 (1 − e −Γ(q)t/τ ) , (D4) with Γ(q) = 1 + τ ω 2 q γ . Finally, differentiating with respect to t, we get the short-time formula for the velocity correlation function: v n (t)v 0 (0) = v 2 0 2π π −π dq cos(qn) Γ(q) e −Γ(q)t/τ .
Perhaps, the best strategy to display the presence of a small t interval where the velocity two-time correlation function for particles separated by a distance n is to study the first time derivative of such a function: where τ K = γ 2K . The integration can now be performed exactly: where I n (x) is the first-kind modified Bessel function of order n. Using standard properties of the Bessel functions, we find the following short-time aproximation: v n (t)v 0 (0) ≈ v n (0)v 0 (0) − 1 n! v 2 0 τ 1 2τ K n t n+1 n + 1 .
(D8) That is to say, the relaxation rate of the pair correlation becomes smaller with increasing separation of the the two particle considered since it is proportional to (t/τ ) to a power equal to (n + 1). This formula coincides with Eq. (19) and explains the formation of plateau regions. These plateau become more and more evident when the particle's separation increases.

Long-time approximation
We define the crossover time, t * q , as the time at which the VACF vanishes. From Eq. (16), it is straightforward to see that t * q = −τ ln(τ ω 2 q /γ)/(1 − τ ω 2 q /γ). For times larger than t * q , it is safe to neglect the contribution of the e −t/τ term in Eq. (D2). The modes with small values of q yield the main contributions to the integrals Eq. (D2) and the result is negative because the constant term in the denominator is larger than the q-dependent term. To evaluate Eq. (D2) in the long-time regime t/τ 1, we first consider the following integral: Finally, the sought correlation is found by differentiating such a formula with respect to t: v n (t)v 0 (0) ≈ − v 2 0 τ 2 Kπ/γ t −3/2 (1 − 3 32 (4n 2 − 1) (K/γ)t + . . . ) .
This equation coincides with Eq. (20) at the first order and the other terms of the expansion are negligible because of the large value of K.

Appendix E: Response function
The impulse response function for the AOUP was obtained by Szamel [56] for the harmonic case and by Caprini et al. [86] in the general case. In the present Appendix, we derive the response of the system by adding a small impulsive forceĥ q (t) =ĥ 0 q δ(t) to the equation describing the evolution of the q-component of the displacement.
Following a standard procedure (see for instance, Ref [87]) to compute the response to an impulsive perturbation acting on a particular q-mode we arrive at the following formula R q (t) = − û q (t) ∂ ∂û q ln P q (t = 0) , where P q is the known steady state distribution of the q mode, explicitly given by The response function, defined as: turns out to be R q (t) = e −(ω 2 q /γ)t , using Eq. (E1) and the linearity of the system. Remarkably, the response function only depends on the frequency ω q but not on the persistence time, τ .