Inhomogeneous entropy production in active crystals with point imperfections

The presence of defects in solids formed by active particles breaks their discrete translational symmetry. As a consequence, many of their properties become space-dependent and different from those characterizing perfectly ordered structures. Motivated by recent numerical investigations concerning the nonuniform distribution of entropy production and its relation to the config-urational properties of active systems, we study theoretically and numerically the spatial profile of the entropy production rate (EPR) when an active solid contains an isotopic mass defect. The theoretical study of such an imperfect active crystal is conducted by employing a perturbative analysis that considers the perfectly ordered harmonic solid as a reference system. The perturbation theory predicts a nonuniform profile of the entropy production extending over large distances from the position of the impurity. The EPR decays exponentially to its bulk value with a typical healing length that coincides with the correlation length of the spatial velocity correlations characterizing the perfect active solids in the absence of impurities. The theory is validated against


I. INTRODUCTION
Active matter comprehends many systems of biological and technological interest such as bird flocks, cell colonies, spermatozoa, and Janus particles, to mention just a few of them [1,2].All these systems are capable of self-propulsion, namely a mechanism that converts energy from the environment into directed and persistent motion and drives them out of equilibrium.Based on experimental evidence, such a self-propulsion or active force is represented at a coarse-grained level by a stochastic process with memory.In other words, the value of the active force acting on a given particle at a given instant is correlated with the values it took in the past.
In this study, we shall focus on the properties of active matter at high density, a regime characterizing several systems ranging from biological tissues and cell monolayers [3][4][5] populating our skin, to dense colonies of bacteria [6][7][8] capable of self-organizing into active twodimensional crystals of rotating cells [9].Moreover, solidlike configurations have been observed in systems of active Janus colloids [10][11][12][13] and active granular particles [14][15][16][17][18]. Several numerical and theoretical studies investigate the effect of the active force on high-density phases of active matter, such as liquid, hexatic, and solid [19][20][21][22][23][24][25][26][27].In two dimensions, the activity shifts the liquid-hexatic and hexatic-solid transition to larger values of the density [28,29], somehow, increasing the effective temperature of the system and broadens the size of the hexatic region that in the passive case is quite narrow [28].More recently, it has been found that dense phases of active matter display spatial velocity correlations [30][31][32][33][34][35] a feature absent in equilibrium systems but observed in active glasses [36][37][38].This is a phenomenon of dynamical origin determined by the tendency of the particle velocities to align spontaneously even in the ab-sence of direct alignment force [39]: it results from the combined action of the persistence of the direction of motion and the steric repulsion among the particles, while attractive interactions can even induce a flocking transition [40].
To mark the difference between an active and an equilibrium solid with similar structural properties, one can apply the tools of stochastic thermodynamics [41,42].In particular, the so-called entropy production rate (EPR) provides a quantitative measure of the distance of a system from equilibrium [43][44][45][46].This analysis discriminates between non-equilibrium steady states which produce entropy [47][48][49][50][51] and truly equilibrium states whose EPR is zero.The non-vanishing of the EPR is a universal feature of non-equilibrium systems and occurs when their dynamics break the time-reversal symmetry, i.e. the detailed balance condition is violated.In the present problem, the system produces both entropy and steady probability currents, a situation that never occurs under equilibrium conditions.
Being intrinsically out of equilibrium, active matter is an ideal platform to investigate entropy production and shed light on several general properties of non-equilibrium systems.However, except for specific cases [52,53], such as non-interacting active particles [54,55] and harmonically confined systems [56,57], analytical results for the EPR are difficult to achieve: in general, the EPR for non-linear confining forces or interacting systems has been obtained numerically [58,59].Other numerical investigations focused on the study of the EPR of systems spatially inhomogeneous through particleresolved simulations [60][61][62] or using field-theoretical descriptions [63][64][65][66].
Recently, we have studied active crystals and found that their vibrational excitations are of two different kinds: the first is identified with the conventional col-lective oscillatory modes, known as phonons, and the second describes additional vibrational excitations, absent at equilibrium and termed entropons because are the modes associated with the entropy production of the system [67].Under small deviations from equilibrium conditions, entropons coexist without interfering with the conventional phonons, the equilibrium-like excitations.Entropons vanish in equilibrium whereas dominate over phonons when the system is far from equilibrium.While we have a satisfactory description of the EPR in an ideal solid phase much less is known when its order is altered by the presence of imperfections such as surfaces, defects, or other departures from the perfect periodic arrangement of the active particles.
The aim of this paper is to investigate the EPR in active systems with inhomogeneities and determine its spatial distribution in the non-equilibrium steady-state.We consider a case that lends itself to analytical and numerical scrutiny: the active crystal containing point imperfections that are due to particles with different masses and destroy the periodic order characterizing perfect crystalline structures.This is a classical problem of Solid State Physics [68] where it is studied to understand the localization of the phonon modes near the impurity.We employ this model to shed some light on the EPR in active systems with broken translational invariance by developing a suitable perturbation theory around the perfect crystal state by considering a small defect mass.Our method predicts analytically the spatial profile of the EPR as a function of the distance from the lattice imperfection and relates this feature to the existence of a velocity correlation profile.
The paper is structured as follows: in Sec.II, we present the model to describe a solid formed by active particles and calculate the entropy production employing a path-integral technique.Section III reports the main results of the paper: we derive the perturbative method introduced to calculate the spatial profile of the EPR in the presence of a mass impurity in the solid.We evaluate explicitly the zeroth-order perturbative EPR, i.e. the entropy production rate of a perfectly ordered crystal, and the first-order perturbative correction that describes the effect of the point imperfection that breaks the translational discrete symmetry of the periodic array.Details about the derivations are presented in the appendices to render the exposition.Finally, the conclusions are presented in Sec.IV.

II. MODEL
We investigate a solid formed by N ABP's [39,[69][70][71][72][73][74][75] in two dimensions, in a square box of size L × L and apply periodic boundary conditions.The evolution of the position, x p and velocity, v p = ẋp of each ABP of mass m p (with p = 1, ..., N ) is governed by the following an underdamped stochastic equation where ξ p is a white noise with zero average and unit variance.The coefficients γ and T are the friction coefficient and the temperature of the solvent bath, respectively.For equal masses, the ratio, mγ corresponds to the typical inertial time, τ I , representing the relaxation time of the velocity in equilibrium systems (we remark that in active systems the relaxation of the velocity is determined both by τ I and τ [32]).
The active force, f a p , provides a certain persistence of the particle trajectory and drives the system out of equilibrium.In the absence of any other force but the friction, f a p and for τ I → 0 the ABP's would travel at the swim velocity, v 0 , as shown by the relation: The stochastic vector n = (cos θ p , sin θ p ) is a unit vector, whose orientation is determined by the angle θ p , subject to Brownian motion θp = 2D r η p , where η p is a white noise with zero average and unit variance and D r is the rotational diffusion coefficient.D r determines the persistence time of the dynamics τ = 1/D r , i.e. the average time needed by a particle to change direction [76,77].The analysis of the mean square displacement in the independent particle limit leads to the introduction of the so-called active temperature, T a = v 2 0 γτ , that is an increasing function of both τ and v 2 0 .The force F p represents the inter-particle interaction due to a pairwise potential, U tot = i>j U (|r i − r j |), that we choose as a purely repulsive and given by the shift-and-cut WCA potential for r < 2 1/6 and zero otherwise.The parameters ǫ and d 0 represent the energy scale and the particle diameter, respectively.To consider a solid configuration in numerical simulations, the packing fraction of the system is set to φ = ρd 2 0 π/4 = 1.1 that for the range of parameters explored in this study will result in a solid configuration as in the phase diagram reported in Ref. [30].
Above a certain density, the particles spontaneously arrange themselves to form an almost regular triangular lattice.Assuming that this configuration corresponds to the minimum of the total potential energy of the system, we Taylor expand U tot around it up to second order in the displacements, u p .[30,78].These are defined as the deviations of the particles' coordinates from the perfect lattice positions, r 0 p , through the relation u p = r p − r 0 p .Within this approximation, the force F p acting on the p particle reads:

Mass defect
Active crystal where the sum is restricted to the first neighbor particles and ω E is the Einstein frequency of the solid: ω E depends explicitly on the derivatives of U calculated at x, i.e. the average distance between neighboring particles of the solid (i.e. the lattice constant), which is determined by the packing fraction.

A. Calculation of entropy production
The stochastic thermodynamics [41,42,79,80] is a powerful tool to measure how far from thermodynamic equilibrium is a system, i.e. its degree of irreversibility.Such information is contained in the so-called entropy production rate (EPR), ṡ, which can be determined by considering the probabilities of the trajectories connecting two different states of the system.The EPR is expressed in terms of path-probability by resorting to pathintegral techniques [50,[81][82][83] as where the symbol • represents the steady-state average performed over the realizations of the noise and P and P r are the path probabilities of the forward and backward trajectories of the system, respectively.These probabilities depend on the whole time history of the dynamical variables of the system (x p , v p , f a p ) conditioned to their initial values (x p (0), v p (0), f a p (0)).In the case of equilibrium systems, in virtue of the detailed balance condition, which is tantamount to the probabilistic time-reversal symmetry, this ratio is one and the EPR vanishes.To estimate P and P r , let us remark that the probability of the trajectory of a stochastic system is uniquely determined by the probability of observing a path-trajectory of the stochastic noises, that in our case have a Gaussian distribution.Therefore, one performs a transformation from the noise variables to the dynamical variables by using the equation of motion (1) together with Eq. ( 3).In doing so, we neglect the determinant of the transformation because, in the present case of additive noise, this term does not contribute to the EPR.Applying this procedure, the probability of forward and backward trajectories are expressed as P ∼ e A and P r ∼ e Ar , respectively, where A and A r are actions associated to backward and reverse dynamics.The action A is obtained by expressing the Gaussian distribution of the noise variables, ξ p , in terms of the state variables (x p , v p , f a p ) using the relation between the two sets of variables given by Eq. ( 1) with the result: Here, we have not included in the action the contribution associated with the rotational noise, η p , since it is known that the simple Brownian process of Eq. ( 3) does not generate entropy.The action of the backward trajectory A r can be obtained by applying the time-reversal transformation to the dynamics (1) and considering the parity of the dynamical variables, (x p , v p , f a p ), under this transformation.By denoting with the subscript r the time-reversed variables, we assume: for each particle (above, the particle index has been suppressed for notational convenience).In this way, the backward action, A r , reads where we have again neglected the irrelevant contribution of the angular dynamics, for the same reasons given above.
Performing algebraic calculations, the expression for ṡ can be analytically derived and reads: where ṡp is the entropy production generated by a single particle given by ṡp The expression b.t.means boundary terms, i.e. the additional terms that do not contribute to the steady-state entropy production and vanish in the long-time limit.Such a result holds, in general, for underdamped active particles subject to a persistent active force and in contact with a thermal bath, independently of their density and of the dimension of the system.By introducing the spatial Fourier transform of the dynamical variables denoted by hat-symbols, and by neglecting boundary terms, ṡp can be calculated in the Fourier space, in particular, in the frequency domain, obtaining where f a p (ω) and vp (ω) are the Fourier transforms in the frequency domain of the active force and velocity of the p particle, respectively, defined in appendix A and the symbol "c.c" stands for complex conjugate.The EPR, ṡ, is proportional to the frequency integral of the real part of the cross-correlation between the Fourier components of active force and velocity.We also introduce the spectral entropy σ p (ω) for each particle as the integrand in Eq. ( 13) Note that this term corresponds to the spectral dissipation of the particle p due to the active force, divided by the temperature of the bath.

III. ENTROPY PRODUCTION OF ACTIVE SOLIDS WITH IMPURITIES
A real solid may contain various kinds of imperfections or surfaces which affect the properties of the perfect crystal.For the sake of simplicity, we confine ourselves to isolated defects such as substitutional particles of different mass [68].
To understand the effect of substitutional impurities on the EPR, we modify the mass of one particle, setting m 1 = m and m p = m for all remaining particles.Since this operation breaks the discrete translational symmetry of the lattice, both the displacement and the EPR, ṡp , become position dependent.In the continuum limit, the local EPR, ṡp → ṡ(r), becomes a function of the distance r from the location of the impurity.A sketch of the model is shown in Fig. 1.Notwithstanding the problem described by Eqs.( 1) and ( 5) is linear and one could use a numerical matrix inversion method to determine with great accuracy the solution u p and ṡ(r), we are interested in getting explicit predictions.Therefore, we apply an analytical perturbative approach choosing the mass difference δm = (m 1 − m) as a small parameter.

A. General strategy of the perturbation scheme
To proceed analytically, we choose m 1 = m + δm with |δm| ≪ m and apply a perturbative method by expanding the solution in powers of the small parameter δm/m ≪ 1. Explicitly the modified equation of motion for the imperfect lattice reads where the Kronecker delta function δ p0 selects the particle p = 0 corresponding to the imperfection.The force F p is approximated within the harmonic approximation of Eq. (5).
By introducing the continuous Fourier transforms in the frequency domain ω of the displacement ûp = ûp (ω), the velocity vp = iω ûp (ω), the active force f a p (ω) and the white noise ξp (ω) Eq. ( 15) can be rewritten as where the time-Fourier transform of a generic observable is denoted by the hat symbol and the Einstein summation convention used.The matrix elements L pk (ω) have the following form where the sum over the index j runs over the nearest neighbors (k+j) of the particle p.When δm = 0, Eq. ( 16) corresponds to the dynamics of a perfect lattice and can be rewritten with the help of the lattice Green's function, G pn (ω), in the ω-representation as: where G pn (ω) is expressed with the help of the discrete spatial-Fourier transform: where the sum runs over the dimensionless reciprocal lattice wave vectors q (see appendix C).The quantity ω(q) represents the dispersion relation of the vibrational modes of the lattice and depends on the Einstein frequency, ω E and on the lattice structure.It is obtained by solving the secular equation associated with the lattice harmonic oscillations [84].Since the perturbative method discussed below is independent of the specific lattice structure, being the relative information contained in ω(q), we postpone (see Eq. ( 34)) the presentation of the explicit expression of ω(q) for the case of the triangular lattice in two dimensions and nearest neighbor interactions.By solving Eq. ( 16) when δm = 0, we obtain The entropy production ṡ can now be calculated by using Eq. ( 13), which requires the knowledge of the crossdynamical correlation between the active force and velocity.We multiply Eq. ( 20) by f a n (−ω), take the average over the realizations of the noise, and obtain Multiplying both sides of Eq. ( 21) by (iω) we find a relation for the average vp (ω) f a p (−ω) .To proceed further, we need the explicit expression of the correlation f a p (ω) f a p (−ω) , which is estimated by approximating the ABP dynamics with the AOUP model [85][86][87][88][89].This strategy has been often adopted with success in the literature to get analytical predictions [43,[90][91][92].We obtain and as discussed in Appendix D. By replacing the expression (22) in Eq. ( 21), we finally get the equation describing the dynamical cross correlation between active force and velocity lim t→∞ ) By dividing by the temperature T and taking the real part of Eq. ( 24), we obtain the equation to determine σ p (ω), that reads Equation ( 25) is not closed and cannot be used to determine the entropy production rate because contains the dependence on the unknown correlation û0 (ω) f a p (−ω) .To proceed we employ a perturbative method and expand the solution in powers of λ = δm/m ≪ 1: where the superscript (n) denotes the order of the perturbative solution that is consistent with a perturbative solution of the spectral entropy production: and of its integral over ω, i.e. ṡ, so that ṡp = ṡ(0) p + λ ṡ(1) p + λ 2 ṡ(2) p + ... .
Note that the zeroth-order entropy production is independent of the particle index p, which can be dropped at this order.Instead, the EPR is spatially dependent starting from the first correction in δm/m.In our discrete formalism, this implies that ṡ(0) p = ṡ(0) and σ (0) p (ω) = σ (0) (ω) for every p, while, in general, we expect a dependence on p, or in other words a spatial dependence, only starting from the first correction in δm/m.

B. Zeroth-order result: the homogeneous solid
The zero-order solution of Eq. ( 25) (obtained for δm = 0) corresponds to the spectral entropy, σ (0) (ω), of the homogeneous solid, in the absence of the impurity.
Setting δm = 0 in Eq. ( 25), we obtain the zeroth-order value of the spectral entropy production per particle which is independent of the position.By integrating over ω, we get the zeroth-order expression for the entropy production rate The second equality introduces the propagator, G 00 , whose expression is: The ω-integration in Eq.(E1) leading to Eq.( 31) is reported in appendix E (see Eq.(E1)).In practice, we convert the sum over wave-vectors into an integral over the Brillouin zone of volume Ω by replacing and G 00 → G(0).The prefactor in Eq.( 30) is identified with the EPR of the non-interacting system, ṡfree , according to the relation: and rewrite: The quantity ṡfree is a function of the ratio between the active temperature, v 2 0 γτ , and the thermal temperature.At fixed active temperature, it decreases as the persistence time, τ , and inertial time, τ I increase .Instead, the term G 00 accounts for the interactions among the particles in the solid, is 1 in the non-interacting limit, and G 00 ≤ 1 for the interacting case since ω 2 (q) ≥ 0. In other words, the interaction decreases the value of the EPR with respect to the free-particles case because the interactions characterizing the solid hinder the particle's ability to move with the same speed as free particles so that |v| ≪ v 0 .As a consequence, the entropy production of a solid formed by N active particles is always smaller than the entropy production of N potential-free active particles, ṡ ≤ ṡfree .

Explicit evaluation of the zeroth-order correction
To obtain an explicit expression for the zeroth order EPR of Eq.( 30) we, now, compute analytically G(0) whose value depends on the dimension of the system and the lattice structure.In the case of a triangular lattice, which is the structure found in our simulations at highdensity, the dispersion relation reads: The summand appearing in Eq.( 31) when |q| → 0 becomes of the Ornstein-Zernike form ∝ (1 + ξ 2 q 2 ) −1 and thus we can define a correlation length ξ from the relation: This length coincides with the correlation length of the spatial velocity correlation of an active solid, v(r) • v(0) and, as already discussed in Ref. [23], is an increasing function of τ and of τ I .The integral can be computed exactly as shown in the appendix B where we find where the parameters c, z and k are also given in appendix and K[k] is the complete elliptic integral of the first kind [93].In Fig. 2, the theoretical EPR, ṡ(0) , is compared with the one obtained in numerical simulations of the solid phase.Results are plotted as a function of the rescaled inertial time τ I /τ : ṡ(0) increases from zero, attains a maximum value before vanishing for large values of τ I /τ , a behavior consistent with the Clausius inequality, ṡ ≥ 0, being ṡ = 0 at equilibrium.In fact, when the persistence time, τ , is the shortest time scale of the system, (τ I /τ → ∞), the active force f a p can be assimilated to a Brownian process (persistence time ≈ 0) whose EPR is null: it is easy to verify from Eqs.( 31) and ( 33) that when ξ → 0 also ṡ → 0 because the factor ṡfree vanishes.This situation corresponds to an underdamped colloidal solid under equilibrium conditions, described τ I / τ EPR of a perfect two-dimensional crystal with triangular structure, ṡ, rescaled by the inertial time, τI = m/γ, as a function of the reduced inertial time τI /τ ).Points are obtained by numerical simulations obtained by integrating the dynamics (1) in the absence of defects, i.e. mp = m for every i, while the solid black line is calculated from the theoretical prediction (33).The parameters of the simulations are N = 10 4 and φ = 1.1.
by standard Boltzmann statistics.By decreasing τ I /τ (i.e.increasing the persistence time), the system departs from equilibrium and ṡ increases.For small deviations from equilibrium, the growth of ṡ is essentially determined by the factor ṡfree because G(0) remains close to 1 and the EPR scales as ṡ ≈ ṡfree ∼ τ /(τ + τ I ) ≈ τ /τ I for τ I /τ ≫ 1.Such a linear increase continues up to values of τ I /τ where G(0) sensibly departs from 1.This situation occurs when the size of coherent domains of the velocities becomes relevant, as revealed by the increasing velocity correlations.Indeed, when the correlation length of these domains reaches the size of the particle diameter, ξ ≈ σ, the entropy production rate attains its maximum.A further decrease of τ I /τ leads to the decrease of ṡ because the most relevant contribution to the EPR stems from the boundaries of the domains, whereas the particles (whose velocities are more aligned) located in their interior provide less relevant contributions.As a consequence, ṡ → 0 when τ I /τ → 0, the limit where ξ → ∞.We remark that this decrease is in agreement with the presence of arrested states observed numerically in systems of dense active particles in the infinite persistence time limit [94,95] for which ṡ ≈ 0.

C. First-order result: the effect of an impurity
Imperfections always break the discrete spatial translational symmetry of the crystal.Hence, some observables, including the local entropy production rate, are expected to become spatially dependent on the distance from the defect and take a constant value away from it, the one characterizing the perfect crystal.
After replacing the expression for G 0p (ω) and G p0 (ω) using Eq.( 19) and integrating over the frequency ω, we obtain the first-order perturbative correction to the EPR.Since the calculations are rather lengthy, here, we show only the final results while details of the derivation are reported in Appendix C: where and G p0 = G * 0p .By approximating the sum q by an integral, The above integral depends on the distance from the impurity, r and converges to the value G(0) as r → 0. Thus, the first-order correction to the entropy production, at distance r from the impurity, becomes: .
According to Eq. ( 42) the first order perturbative correction to the EPR at the position of the impurity, r = 0, is We estimated numerically ṡ(1) (0) from the difference ( ṡ(0) − ṡ(0) ), i.e. by subtracting the zeroth-order EPR from the total EPR obtained from the simulations.The resulting value of ṡ(1) (0) together with the zeroth-order s (1) (0) (0) τ . Entropy production rate, ṡ(0) (rescaled by τI = m/γ) calculated at the position of the defect, as a function of the reduced inertial time, τI /τ .The yellow curve denotes the zeroth-order prediction ṡ(0) i.e. the bulk value of the EPR theoretically predicted (see Eq. ( 33)).The blue light curve and the light blue points are the first-order correction of the entropy production, ṡ(1) (0), calculated at the same position.The solid line represents the theoretical predicition (43), while points are obtained from numerical simulations conducted by integrating the dynamics (1) with |δm|/m = 0.2.The values of the blue points, ṡ(1) (0), are estimated from the difference ṡ(0) − ṡ(0) between the theoretical bulk value ṡ(0) and the numerical value of ṡ(0).The parameters of the simulations were fixed at N = 10 4 and φ = 1.1.
In the large persistence regime, τ I /τ ≤ 10, the firstorder correction ṡ(1) (0) displays a decrease similar to ṡ(0) .From Eq. (43) we see that = G(0) τI τ +τI ≤ 1 .The inequality holds because both G(0) and τ I /(τ +τ I ) are less or equal to 1.Moreover, is a decreasing function of τ and an increasing function of τ I .In the limit of infinite τ , the system approaches an arrested state (with decreasing speed, |v|), and, as a consequence, not only the entropy production of the bulk ṡ(0) decreases but also the EPR due to the imperfection at r = 0 does.

Spatial profile of the entropy production generated by the impurity
To obtain explicitly the spatial profile of the entropy production, we consider G(r).Since we could not find an exact expression for it, we evaluated the integral (41) by approximating τ 2 τ I /(τ + τ I )ω 2 (q) ≈ ξ 2 q 2 and extending the integration to the whole reciprocal space.Within this approximation, the resulting formula reads: where K 0 (r/ξ) is the zeroth-order modified Bessel function of the second kind [93] and the length ξ is given by Eq. (35).The divergence at r = 0 is the result of the absence of an upper cutoff in the q-integral caused by the replacement of the Brillouin zone with an infinite integration domain.Since this problem can be easily fixed by recalling the exact evaluation of G(0), Eq. ( 36), we have an estimate of the long distance (r ≫ ξ) exponential behavior of G(r).By using the approximation (44), the first order correction to the profile of ṡ(r) is given by It displays an exponential-like decay towards its bulk value 0, since the EPR of the homogeneous crystal is just ṡ(0) .The spatial profile of ṡ(1) (r) is shown in Fig. 4 for two values of the reduced inertial time τ I /τ and reveals a good agreement between the prediction (45) and simulations where ṡ(1) (r) was estimated from the difference ( ṡ(r) − ṡ(0) ), which is exact except for terms order (δm/m) 2 .The spatial profile of ṡ(r) indicates that ξ/2 is the typical distance beyond which the entropy production is unchanged by the presence of the impurity: therefore, for large values of ξ an impurity affects the entropy production of the crystal even at long distances from its position.We recall that ξ increases as τ for τ /τ I ≪ 1 (in the underdamped regime), and as ξ ∼ √ τ for τ /τ I ≫ 1 (in the overdamped regime), while in general the value of ξ decreases as the inertia is increased (with the growth of τ I ).Interestingly, ξ/2 coincides with half of the correlation length of the spatial velocity correlations characterizing active crystals.The correlation between the velocities of two distant particles means a transfer of information that is reflected in the value of the entropy production.
Let us remark that the sign of the first-order correction to ṡ(r) depends on the sign of δm: a defect of the crystal with a mass larger than the one of the remaining particles (heavy impurity) decreases the total entropy production of the system, while an impurity with smaller mass (light impurity) increases the total entropy production.This phenomenon occurs because a light impurity means more freedom to move for the particles close to the imperfection, whereas a heavy impurity reduces the amplitude of their displacements.Regarding higher order terms in the perturbative expansion, it is possible to carry on the calculation using the same method as reported in appendix E 3 where we have derived a theoretical formula valid for two imperfections up to order (δm/m) 2 .However, its application to the present problem is complicated because one would need higher statistics in the numerical simulations.This task remains to be carried out in future work.

IV. CONCLUSIONS
As shown in the recent literature the EPR is an important theoretical tool to characterize the physics of active and living systems which often involve the presence of many degrees of freedom.However, in the case of nonuniform systems such a tool becomes much more accurate if instead of employing a global description in terms of the total EPR we use its local version, the space dependent EPR.In general, this goal is complicated to achieve because of the high dimensionality of the state space.In this paper, we have analytically and numerically investigated the space dependent EPR in the case of a active solid with point imperfections These lattice defects break the discrete translational symmetry of the crystal and induce a spatial dependence of the observable quantities.This work provides an explicit representation of the spatial variation of the EPR in a very simple model of a nonuniform system.Our analytical results for the entropy production in the non-equilibrium steady-state of an active harmonic solid, derived utilizing a perturbative expansion in powers of the mass of the impurity, agree fairly well with the numerical data obtained by computer simulation of the active solid with a soft repulsive shift and cut WCA interparticle potential and ABP dynamics.The non-uniform profile of the EPR has been calculated as a function of the distance from the impurity position.The zeroth-order result of the perturbation theory gives the value of the entropy production of a homogeneous solid [67], whose explicit expression has been reported in the case of a two-dimensional active crystal, while the first-order result accounts for the spatial dependence of the entropy production.Our study demonstrates that, in regimes of large persistence, an impurity affects the properties of the crystal also at a large distance from the impurity.Interestingly, the typical length scale that rules the spreading of information, described by the EPR, is the same length that controls the extent of the spatial velocity correlations in active crystals [30,78].Finally, we remark that the lattice imperfection is the cause of the non uniform EPR, but would not cause any velocity correlation profile in the case where only thermal noise would be present.This is in agreement with the argument that v p • v j = 2T δ ij for any equilibrium system.
Future work should extend our approach to binary and polydisperse crystals, to active crystals in three dimensions and to different kind of defects (dislocations, disclinations, grain boundaries) [96].and allows to find the explicit solution for f a n (ω) as f a n (ω) = v 0 γ √ 2τ 1 + iωτ ζn (ω) . (D5) By multiplying by f a m (−ω), taking the average over the realizations of the noise, dividing by t, and applying the limit t → ∞, we obtain Eq. (22).

Appendix E: Integrals over the frequency domain
This appendix contains some details about the derivation of the formula for first order correction to the EPR.We first compute the following integral over frequencies which appears in Eq.(??) 1. Zeroth order EPR Using the result (E1) we evaluate the zeroth order entropy production: T ṡ(0) = v 2 0 τ γ We arrive at T ṡ(0) = v 2 0 τ γ

First order EPR
According to the second equality in Eq.( 38) We now consider the first order correction.We need to perform the following frequency integral (E2) We use the following identity Gathering all together, we find the following nice formula for the EPR: The propagators G aa and G bb are equal but the prefactors depend on the signs and amplitudes of the mass ratios.Switching to continuous notation G ab → G(r ab ).It decays exponentially with the separation r ab between two lattice points and a characteristic length ξ.Perhaps, the most interesting term is the last because describes the local value of the EPR (i.e. at r p ) due to the combined effect of two imperfections sitting at r a and r b , respectively.

Figure 1 .
Figure 1.A typical solid-like snapshot of the ABP's.The red-yellow particle has a mass m+δm, whereas the remaining grey-black particles have mass m.The different colors of each half-disk represent the instantaneous orientation, along the normal to the diameter, of the active force acting on each particle.The hexagons have been drawn to emphasize the presence of the triangular lattice.
Figure 2.EPR of a perfect two-dimensional crystal with triangular structure, ṡ, rescaled by the inertial time, τI = m/γ, as a function of the reduced inertial time τI /τ ).Points are obtained by numerical simulations obtained by integrating the dynamics (1) in the absence of defects, i.e. mp = m for every i, while the solid black line is calculated from the theoretical prediction(33).The parameters of the simulations are N = 10 4 and φ = 1.1.